eLife Assessment

This manuscript describes the identification and characterization of 12 specific phosphomimetic mutations in the recombinant full-length human tau protein that trigger tau to form fibrils. This fundamental study will allow in vitro mechanistic investigations. The presented evidence is convincing. This manuscript will be of interest to all scientists in the amyloid formation field.

Reviewer #1 (Public review):

Summary and Strengths:

The very well-written manuscript by Lövestam et al. from the Scheres/Goedert groups entitled "Twelve phosphomimetic mutations induce the assembly of recombinant full-length human tau into paired helical filaments" demonstrates the in vitro production of the so-called paired helical filament Alzheimer's disease (AD) polymorph fold of tau amyloids through the introduction of 12 point mutations that attempt to mimic the disease-associated hyper-phosphorylation of tau. The presented work is very important because it enables disease-related scientific work, including seeded amyloid replication in cells, to be performed in vitro using recombinant-expressed tau protein.

Comments on revised version:

The manuscript is significantly improved, as also indicated by Reviewer 2, with the 100% formation of the PHF and the additional experiments to elucidate on the potential mechanism by the PTMs. This is a great work.

Reviewer #2 (Public review):

Summary:

This manuscript addresses an important impediment in the field of Alzheimer's disease (AD) and tauapathy research by showing that 12 specific phosphomimetic mutations in full-length tau allow the protein to aggregate into fibrils with the AD fold and the fold of chronic traumatic encephalopathy fibrils in vitro. The paper presents comprehensive structural and cell based seeding data indicating the improvement of their approach over previous in vitro attempts on non-full-length tau constructs. The main weaknesses of this work results from the fact that only up to 70% of the tau fibrils form the desired fibril polymorphs. In addition, some of the figures are of low quality and confusing.

Strengths:

This study provides significant progress towards a very important and timely topic in the amyloid community, namely the in vitro production of tau fibrils found in patients.

The 12 specific phosphomimetic mutations presented in this work will have an immediate impact in the field since they can be easily reproduced.

Multiple high-resolution structures support the success of the phosphomimetic mutation approach.

Additional data show the seeding efficiency of the resulting fibrils, their reduced tendency to bundle, and their ability to be labeled without affecting core structure or seeding capability.

Comments on revised version:

Generally, I am satisfied with the revisions. Specifically, the new results showing 100% formation of PHF is a significant improvement.

Author response:

The following is the authors’ response to the previous reviews

Reviewer #1 (Public review): 

Summary and Strengths:

The very well-written manuscript by Lövestam et al. from the Scheres/Goedert groups entitled "Twelve phosphomimetic mutations induce the assembly of recombinant fulllength human tau into paired helical filaments" demonstrates the in vitro production of the so-called paired helical filament Alzheimer's disease (AD) polymorph fold of tau amyloids through the introduction of 12 point mutations that attempt to mimic the disease-associated hyper-phosphorylation of tau. The presented work is very important because it enables disease-related scientific work, including seeded amyloid replication in cells, to be performed in vitro using recombinant-expressed tau protein. 

Weaknesses: 

The following points are asked to be addressed by the authors:

(i) In the discussion it would be helpful to note the findings that in AD the chemical structure tau (including phosphorylation) is what defines the polymorph fold and not the buffer/cellular environment. It would be further interesting to discuss these findings in respect to the relationship between disease and structure. The presented findings suggest that due to a cellular/organismal alteration, such as aging or Abeta aggregation, tau is specifically hyper-phosphorylated which then leads to its aggregation into the paired helical filaments that are associated with AD. 

We have added an extra sentence to the Introduction to emphasise this possibility: “Besides the cellular environment in which they assemble, different tau folds may also be determined by chemical modifications of tau itself.”

In addition, the last paragraph of the Discussion now reads: “It could be that, besides different cellular environments in which the filaments assemble, different posttranslational modification patterns are also important for the assembly of tau into protofilament folds that are specific for the other tauopathies.”

(ii) The conditions used for each assembly reaction are a bit hard to keep track of and somewhat ambiguous. In order to help the reader, I would suggest making a table to show conditions used for each type of assembly (including the diameter / throw of the orbital shaker) and the results (structural/biological) of those conditions. For example, presumably the authors did not have ThT in the samples used for cryo-EM but the methods section does not specify this. Also, the presence of trace NaCl is proposed as a possible cause for the CTE fold to appear in the 0N4R sample (page 4) but no explanation of why this particular sample would have more NaCl than the others. Furthermore, it appears that NaCl was actually used in the seeded assembly reactions that produced the PHF and not the CTE fold. This would seem to indicate the CTE structure of 0N4RPAD12 is not actually induced by NaCl (like it was for tau297-391). In order for the reader to better understand the reproducibility of the polymorphs, it would be helpful to indicate in how many different conditions and how many replicates with new protein preparations each polymorph was observed (could be included in the same table)  

We have added a new table (Table 1) with the buffer conditions, protein concentration and shaking speed and time, for all structures described in this paper. We never added ThT to assembly reactions that were used for cryo-EM.

We did not use NaCl in the seeded assembly reactions (we used sodium citrate). We don’t really know why 0N4R PAD12 tau more readily forms the CTE fold. The observation that it does so prompted us to use 0N3R for all ensuing experiments. 

(iii) It is not clear how the authors calculate the percentage of each filament type. In Figure 1 it is stated "discarded solved particles (coloured) and discarded filaments in grey" which leaves the reviewer wondering what a "discarded solved particle" is and which filaments were discarded. From the main text one guesses that the latter is probably false positives from automated picking but if so, these should not be referred to as filaments. Also, are the percentages calculated for filaments or segments? In any case, it would be more helpful in such are report to know the best estimate of the ratio of identified filament types without confusing the reader with a measure of the quality of the picking algorithm. Please clarify. Also, a clarification is asked for the significance of the varying degrees of PHF and AD monomer filaments in the various assembly conditions. It could be expected that there is significant variability from sample to sample but it would be interesting to know if there has been any attempt to reproduce the samples to measure this variability. If not, it might be worth mentioning so that the % values are taking with the appropriate sized grain of salt. Finally, the representation of the data in Figure 1 would seem to imply that the 0N3R forms less or no monofilament AD fold because no cross-section is shown for this structure, however it is very similar to (or statistically the same as) the 1:1 mix of 0N3R:0N4R.

In the revised manuscript, we have used bi-hierchical clustering of filaments, where each segment (or particle) is classified based on both 2D class assignment and to which filament it belongs (this method is based on [Porthula et al (2019), Ultramicroscopy 203, 132-138] and was further developed in [Lövestam et al (2024) Nature 7993, 119-125]. Based on the assumption that filament type does not change within a single filament type, we have observed that this gives excellent classification results, and that this approach allows classification of many, even small minority, filament types. Using this approach, we now quantify the different filament types on the number of segments extracted from filaments classified in this way. 

Moreover, we have also addressed the problem of having singlets among the PHF preparation: it turns out that waiting longer, just by transferring samples out of the shaker after one week and incubating it quiescently at 37 ºC for two more weeks, the singlets disappear and only PHFs remain. Filaments made for the fluorophore labelling in the revised Figure 3 were also done using the new protocol. In total, we have N=7 replicates with a mean of 95.3% PHFs and a standard deviation of 9.4%. The revised text in the Results section reads:

“To further increase the proportions of PHFs-to-singlet ratio, we removed the plate from the shaker after one week and incubated it quiescently at 37 ºC for two more weeks. This resulted in 100% PHFs formed (Figure 1 – figure supplement 4). When repeated seven times, on average 95.3% PHFs formed, with 25% of singlets formed in a single outlier (Figure 1 – figure supplement 5)” 

(iv) The interpretation of the NMR data on soluble tau that the mutations on the second site are suppressing in part long range dynamic interaction around the aggregationinitiation site (FIA) is sound. It is in particular interesting to find that the mutations have a similar effect as the truncation at residue 391. An additional experiment using solvent PREs to elaborate on the solvent exposed sequence-resolved electrostatic potential and the intra-molecular long range interactions would likely strengthen the interpretation significantly (Iwahara, for example, Yu et al, in JACS 2024). Figure 6D Figure supplement shows the NMR cross peak intensities between tau 151-391 and PAD12tau151-391. Overall the intensities of the PAD12 tau construct are more intense which could be interpreted with less conformational exchange between long range dynamic interactions. There are however several regions which do not show any intensity anymore when compared with the corresponding wildtype construct such as 259-262, 292-294 which should be discussed/explained. 

While long-range intramolecular interactions of tau have previously been reported through the use of spin labels (Mukrasch et al 2009 PLoS Biol 7(2): e1000034), we have been hesitant to introduce paramagnetic agents into our samples for two reasons. First, the bulky size of the spin label may affect filament formation or influence the dynamic properties of the protein. Second, covalent addition of the spin label requires mutation of the primary sequence to both remove native cysteine residues and add cysteines at the desired label location. We have previously shown that mutation of cysteine 322 to alanine leads to the formation of tau filaments with a structure that is different from the PHF (Santambrogio et al (2025) bioRxiv 2025.03.29.646137). 

Instead, we have included in the revised manuscript new NMR and cryo-EM data that provide further support for the model that a FIA-like interaction between residues ₃₉₂IVYK₃₉₅ and residues ₃₀₆VQIVYK₃₁₁ has an inhibiting effect on filament nucleation in unmodified full-length tau. A mutant of tau297-441 where residues ₃₉₂IVYK₃₉₅ have been deleted and that does not contain the four PAD12 mutations in the carboxy-terminal domain behaves similarly in the NMR experiment as the tau297-441 construct with those four PAD12 mutations. Moreover, full-length 0N3R tau with the eight PAD12 mutations in the amino-terminal fuzzy coat and with the deletion of₃₉₂IVYK₃₉₅, but without the four PAD12 mutations in the carboxy-terminal domain, assembles readily into amyloid filaments (of which we also solved a cryo-EM structure, see the revised Figure 6B). These observations provide mechanistic insights into the previously proposed paper-clip model [Jeganathan (2008), J Biol Chem 283, 32066-32076], where interactions between the fuzzy coat inhibit filament formation of unmodified full-length tau, and phosphorylation in the fuzzy coat interferes with these interactions, thus leading to filament nucleation. Of course, the identification of residues ₃₉₂IVYK₃₉₅ for this interaction also explain why truncation of tau at residue 391 leads to spontaneous assembly. We have introduced a new Figure 7 to the revised manuscript to explain this model in more detail. The corresponding new section in the Results reads:

“To investigate this further, we also tested a tau construct comprising residues tau297-441 without the phosphomimetic mutations, but with a deletion of residues (Δ392-395). Filaments formed rapidly and the cryo-EM structure showed that the ordered core consisted of the amino-terminal part of the construct spanning residues 297-318 (Figure 6B). NMR analysis (Figure 6 – figure supplement 5B) showed that the tau297441 Δ392-395 construct exhibited similar backbone rigidity properties to the tau297-441 PAD12 construct, despite peak locations and local secondary structural propensities being more similar to the wildtype tau297-441 (Figure 6 – figure supplement 5A; Figure 6 – figure supplement 6). HSQC peak intensities in the 297-319 and 392-404 regions of tau297-441 Δ392-395 (Figure 6A, expanded from Figure 6 - figure supplement 5C) were like those in the tau297-441 PAD12. These data suggest that the IVYK deletion has a similar effect as the phosphomimetics on residues 396, 400, 403 and 404 on disrupting an intra-molecular interaction between the FIA core region and the carboxy-terminal domain, which may therefore be mediated by interactions between the two IVYK motifs that are similar to those observed in the FIA (Lövestam et al, 2024).”

A new section in the Discussion now reads:

“Our NMR data provide insights into the mechanism by which phosphorylation in the fuzzy coat of tau, or truncations of tau, lead to the formation of filaments with ordered cores of residues that are themselves not phosphorylated. HSQC peak intensity differences between unmodified tau 297-441, PAD12 tau 297-441 and tau297-391 suggest that phosphorylation of the fuzzy coat, particularly near the ₃₉₂IVYK₃₉₅ motif in the carboxy-terminal domain, a7ects the conformation of the residues of tau that become ordered in the FIA (Lövestam et al., 2024). Removal of residues ₃₉₂IVYK₃₉₅ in the carboxyterminal domain of tau 297-441 led to rapid filament formation in the absence of phosphomimetics, while HSQC peak intensity di7erences for this construct indicate similar backbone rigidity compared to tau 297-441 without the deletion, but with the four PAD12 mutations in the carboxy-terminal domain. Combined, these observations support a model where the ₃₉₂IVYK₃₉₅ motif in unmodified full-length tau monomers interacts with the ₃₀₈IVYK₃₁₁ motif, thus inhibiting filament formation by preventing the formation of the nucleating species, the FIA. Phosphorylation of nearby residues 396, 400, 403 and 404, or truncation at residue 391, disrupt this interaction and lead to filament formation. This model agrees with the previously proposed hairpin-like model of tau (Jeganathan et al., 2008), although the corresponding interaction between the aminoterminal domain of tau and the core-forming region remains unknown (Figure 7).”

Due to the challenging nature of the assignment, it was not possible to assign all residues in the HSQC of the tau151-391 and the PAD12 tau151-391 samples, including residues 259-262 and 292-294 for PAD12 tau151-391. To make this clearer, we have marked residues that are not assigned with an asterisk in the revised version of Figure 6 – figure supplement 1.  

(v) Concerning the Cryo-EM data from the different hyper-phosphorylation mimics, it would seem that the authors could at least comment on the proportion of monofilament and paired-filaments even if they could not solve the structures. Nonetheless, based on their previous publications, one would also expect that they could show whether the nontwisted filaments are likely to have the same structure (by comparing the 2D classes to projections of non-twisted models). Also, it is very interesting to note that the twist could be so strongly controlled by the charge distribution on the non-structured regions (and may be also related to the work by Mezzenga on twist rate and buffer conditions). Is the result reported in Figure 2 a one-oT case or was it also reproducible?

As also indicated in the main text, the assembly conditions for the PAD12+4, PAD12-4 and PAD12+/-4 constructs were kept the same as those for the PAD12 construct. It is possible that further optimisation of the conditions could again lead to twisting filaments, but we chose not to pursue this route. With unlimited resources and time, one could assess in detail which of the PAD12 mutations are required and which ones could be omitted to form PHFs. However, this would require a lot of work and cryo-EM time. For now, we chose to prioritise reporting conditions that do work to reproducibly make PHFs in the laboratory (using the PAD12 construct) and leave the more detailed analysis of other constructs for future studies. 

Reviewer #2 (Public review): 

Summary: 

This manuscript addresses an important impediment in the field of Alzheimer's disease (AD) and tauapathy research by showing that 12 specific phosphomimetic mutations in full-length tau allow the protein to aggregate into fibrils with the AD fold and the fold of chronic traumatic encephalopathy fibrils in vitro. The paper presents comprehensive structural and cell based seeding data indicating the improvement of their approach over previous in vitro attempts on non-full-length tau constructs. The main weaknesses of this work results from the fact that only up to 70% of the tau fibrils form the desired fibril polymorphs. In addition, some of the figures are of low quality and confusing. 

As also explained in our response to reviewer #1, we have performed better quantification of filament types in the revised manuscript, and we have investigated how to get rid of the singlets. In the revised manuscript, we report that singlets disappear as time passes and that one can obtain 100% pure PHFs by quiescently incubating samples for another two weeks, after shaking for a week.

Strengths: 

This study provides significant progress towards a very important and timely topic in the amyloid community, namely the in vitro production of tau fibrils found in patients.

The 12 specific phosphomimetic mutations presented in this work will have an immediate impact in the field since they can be easily reproduced.

Multiple high-resolution structures support the success of the phosphomimetic mutation approach. Additional data show the seeding efficiency of the resulting fibrils, their reduced tendency to bundle, and their ability to be labeled without affecting core structure or seeding capability.

Weaknesses: 

Despite the success of making full-length AD tau fibrils, still ~30% of the fibrils are either not PHF, or not accounted for. A small fraction of the fibrils are single filaments and another ~20% are not accounted for. The authors mention that ~20% of these fibrils were not picked by the automated algorithm. However, it would be important to get additional clarity about these fibrils. Therefore, it would improve the impact of the paper if the authors could manually analyze passed-over particles to see if they are compatible with PHF or fall into a different class of fibrils. In addition, it would be helpful if the authors could comment on what can be done/tried to get the PHF yield closer to 90-100%

As mentioned above, in the revised manuscript we show that the singlets disappear over time and we now include a description of a method that leads to 100% PHF formation.

Reviewer #1 (Recommendations for the authors):

Minor points: 

(a) In Figure 6 the dashed purple vertical lines overlap with the black bars, rendering a grey color which is confusing because the grey bars used for the shorter construct. It is suggested to improve the colors (remove transparency on the purple?)

We thank the reviewers for their suggestions for improving the visualisation of our data. We have recoloured the tau297-391 data from grey to gold and moved the dashed lines to the back of image to remove the apparent colour changes.  

(b) Is there any support for the suggestion that "part of the second microtubule-binding repeat is ordered" being "related to this construct forming filaments with only a single protofilament"? It seemed to have come out of nowhere.

There is no further support for this statement, but we thought it would be worth hypothesizing about this observation. 

(c) Figures 1 and 4 E is better described as a "main chain trace" or "backbone trace" although the latter usually refers to only CA positions. Ribbon usually refers to something else in representations of protein structures. 

This has been changed into “main chain trace” in Figures 1 and 4. 

(d) Figure 1 Supplement 3: Panel letters in the legend do not match. 

This has been fixed.

Reviewer #2 (Recommendations for the authors): 

The introduction is a bit lengthy (e.g. 3rd paragraph of introduction) and could benefit by focusing specific question the manuscript addresses. 

We have shortened the Introduction. It now contains ~1150 words, which we hope provides a better compromise between length and sufficient background information.

Figure captions are generally not helpful in conveying a message to the reader.

Figure 1 - figure supplement 3 is quite confusing. The 4 structures in A) do not correspond to the grids in B-E. What is this figure supposed to show?

This confusion was probably the result of incorrect labelling of panels in the legend, which was also pointed out by reviewer #1. This has been fixed in the revised manuscript.

Page 11: Although I know what you mean, 'linear increase of ThT fluorescence' is not the correct term. 

We have replaced “linear” with “rapid”.

Page 15: Although line shape and peak intensity can be related you are not reporting on line shape or width but simply on peak intensity. Therefore, I wouldn't talk about the result of a 'line shape analysis'.

We have changed the wording accordingly. 

Figure 6 (and supplement 1) are confusing and too small to be readable in print. It might be sufficient to show the CSP and upload the remaining data to the BMRB. 

We have made a clearer version of the main NMR Figure 6 in the revised manuscript showing the most pertinent NMR data and have moved the previous version into the figure supplements. We designed these figures to be viewed as full page A4 panels, ideally seen in one image as they show multiple comparisons of different experiments and constructs.

As such we feel these will be best viewed on screen as part of the eLife web document. We have uploaded HSQC spectra and assignments to the BMRB (see below).

Figure 6 supplement 3 might benefit from pointing out key residues in the overlay.

We have added the labels (this is now Figure 6 supplement 4).

Data availability: Please upload the assignments to the BMRB together with key spectra (e.g. HSQCs). 

We have uploaded HSQC data along with our assignments to the BMRB, the accession codes are 52694 – tau297-441 wt; 52695 – tau297-441 PAD-12; 52696 – tau151-391 wt; 52697 – tau151-391 PAD-12; and 53230 – tau297-441 delta392-395.  These accession codes have been added to the manuscript. 

The quality of some of the figures (specifically Figure 1 - supplement 3 and Figure 6) is not suitable for publication. 

For the original submission to bioRxiv, we produced a single PDF with a manageable file size. We will liaise with the eLife staff to ensure the images used in the version of record will be suitable for publication.

eLife Assessment

This important work presents a stochastic branching process model of tumour-immune coevolution, incorporating stochastic antigenic mutation accumulation and escape within the cancer cell population. They then used this model to investigate how tumour-immune interactions influence tumour outcome and the summary statistics of sequencing data of bulk and single-cell sequencing of a tumour. The evidence is compelling and the work will be of interest to cancer-immune biology fields.

Reviewer #1 (Public review):

Summary:

The topic of tumor-immune co-evolution is an important, understudied topic with, as the authors noted, a general dearth of good models in this space. The authors have made important progress on the topic by introduced a stochastic branching process model of antigenicity / immunogenicity and measuring the proportion of simulated tumors which go extinct. The model is extensively explored and authors provide some nice theoretical results in addition to simulated results, including an analysis of increasing cancer/immune versus cyclical cancer/immune dynamics. The analysis appropriately builds upon the foundation of other work in the field of predicting site frequency spectrum, but extends the results into cancer-immune co-evolution in an intuitive computational framework.

Author response:

The following is the authors’ response to the original reviews.

Reviewer #1 (Public review): 

Summary: 

The topic of tumor-immune co-evolution is an important, understudied topic with, as the authors  noted, a general dearth of good models in this space. The authors have made important progress on the topic by introducing a stochastic branching process model of antigenicity/immunogenicity and measuring the proportion of simulated tumors that go extinct. The model is extensively explored, and the authors provide some nice theoretical results in addition to simulated results. 

We thank the reviewer for the positive comments on our work.

Major comments 

The text in lines 183-191 is intuitively and nicely explained. However, I am not sure all of it follows from the figure panels in Figure 2. For example, the authors refer to a mutation that has a large immunogenicity, but it's not shown how many mutations, or the relative size of the mutations in Figure 2. The same comment holds true for the claim that spikes also arise for mutations with low antigenicity. 

We thank the reviewer for helping us to further specify this statement in our original submission. We now added muller plots in a new Appendix Figure (Figure A3) presenting the relative abundances of different types of effector cells in the population over time. Each effector type is colour-coded with its antigenicity and immunogenicity. To align with this Appendix Figure (Figure A3), we also updated our Figure 2 generated under the same realisation as Figure A3. We can now see clearly that the spikes in the mean values of the antigenicity and immunogenicity over the whole effector populations in new Figure 2B&2D indeed correspond to the expansion of single or several antigenic mutations recruiting the specific effector cell types. For example, in Figure 2B, we can see that the spikes of low average antigenicity and high immunogenicity (around time 11) happen at the same time when an effector type in Figure A3 with such a trait (coloured in green) arises and takes over the population. We have rewritten our Results section related (Line 192 - Line 222 in main text and Appendix A6).

Reviewer #2 (Public review): 

Summary: 

In this work, the authors developed a model of tumour-immune dynamics, incorporating stochastic antigenic mutation accumulation and escape within the cancer cell population. They then used this  model to investigate how tumour-immune interactions influence tumour outcome and summary  statistics of sequencing data. 

Strengths: 

This novel modeling framework addresses an important and timely topic. The authors consider the useful question of how bulk and single-cell sequencing may provide insights into the tumourimmune interactions and selection processes. 

We thank the reviewer for the positive comments.

Weaknesses: 

One set of conclusions presented in the paper is the presence of cyclic dynamics between effector/cancer cells, antigenicity, and immunogenicity. However, these conclusions are supported in the manuscript by two sample trajectories of stochastic simulations, and these provide mixed support for the conclusions (i.e. the phasing asynchrony described in the text does not seem to apply to Figure 2C). 

We have now developed a method to quantify the cyclic dynamics in our system (Appendix A7), where can track the directional changes phase portrait of the abundances of the cancer and effector cells. We first tested this method in a non-evolving stochastic predator-prey system, where our method can correctly capture the number of cycles in this system (Figure A7). We then use this method to quantify the number of cycles we observed between cancer and effector cells under different mutation rates (Figure A5) as well as whether they are counter-clockwise or clockwise cycles (Figure A6). Our results showed that the cyclic dynamics are more often to be observed when mutation rates are higher, and the majority of those cycles are counter-clockwise. When the mutation rate is high, we observe an increase of clockwise cycles, which have been observed in predator-prey systems and explained through coevolution. However, even under high mutation rates, counter-clockwise cycles are still the more frequent type. 

In our simulations, we observed rarely out-of-phase cycles, which was by chance present in our original Figure 2. We have now removed that statement about out-of-phase cycles and replaced by more systematic analysis of the cyclic dynamics as described above (Line 192 to 207 in the revised version). We thank the constructive comment of the reviewer, which motivated us to improve our analysis significantly. 

Similarly, the authors also find immune selection effects on the shape of the mutational burden in Figure 5 D/H using a qualitative comparison between the distributions and theoretical predictions in  the absence of immune response. However the discrepancy appears quite small in panel D, and  there are no quantitative comparisons provided to evaluate the significance. An analysis of the robustness of all the conclusions to parameter variation is missing. 

We have now added statistical analysis using Wasserstein distance between the simulated mutation burden distribution and theoretical (neutral) expectation in Figure 5 C, D, G, H as well as in Figure A11 C&D when there is no cancer-immune interaction. We can see that the measurements of the  Wasserstein distance agrees with our statement, that the higher immune effectiveness leads to larger deviation from the neutral expectation.

Lastly, the role of the Appendix results in the main messages of the paper is unclear. 

We agree with the review and have now removed the Appendix sections “Deterministic Analysis”. 

Reviewing Editor Comments: 

I find the abstract too long. For example, "Knowledge of this coevolutionary system and the selection taking place within it can help us understand tumour-immune dynamics both during tumorigenesis but also when treatments such as immunotherapies are applied." can be shortened to: "Knowledge of this coevolutionary system can help us understand tumour-immune dynamics both during tumorigenesis and during immunotherapy treatments." 

We agree and have taken the suggestion of the reviewer to shorten our abstract.

Reviewer #1 (Recommendations for the authors): 

The discussion at lines 134-140, centered around Figure A1, is an important and nicely constructed feature of the model. 

Reviewer #2 (Recommendations for the authors): 

I suggest that the authors conduct a more in-depth analysis of their conclusions on cyclic dynamics over a large set of sample paths.

Done and please see our detailed response to the reviewer 2 above.

In addition, statistical comparisons between the observed mutational burden distribution and  theoretical predictions in the absence of immune selection should be carried out to support their conclusions. In all cases, conclusions should be tested extensively for robustness/sensitivity to parameters. 

Done and please see our detailed response to the reviewer 2 above.

Here are some specific suggestions/comments: 

(1) Please provide a precise mathematical description of the model to complement Figure 1. 

We have significantly revised our “Model” section to provide a precise mathematical description of our model (Line 138 - 148). Please also see our document showing the difference between the revised version and original submission.

(2) Section on "Interactions dictate outcome of tumour progress" and Figure 3: please define 'tumour outcome' - are the heatmaps produced in Figure 3 tumor size reflecting whether or not the population has reached level K before a particular time? Also, I do not see a definition for the 'slowgrowing' tumour proportion plotted in Figure 3CF or in the accompanying text. 

We have now added the definition of “tumour outcome” in our “Model” section (line 171 to 176), where we explain our model parameters and quantities measured in the following “Results” section.

(3) Figure 5C/G: the green dotted vertical line is difficult to see. 

We have now changed the mean of the simulations to solid red lines instead of using the green dotted vertical lines previously.

(4) Appendix A1 text under (A2) should U/N be U/C? N does not appear to be defined. 

We have more removed the previous A1 section. Please see our response to reviewer 2 as well.

(5) Text under (A5): it is unclear what is meant by "SFS must be heavy tailed (that is, more heterogeneous)" -- a more precise statement regarding tail decay rate and associated consequences would be more helpful. 

We have more removed the previous A section, where the original text "...SFS must be heavy-tailed" was.

(6) Section A4 and Figure A1: can these calculations be compared to simulations? 

We have more removed the previous A section on the deterministic analysis as they are not so  relevant to our stochastic simulations indeed. Please see our response to reviewer 2 as well.

(7) Also, in general, please clarify how the results in the Appendix are used in the main text conclusions or provide insights relevant to these conclusions. If they are not, one can consider removing them.  

We have more removed the previous A section on the deterministic analysis. The remaining sections are about stochastic simulations and extended figures which support our main figures.  

(8) Figure A2: the two lines are difficult to tell apart on each panel. Please consider different styles.

We have changed one of the dotted lines to be solid. This figure is now Figure A1 in our revision.

eLife Assessment

This important study introduces a new class of spectrally tunable, dye-based calcium sensors optimized for imaging in organelles with high calcium concentrations, such as the endoplasmic reticulum and mitochondria. The experimental evidence supporting the applicability of these sensors is convincing, with thorough validation in cultured cells and neurons. The work will be of high interest to researchers studying calcium signaling dynamics in subcellular compartments.

Reviewer #1 (Public review):

Summary:

The manuscript by Moret et al. details the development and characterisation of novel ER- and mitochondria-targeted genetically encoded chemogenic Ca2+ sensors.

Strengths:

Compared to existing probes, these sensors exhibited superior responsiveness, brightness, and photostability within the red and far-red emission spectrum, enabling triple compartment Ca2+ measurements (ER, mitochondria, cytosol) and the detection of Ca2+ dynamics in axons and dendrites.

Weaknesses:

The data are robust and convincing, although the manuscript text lacks precision.

Reviewer #2 (Public review):

Summary:

Moret et al. present an engineered family of fluorescent calcium indicators based on HaloCamp, a HaloTag-based sensor system that utilizes Janelia Fluorophores (JF dyes) to report calcium dynamics. By introducing single or multiple amino acid substitutions, the authors reduce HaloCamp's calcium affinity, making these low-affinity variants well-suited for imaging calcium transients in high-calcium environments such as the endoplasmic reticulum (ER) and mitochondria. The study validates the sensors' dissociation constants (Kd), spectra, and multiplex capabilities. It demonstrates improved performance compared to existing tools when targeted to subcellular compartments in mammalian cells and cultured neurons. The sensors can be tuned across the red-to-far-red spectrum via JF585 and JF635 labeling, enabling flexible multiplexed imaging. For example, the authors show that HaloCamp can be targeted to mitochondria and used alongside other green and red sensors, allowing simultaneous imaging of calcium dynamics in the cytosol, ER, and mitochondria. Overall, they achieve their goals, and the data demonstrate that HaloCamp variants are effective for detecting ER and mitochondrial calcium changes under physiological conditions. The presented experiments support the conclusions. However, some key aspects, such as sensor kinetics and axonal validation, would benefit from further analysis.

This work is likely to have an important impact on the fields of calcium imaging and organelle physiology. The modular design of HaloCamp and its compatibility with a wide range of fluorophores offer a broad application range for cell biologists and neuroscientists.

Strengths:

(1) The authors introduce the first tunable, dye-based, low-affinity HaloTag calcium sensors for subcellular imaging, addressing a significant unmet need for ER and mitochondrial calcium detection.

(2) The ability to pair HaloCamp with JF585 and JF635 extends the spectral range, facilitating multiplexed imaging with existing calcium indicators.

(3) The sensors are validated in a range of subcellular compartments (ER, mitochondria, cytosol) in both mammalian cells and neurons.

(4) The authors successfully demonstrate simultaneous imaging of three compartments using orthogonal sensors, a technically impressive feat.

(5) Kd values are measured, and fluorescent responses are tested under physiologically relevant stimulation.

Weaknesses:

(1) The authors do not quantify the kinetics (e.g., decay tau or off-rate) of the fluorescent signals, particularly after stimulation. For example, in the ER imaging experiments in neurons, the decay of the HaloCamp fluorescence after field stimulation (20 APs @ 20 Hz) is not analyzed or compared to ER-GCaMP6-210 or R-CEPIer.

(2) It remains unclear whether the observed decay represents the sensor's off-kinetics or actual physiological calcium clearance from the ER. A comparison between sensors or an independent measurement of ER clearance rates in vitro would clarify this.

(3) The choice of 20 APs at 20 Hz is not justified. Specifically, single APs or low-frequency stimulations are not tested, leaving unclear what the detection threshold of the new sensors is.

(4) In neuron experiments, the authors report measuring ER calcium in axons based presumably on morphology, but no specific justification for selection, markers, or post hoc labeling is described.

(5) Figure 5 assumes that all three indicators (cytosolic, ER, and mitochondrial) are fast enough to report calcium dynamics in response to histamine. This assumption is not fully validated. Cross-controls (e.g., expressing GCaMP6-210 in mitochondria and HaloCamp in the ER) would strengthen confidence that the sensors are correctly reporting dynamic changes.

(6) It is not clear why Thapsigargin leads to depletion in HeLa cells and neurons in experiments shown in Figure 1E, but not in 2B upon field stimulation.

eLife Assessment

This study presents useful findings on the molecular mechanisms driving female-to-male sex reversal in the ricefield eel (Monopterus albus) during aging, which would be of interest to biologists studying sex determination. The manuscript describes an interesting mechanism potentially underlying sex differentiation in M. albus. However, the current data are incomplete and would benefit from more rigorous experimental approaches.

Reviewer #1 (Public review):

Summary:

This study investigates the molecular mechanism by which warm temperature induces female-to-male sex reversal in the ricefield eel (Monopterus albus), a protogynous hermaphroditic fish of significant aquacultural value in China. The study identifies Trpv4 - a temperature-sensitive Ca²⁺ channel - as a putative thermosensor linking environmental temperature to sex determination. The authors propose that Trpv4 causes Ca²⁺ influx, leading to activation of Stat3 (pStat3). pStat3 then transcriptionally upregulates the histone demethylase Kdm6b (aka Jmjd3), leading to increased dmrt1 gene expression and ovo-testes development. This work aims to bridge ecological cues with molecular and epigenetic regulators of sex change and has potential implications for sex control in aquaculture.

Strengths:

(1) This study proposes the first mechanistic pathway linking thermal cues to natural sex reversal in adult ricefield eel, extending the temperature-dependent sex determination paradigm beyond embryonic reptiles and saltwater fish.

(2) The findings could have applications for aquaculture, where skewed sex ratios apparently limit breeding efficiency.

Weaknesses:

(A) Scientific Concerns:

(1) There is insufficient replication and data transparency. First, the qPCR data are presented as bar graphs without individual data points, making it impossible to assess variability or replication. Please show all individual data points and clarify n (sample size) per group. Second, the Western blotting is only shown as single replicates. If repeated 2-3 times as stated, quantification and normalization (e.g., pStat3/Stat3, GAPDH loading control) are essential. The full, uncropped blots should be included in the supplementary data.

(2) The biological significance of the results is not clear. Many reported fold changes (e.g., kdm6b modulation by Stat3 inhibition, sox9a in S3A) are modest (<2-fold), raising concerns about biological relevance. Can the authors define thresholds of functional relevance or confirm phenotypic outcomes in these animals?

(3) The specificity of key antibodies is not validated. Key antibodies (Stat3, pStat3, Foxl2, Amh) were raised against mammalian proteins. Their specificity for ricefield eel proteins is unverified. Validation should include siRNA-mediated knockdown with immunoblot quantification with 3 replicates. Homemade antibodies (Sox9a, Dmrt1) also require rigorous validation.

(4) Most of the imaging data (immunofluorescence) is inconclusive. Immunofluorescence panels are small and lack monochrome channels, which severely limits interpretability. Larger, better-contrasted images (showing the merge and the monochrome of important channels) and quantification would enhance the clarity of these findings.

(B) Other comments about the science:

(1) In S3A, sox9a expression is not dose-responsive to Trpv4 modulation, weakening the causal inference.

(2) An antibody against Kdm6b (if available) should be used to confirm protein-level changes.

In sum, the interpretations are limited by the above concerns regarding data presentation and reagent specificity.

Reviewer #2 (Public review):

Summary:

This study presents valuable findings on the molecular mechanisms driving the female-to-male transformation in the ricefield eel (Monopterus albus) during aging. The authors explore the role of temperature-activated TRPV4 signaling in promoting testicular differentiation, proposing a TRPV4-Ca²⁺-pSTAT3-Kdm6b axis that facilitates this gonadal shift.

Strengths:

The manuscript describes an interesting mechanism potentially underlying sex differentiation in M. albus.

Weaknesses:

The current data are insufficient to fully support the central claims, and the study would benefit from more rigorous experimental approaches.

(1) Overstated Title and Claims:

The title "TRPV4 mediates temperature-induced sex change" overstates the evidence. No histological confirmation of gonadal transformation (e.g., formation of testicular structures) is presented. Conclusions are based solely on molecular markers such as dmrt1 and sox9a, which, although suggestive, are not definitive indicators of functional sex reversal.

(2) Temperature vs Growth Rate Confounding (Figure 1E):

The conclusion that warm temperature directly induces gonadal transformation is confounded by potential growth rate effects. The authors state that body size was "comparable" between 25{degree sign}C and 33{degree sign}C groups, but fail to provide supporting data. In ectotherms, growth is intrinsically temperature-dependent. Given the known correlation between size and sex change in M. albus, growth rate-rather than temperature per se-may underlie the observed sex ratio shifts. Controlled growth-matched comparisons or inclusion of growth rate metrics are needed.

(3) TRPV4 as a Thermosensor-Insufficient Evidence:

The characterisation of TRPV4 as a direct thermosensor lacks biophysical validation. The observed transcriptional upregulation of Trpv4 under heat (Figure 2) reflects downstream responses rather than primary sensor function. Functional thermosensors, including TRPV4, respond to heat via immediate ion channel activity-typically measurable within seconds-not mRNA expression over hours. No patch-clamp or electrophysiological data are provided to confirm TRPV4 activation thresholds in eel gonadal cells. Additionally, the Ca²⁺ imaging assay (Figure 2F) lacks essential details: the timing of GSK1016790A/RN1734 administration relative to imaging is unclear, making it difficult to distinguish direct channel activity from indirect transcriptional effects.

(4) Cellular Context of TRPV4 Activity Is Unclear:

In situ hybridisation suggests TRPV4 expression shifts from interstitial to somatic domains under heat (Figures. 2H, S2C), implying potential cell-type-specific roles. However, the study does not clarify: (i) whether TRPV4 plays the same role across these cell types, (ii) why somatic cells show stronger signal amplification, or (iii) the cellular composition of explants used in in vitro assays. Without this resolution, conclusions from pharmacological manipulation (e.g., GSK1016790A effects) cannot be definitively linked to specific cell populations.

(5) Rapid Trpv4 mRNA Elevation and Channel Function:

The authors report a dramatic increase in Trpv4 mRNA within one day of heat exposure (Figures 4D, S2B). Given that TRPV4 is a membrane channel, not a transcription factor, its rapid transcriptional sensitivity to temperature raises mechanistic questions. This finding, while intriguing, seems more correlational than functional. A clearer explanation of how TRPV4 senses temperature at the molecular level is needed.

(6) Inconclusive Evidence for the Ca²⁺ -pSTAT3-Kdm6b Axis:

Although the authors propose a TRPV4-Ca²⁺ -pSTAT3-Kdm6b-dmrt1 pathway, intermediate steps remain poorly supported. For example, western blot data (Figures 3C, 4B) do not convincingly demonstrate significant pSTAT3 elevation at 34{degree sign}C. Higher-resolution and properly quantified blots are essential. The inferred signalling cascade is based largely on temporal correlation and pharmacological inhibition, which are insufficient to establish direct regulatory relationships.

(7) Species-Specific STAT3-Kdm6b Regulation Is Unresolved:

The proposed activation of Kdm6b by pSTAT3 contrasts with findings in the red-eared slider turtle (Trachemys scripta), where pSTAT3 represses Kdm6b. This divergence in regulatory direction between the two TSD species is surprising and demands further justification. Cross-species differences in binding motifs or epigenetic context should be explored. Additional evidence, such as luciferase reporter assays (using wild-type and mutant pSTAT3 binding motifs in the Kdm6b promoter) is needed to confirm direct activation. A rescue experiment-testing whether Kdm6b overexpression can compensate for pSTAT3 inhibition-would also greatly strengthen the model.

(8) Immunofluorescence-Lack of Structural Markers:

All immunofluorescence images should include structural markers to delineate gonadal boundaries. Furthermore, image descriptions in the figure legends and main text lack detail and should be significantly expanded for clarity.

(9) Pharmacological Reagents-Mechanisms and References:

The manuscript lacks proper references and mechanistic descriptions for the pharmacological agents used (e.g., GSK1016790A, RN1734, Stattic). Established literature on their specificity and usage context should be cited to support their application and interpretation in this study.

(10) Efficiency of Experimental Interventions:

The percentage of gonads exhibiting sex reversal following pharmacological or RNAi treatments should be reported in the Results. This is critical for evaluating the strength and reproducibility of the interventions.

?

eLife Assessment

This important work advances our understanding of DNA methylation and its consequences for susceptibility to DNA damage. This work presents evidence that DNA methylation can accentuate the genomic damage propagated by DNA damaging agents as well as potentially being an independent source of such damage. The experimental results reported are sound. The evidence presented to support the conclusions drawn is convincing and alternative interpretations are considered. The work will be of broad interest to biochemists, cell and genome biologists.

Reviewer #1 (Public review):

Summary:

The manuscript titled "Introduction of cytosine-5 DNA methylation sensitizes cells to oxidative damage" proposes that 5mC modifications to DNA, despite being ancient and wide-spread throughout life, represent a vulnerability, making cells more susceptible to both chemical alkylation and, of more general importance, reactive oxygen species. Sarkies et al take the innovative approach of introducing enzymatic genome-wide cytosine methylation system (DNA methyltransferases, DNMTs) into E. coli, which normally lacks such a system. They provide compelling evidence that the introduction of DNMTs increases the sensitivity of E. coli to chemical alkylation damage. Surprisingly they also show DNMTs increase the sensitivity to reactive oxygen species and propose that the DNMT generated 5mC presents a target for the reactive oxygen species that is especially damaging to cells. Evidence is presented that DNMT activity directly or indirectly produces reactive oxygen species in vivo, which is an important discovery if correct, though the mechanism for this remains obscure.

I am satisfied that the points #2, #3 and #4 relating to non-addativity, transcriptional changes and ROS generation have been appropriately addressed in this revised manuscript. The most important point (previously #1) has not been addressed beyond the acknowledgement in the results section that: "Alternatively, 3mC induction by DNMT may lead to increased levels of ssDNA, particularly in alkB mutants, which could increase the risk of further DNA damage by MMS exposure and heighten sensitivity." This slightly miss-represents the original point that 5mC the main enzymatic product of DNMTs rather or in addition to 3mC is likely to lead to transient damage susceptible ssDNA, especially in an alkB deficient background. And more centrally to the main claims of this manuscript, the authors have not resolved whether methylated cytosine introduced into bacteria is deleterious in the context of genotoxic stress because of the oxidative modification to 5mC and 3mC, or because of oxidative/chemical attack to ssDNA that is transiently exposed in the repair processing of 5mC and 3mC, especially in an alkB deficient background. This is a crucial distinction because chemical vulnerability of 5mC would likely be a universal property of cytosine methylation across life, but the wide-spread exposure of ssDNA is expected to be peculiarity of introducing cytosine methylation into a system not evolved with that modification as a standard component of its genome.

These two models make different predictions about the predominant mutation types generated, in the authors system using M.SssI that targets C in a CG context - if oxidative damage to 5mC dominates then mutations are expected to be predominantly in a CG context, if ssDNA exposure effects dominate then the mutations are expected to be more widely distributed - sequencing post exposure clones could resolve this.

Strengths:

This work is based on an interesting initial premise, it is well motivated in the introduction and the manuscript is clearly written. The results themselves are compelling.

Weaknesses:

I am not currently convinced by the principal interpretations and think that other explanations based on known phenomena could account for key results. Specifically the authors have not resolved whether oxidative modification to 5mC and 3mC, or chemical attack to ssDNA that is transiently exposed in the repair processing of 5mC and 3mC is the principal source of the observed genotoxicity. The authors acknowledge this potential alternative model in their discussion of the revised manuscript.

Reviewer #2 (Public review):

5-methylcytosine (5mC) is a key epigenetic mark in DNA and plays a crucial role in regulating gene expression in many eukaryotes including humans. The DNA methyltransferases (DNMTs) that establish and maintain 5mC, are conserved in many species across eukaryotes, including animals, plants, and fungi, mainly in a CpG context. Interestingly, 5mC levels and distributions are quite variable across phylogenies with some species even appearing to have no such DNA methylation.

This interesting and well-written paper discusses continuation of some of the authors' work published several years ago. In that previous paper, the laboratory demonstrated that DNA methylation pathways coevolved with DNA repair mechanisms, specifically with the alkylation repair system. Specifically, they discovered that DNMTs can introduce alkylation damage into DNA, specifically in the form of 3-methylcytosine (3mC). (This appears to be an error in the DNMT enzymatic mechanism where the generation 3mC as opposed to its preferred product 5-methylcytosine (5mC), is caused by the flipped target cytosine binding to the active site pocket of the DNMT in an inverted orientation.) The presence of 3mC is potentially toxic and can cause replication stress, which this paper suggests may explain the loss of DNA methylation in different species. They further showed that the ALKB2 enzyme plays a crucial role in repairing this alkylation damage, further emphasizing the link between DNA methylation and DNA repair.

The co-evolution of DNMTs with DNA repair mechanisms suggest there can be distinct advantages and disadvantages of DNA methylation to different species which might depend on their environmental niche. In environments that expose species to high levels of DNA damage, high levels of 5mC in their genome may be disadvantageous. This present paper sets out to examine the sensitivity of an organism to genotoxic stresses such as alkylation and oxidation agents as the consequence of DNMT activity. Since such a study in eukaryotes would be complicated by DNA methylation controlling gene regulation, these authors cleverly utilize Escherichia coli (E.coli) and incorporate into it the DNMTs from other bacteria that methylate the cytosines of DNA in a CpG context like that observed in eukaryotes; the active sites of these enzymes are very similar to eukaryotic DNMTs and basically utilize the same catalytic mechanism (also this strain of E.coli does not specifically degrade this methylated DNA) .

The experiments in this paper more than adequately show that E. coli expression of these DNMTs (comparing to the same strain without the DNMTS) do indeed show increased sensitivity to alkylating agents and this sensitivity was even greater than expected when a DNA repair mechanism was inactivated. Moreover, they show that this E. coli expressing this DNMT is more sensitive to oxidizing agents such as H2O2 and has exacerbated sensitivity when a DNA repair glycosylase is inactivated. Both propensities suggest that DNMT activity itself may generate additional genotoxic stress. Intrigued that DNMT expression itself might induce sensitivity to oxidative stress, the experimenters used a fluorescent sensor to show that H2O2 induced reactive oxygen species (ROS) are markedly enhanced with DNMT expression. Importantly, they show that DNMT expression alone gave rise to increased ROS amounts and both H2O2 addition and DNMT expression has greater effect that the linear combination of the two separately. They also carefully checked that the increased sensitivity to H2O2 was not potentially caused by some effect on gene expression of detoxification genes by DNMT expression and activity. Finally, by using mass spectroscopy, they show that DNMT expression led to production of the 5mC oxidation derivatives 5-hydroxymethylcytosine (5hmC) and 5-formylcytosine (5fC) in DNA. 5fC is a substrate for base excision repair while 5hmC is not; more 5fC was observed. Introduction of non-bacterial enzymes that produce 5hmC and 5fC into the DNMT expressing bacteria again showed a greater sensitivity than expected. Remarkedly, in their assay with addition of H2O2, bacteria showed no growth with this dual expression of DNMT and these enzymes.

Overall, the authors conduct well thought-out and simple experiments to show that a disadvantageous consequence of DNMT expression leading to 5mC in DNA is increased sensitivity to oxidative stress as well as alkylating agents.

Again, the paper is well-written and organized. The hypotheses are well-examined by simple experiments. The results are interesting and can impact many scientific areas such as our understanding of evolutionary pressures on an organism by environment to impacting our understanding about how environment of a malignant cell in the human body may lead to cancer.

In a new revised version of the paper, the authors have adequately addressed issues put forth by other reviewers.

Reviewer #3 (Public review):

Summary:

Krwawicz et al., present evidence that expression of DNMTs in E. coli results in (1) introduction of alkylation damage that is repaired by AlkB; (2) confers hypersensitivity to alkylating agents such as MMS (and exacerbated by loss of AlkB); (3) confers hypersensitivity to oxidative stress (H2O2 exposure); (4) results in a modest increase in ROS in the absence of exogenous H2O2 exposure; and (5) results in the production of oxidation products of 5mC, namely 5hmC and 5fC, leading to cellular toxicity. The findings reported here have interesting implications for the concept that such genotoxic and potentially mutagenic consequences of DNMT expression (resulting in 5mC) could be selectively disadvantageous for certain organisms. The other aspect of this work which is important for understanding the biological endpoints of genotoxic stress is the notion that DNA damage per se somehow induces elevated levels of ROS.

Strengths:

The manuscript is well-written, and the experiments have been carefully executed providing data that support the authors' proposed model presented in Fig. 7 (Discussion, sources of DNA damage due to DNMT expression).

Weaknesses:

(1) The authors have established an informative system relying on expression of DNMTs to gauge the effects of such expression and subsequent induction of 3mC and 5mC on cell survival and sensitivity to an alkylating agent (MMS) and exogenous oxidative stress (H2O2 exposure). The authors state (p4) that Fig. 2 shows that "Cells expressing either M.SssI or M.MpeI showed increased sensitivity to MMS treatment compared to WT C2523, supporting the conclusion that the expression of DNMTs increased the levels of alkylation damage." This is a confusing statement and requires revision as Fig. 2 does ALL cells shown in Fig. 2 are expressing DNMTs and have been treated with MMS. It is the absence of AlkB and the expression of DNMTs that that causes the MMS sensitivity.

(2) It would be important to know whether the increased sensitivity (toxicity) to DNMT expression and MMS is also accompanied by substantial increases in mutagenicity. The authors should explain in the text why mutation frequencies were not also measured in these experiments.

(3) Materials and Methods. ROS production monitoring. The "Total Reactive Oxygen Species (ROS) Assay Kit" has not been adequately described. Who is the Vendor? What is the nature of the ROS probes employed in this assay? Which specific ROS correspond to "total ROS"?

(4) The demonstration (Fig. 4) that DNMT expression results in elevated ROS and its further synergistic increase when cells are also exposed to H2O2 is the basis for the authors' discussion of DNA damage-induced increases in cellular ROS. S. cerevisiae does not possess DNMTs/5mC, yet exposure to MMS also results in substantial increases in intracellular ROS (Rowe et al, (2008) Free Rad. Biol. Med. 45:1167-1177. PMC2643028). The authors should be aware of previous studies that have linked DNA damage to intracellular increases in ROS in other organisms and should comment on this in the text.

Comments for the revised manuscript:

In this revised manuscript, the authors have satisfactorily addressed the issues raised in the review of the original submission and have significantly improved these studies.

Author response:

The following is the authors’ response to the previous reviews

Reviewer #1 (Public review):

I am not currently convinced by the principal interpretations and think that other explanations based on known phenomena could account for key results. Specifically the authors have not resolved whether oxidative modification to 5mC and 3mC, or chemical attack to ssDNA that is transiently exposed in the repair processing of 5mC and 3mC is the principal source of the observed genotoxicity.

(1) Original query which still stands: As noted in the manuscript, AlkB repairs alkylation damage by direct reversal (DNA strands are not cut). In the absence of AlkB, repair of alklylation damage/modification is likely through BER or other processes involving strand excision and resulting in single stranded DNA. It has previously been shown that 3mC modification from MMS exposure is highly specific to single stranded DNA (PMID:20663718) occurring at ~20,000 times the rate as double stranded DNA. Consequently the introduction of DNMTs is expected to introduce many methylation adducts genome-wide that will generate single stranded DNA tracts when repaired in an AlkB deficient background (but not in an AlkB WT background), which are then hyper-susceptible to attack by MMS. Such ssDNA tracts are also vulnerable to generating double strand breaks, especially when they contain DNA polymerase stalling adducts such as 3mC. The generation of ssDNA during repair is similarly expected follow the H2O2 or TET based conversion of 5mC to 5hmC or 5fC neither of which can be directly repaired and depend on single strand excision for their removal. The potential importance of ssDNA generation in the experiments has not been [adequately] considered.

We thank the reviewer for expanding on their previous comment.  We completely agree with the possibility that they raise and have added an extra paragraph in the discussion to expand on our consideration of the role of ssDNA in DNMT-induced DNA damage, which we reproduce here:

"The observation that TET overexpression sensitizes cells expressing DNMTs to oxidative stress strongly suggests that the site of DNA damage is the modified cytosine itself.  However, we do not currently have definitive evidence supporting this.  As mentioned in the results section, the presence of unrepaired 3mC may lead to increased levels of ssDNA; it is also possible that 5mC itself may increase ssDNA levels.  Loss of alkB would be expected to increase the amount of ssDNA.  Thus DNA damage surrounding modification sites, but not specifically localised to it, might be the cause of the increased sensitivity.  These two different models make different predictions.  If modified cytosines are the source of the damage, mutations arising would be predominantly located at CG dinucleotides.  Alternatively, ssDNA exposure would result in distributed mutations that would not necessarily be located at CG sites.  The highly biased spectrum of mutations that can be screened through the Rif resistance assay does not allow us to address this currently.  However, future experiments to create mutation accumulation lines could allow us to address the question systematically on a genome-wide level. "

Okay so I guess the right question to consider is: now I only give myself 2 weeks for this task, then how to make the most out of it. This is the right problem frame! Now we proceed to tackling it. I believe the Parento 20/80 law is one of the correct answer. Any others?

this is super important

What is this? It's might be a very important idea to overcome the problem of evolving tree structure!

A lot of vague words are used. Unlike the other speech, Trump's speech doesn't use history for context or draw from any specifics. It's a lot of vague words and promises.

strawman

Begging the question.

Really adds to how much hatred and violence there is.

This shows how in all of history black resident are placed in urban areas while white people controlled bigger cities and live in suburbs. A continuing cycle.

Enfranchisement=the right to vote

Greenbacks=paper issued by the government (our current dollar bills!)

Sharecropping= white owners keeping freed slaves under contract to basically keep them as slaves again

Really hones in on the idea of how after the war the United States was still struggling.

This shows me how the Civil war never really stopped after the surrender because the issues were still there. Only the war stopped but the problems were never fixed, not until much later.

Really talking about how the Civil war barely started the conversation about race and equality. Black Americans everywhere but especially in the south had to fight for their freedom for many centuries after the Civil war.

This comment is mostly about enslaved people and their future. Especially their rights moving on and people reevaluating what the constitution really means.

I think this is how most of the United States felt after the civil war. There was much uncertainty about the future of enslaved people and the future of the states. This is the reconstruction era.

This is were reconstruction actually begins to show its positive progress of integrating African Americans as free citizens.

So the two parties and their beliefs were reversed in the past (?)

identification term

Because of newfound "freedom", states found loopholes to keep African Americans under any and all restrictions to still hold power over them.

identification term

They were the exception to the rule basically

Establishing new founded equality and citizenship for those endured slavery would take a long time until they began to be more and more accepted and supported.

Would they return to the way things were? Or would former slaves be given stepping stones toward equality and recognition as an American Citizen.

should be annotate-able

experiment.with-short.names

// moving forward by using public annotations on the margin on open web accessible working document

on a personal note I've been scattering ideas in annotations for years, assuming that one they all that will be hyper mapped

now with this arrangement detailed bellow I can reflect annotate and provide links to My own MEMEplEXes and working docments as I bootstrap the IndyWeb

since hyperpost.peergos.me is now accessible by using via.hypothes.is we an do that even before the via.indy.web.annotation capability is ready to launched

不一致就不能平行运算吗？

I personally believe that access to information is key to a more well informed and a smarter society. Without open information then people can be pushed into narrow views and cannot form their own opinions.

The idea of the Open Pedagogy is an interesting one in the fact that rather than just showing you the definition and having you absorb the information and nothing else, it encourages the student to engage with it and write things in their own creative and interesting ways rather than the way a lot of institutions intend for it to be read. This gets me excited to put this into action with future assignments.

One question I have with the Open Pedagogy is how exactly would we be able to properly document all of these new ideas students present with this method? cause I would definitely like to see the different ideas and outlooks documented in one place

I though the idea that "Open Pedagogy" has where students are the creators of knowledge rather than passive consumers of it is amazing because this is the first time I have ever heard of this happening with my time in school and I'm excited for the book. I'm also looking forward to comparing how the students wrote their book compared to a traditional way. In the terms of "Open Pedagogy", I wonder how this idea would work on other subject books?

The “family room” was made to bring everyone together. It is kind of like how we hang out in the living room today to relax or talk.

That is a huge number and shows how badly support was needed after the war.

These movies showed men who were struggling emotionally or with addiction. It’s interesting that TV and film in the 1950s didn’t just show perfect families, they also showed serious problems inside the home.

After the war, people were confused about what men and women were supposed to do at home. Families were expected to look perfect, but real life didn’t always match that.

Field Access: Structs are NOT the same as Maps

The comment "(same as maps)" is misleading.

• Maps: Support both map.key and map[:key] access
• Structs: Only support struct.key dot notation
• Bracket notation struct[:key] is not work with structs

Proof

It says %User{name: "Alice"} fails at compile time with a missing required field error, but when I tested it, it works fine and just sets the other fields to nil. There doesn't seem to be any required field validation happening here.

This is pointing to the required fields instead of the fields with defaults

I had no Idea you could plagarize yourself. Each assignment you complete needs to be completly unique for the class.

In this passage, Gilgamesh dreams of a figure who will be his equal, and interpreters tell him that “heroes kiss his feet” and that he will lead Gilgamesh. This prophetic dream frames Enkidu as Gilgamesh’s destined counterpart: not only a rival but also a partner who will shape his heroic identity. The imagery of “an axe in the street” and “heroes kiss his feet” reflects how masculinity is tied to symbols of power and violence, yet also reverence. The text suggests that Gilgamesh’s greatness requires balance. Gender politics emerge through the absence of women in this dream: the hero’s destiny is mediated entirely through male bonds. The translation describe Enkidu as Gilgamesh’s “likeness,” collapsing rivalry into mirror-image intimacy. Gilgamesh’s heroic identity is forged in masculine struggle and mutual recognition.

Camron Newcomb

CC BY-NC-SA 4.0

This moment in the Old Babylonian version underscores how gendered power is central to the hero making process in early Mesopotamian culture. Shamhat, the unnamed "harlot," initiates Enkidu's transformation from wild beast to man, and then from man to hero, not through brute force, but by teaching him to conform to gendered norms of civilization.

Importantly, civilization here is gendered male: Enkidu must learn to eat bread, drink milk, wear clothes, and accept hierarchy, including the authority of the male king, Gilgamesh. This socialization is mediated by a woman, but it ultimately renders women peripheral once male heroism is established. Even when Enkidu and Gilgamesh bond, it is through violent competition and mutual respect, culminating in a moment where Enkidu prevents Gilgamesh from pursuing the goddess Išhara, framing love or femininity as a threat to masculine heroic purpose.

Clay and Jastrow’s 1920 translation reflects early 20th century ideas about gender and morality. Their diction treats the "harlot" with subtle moral judgment, while placing more noble framing around the “mighty hunter” Enkidu. The translation also shows a preference for structured, formalized syntax, which reinforces the patriarchal lens through which the epic was interpreted at the time.

If you continue to practice/review material over an extended period of time even after a test, you will retain more information as well as relearn the material quicker after some time. (More time spent on subject = more mastery)

This is quite similar to the concept of schemas.

This might be off topic, but this reminds me of how studies have shown that students who chew gum or listen to a specific playlist while studying tend to perform better on exams when they chew the same gum or listen to the same music during the test as well. Just an interesting thought that might correlate to the study.

We should try to obtain observable data and methods that will allow us to specifically measure the process of learning and its effectiveness, instead of just looking at the results from tests and exams.

Countering the common rebuttal of erasing history to introduce the idea of a better future. He's reframing the debate. Rather than tearing down the statues to erase history, instead it becomes tearing them down out of an understanding of history.

The image of a black child passing by a Civil War monument on his way to school creates pathos in Landrieu's argument. Children's innocence is a stark opposition to the violent history represented by the monument.

An example of how diversity not only makes things better, but is the reality of America. In mentioning "everything in the pot" he is gesturing at the "melting pot" story of America, which brings in some Patriotism. He continues to comment on America's greatness perhaps to balance his argument and protect it from those would suggest he is simply being unpatriotic.

"they called you everything but a child of God" is a variation of a southern saying. "They called me everything but a child of God" or, the more rocky "They called me everything but a white man." Landrieu is building his ethos. He is from here. I see a lot of this within his speech. It is an effective way to argue an unpopular and polarizing topic.

Landrieu's speech is coming from a person from this community. The ease of fighting online has created a highly polarized us vs. them mentality. I view this line as a way Landrieu to continue this idea that he is from here. He is not from away. There are many different religions in New Orleans and there are many people "from away" that practice Christianity, but there is a significant amount of Christian ethos in the south.

This image of Phaedra and Hippolytus reflects the central conflict of Euripides’ tragedy: Phaedra’s desire and Hippolytus’ resistance. Phaedra embodies passion, shame, and transgression. Hippolytus, in contrast, who represents purity, self-control, and loyalty especially to Artemis. Phaedra’s speech is described in terms of “madness,” “disease,” or “frenzy,” while Hippolytus’ refusal is couched in terms of “virtue” and “nobility.” The politics of language preserve a worldview where male strength lies in resisting women, casting the hero as morally elevated only through female exclusion.

© 2025 Melinessa Louis Douze. Licensed under Creative Commons Attribution 4.0 International (CC BY 4.0)

In Phaedra, gender roles are central to the tragedy, especially through the contrast between Phaedra's forbidden desire and Hippolytus's proud purity. Phaedra is portrayed as emotionally and sexually unstable, her desire treated as both dangerous and shameful. Her downfall reinforces patriarchal views where female sexuality must be hidden or punished. Meanwhile, Hippolytus's heroism lies in his control and rejection of passion, fitting the Greek ideal of masculine virtue. rational, proud, and emotionally restrained. Compared to Sita Sings the Blues, Phaedra is a woman destroyed by her feelings, while Sita is a woman silenced by social duty but both are trapped in male-dominated systems that define a hero through emotional suppression or moral superiority. Sita, especially in Paley's version, is allowed to speak back, while Phaedra's voice leads to her ruin.

© 2025 Melinessa Louis Douze. Licensed under Creative Commons Attribution 4.0 International (CC BY 4.0).

this is surprising to me because i feel like the number would be lower than 60% i feel like technology had spread all over the world.

I am actually nervous that AI intelligent robots may take over most of the jobs like uber, grocery stores. We already have the Self driving car called the tesla.

I wonder how these people are impacted throughout there daily lives without the internet how do they learn?

By crafting Enkidu to match Gilgamesh’s “stormy heart,” the gods frame male power as something wild, aggressive, and potentially dangerous unless checked by another man of equal force. The word “stormy” conveys emotional turbulence, suggesting that admired manhood in Mesopotamian culture was intense, unpredictable. Peace in Uruk is imagined not as communal cooperation but as the result of two men clashing until balance is achieved. This emphasis on physical struggle reflects a patriarchal worldview where masculinity is proven by combat and domination.

© 2025 Melinessa Louis Douze. Licensed under Creative Commons Attribution 4.0 International (CC BY 4.0)

The harlot’s invitation is an important turning point because it shows a woman actively guiding the male hero’s path rather than existing only as a passive figure. Shamhat uses sexuality as a form of persuasion, but the translation’s choice of the word “harlot” colors her power with moral suspicion, echoing patriarchal anxieties about female influence. Instead of being merely an object of desire, she functions as a bridge between wilderness and civilization, embodying beauty, culture, and religious order. This suggests that female sexuality is not only potent but also necessary for shaping male strength into socialized heroism. Enkidu’s willingness to follow her into Uruk and challenge Gilgamesh shows that the epic constructs heroism as relational male power defined in response to both female influence and urban culture. Unlike the Ramayana, where Sita embodies loyalty and sacrifice, Shamhat’s role is active and influential, showing how Mesopotamian traditions allowed women to act as agents of transformation, even if through sexuality framed as “dangerous.” The language of translation here is crucial: by choosing “harlot,” the text imposes judgment on Shamhat, reinforcing a patriarchal reading that might not fully capture her cultural role as a temple courtesan.

© 2025 Melinessa Louis Douze. Licensed under Creative Commons Attribution 4.0 International (CC BY 4.0).

This information shows that specialization and trade rely mainly in banks and other financial helpers. It's kind of like when people borrow money and don't know how to take care of it. When people trust each other, it is so much easier to do their jobs. However, my question is, why is trust in banks more important than keeping a specialized economy running? I feel like in our economy today, those ideas on specialization are really important.

I think that this section shows the sign of "specialization" because Smith shows the readers that when workers focused on specific tasks that were given, every single day they became more productive than when they tried to do everything themselves. This examples shows the importance of how even task increases efficiency and output.

Adam Smith says, “the number of people of whose industry a part has been employed in procuring him this accommodation, exceeds all computation.” This means even a simple worker depends on many other people. The first society, families do everything themselves grow food, make clothes, and build shelter. In a modern economy, no one can do everything alone. Most of what you use is made by many other people.

How does depending on many people change how we see the value of different jobs compared to a self-sufficient society?

In this paragraph the author says, “economists have lost the art of critical thinking” and “The critical thinker always asks, ‘How do you know that?’” He is criticizing how economists think and the assumptions they make. In the Stanford passage, the author says, “An economy should produce enough goods and services to support its citizens” and “Economic progress requires us to think continuously about how to make our work more effective and useful.” He focuses on practical results and improving everyday life.

Comparison: Both authors agree economics should connect to real life. The Kling author looks at how economists think, while the Stanford author looks at what the economy does for people.

This is our first look into the main theme of the text which is racism can lead to diminishing of ones identity

In this part of the text, it is very evident where and how philosophy plays a role in physics and how it came to be. We can see that this " human endeavor" is due to our curiosity in the way the world/ universe works. It is based off our observations and how we question. For example, questions like what is time? What is space? or do particles exist independently of observation? Are as much philosophical as they are physical.

why ?

The Mercantilists in 1600s Britain believed countries grow richer by exporting more than they import. They supported strong governments to control trade and colonies. Their ideas still influence modern “export-led growth” in countries like Germany, Korea, and China.

In the late 1970s, the “Golden Age” ended, and capitalism became more aggressive. Businesses and investors opposed taxes, social programs, unions, and regulations. Falling profits, high inflation, and global instability, including left-wing revolutions, led them to push for a stricter approach and they got it.

After World War II, capitalism entered a very successful period called the “Golden Age.” For about thirty years, wages and living standards doubled in developed countries. Businesses invested a lot, governments spent more, unemployment was low, and productivity grew, while profits stayed strong. During the Cold War, governments and business leaders also supported more equality and security to show that capitalism worked better than communism.

Not all capitalist economies are the same. The “Anglo-Saxon” model has high inequality, a small government, and a big financial sector. Other countries, like Nordic or Asian ones, give better outcomes for workers. This shows countries can make their economies fairer, even in capitalism, though global pressures still matter.

Economic progress requires us to keep finding better ways to do our work. This is called “innovation” and includes creating new products and improving how we make them. An economy should support innovation, or it will run out of new ideas. The author shows that innovation is important, which makes you curious about how different countries or systems help people be creative and what policies or practices encourage new ideas.

The economy is simultaneously mystifying and straightforward. The author suggests it can be complex but also understandable in everyday life. Everyone has experience with the economy and participates in it, which assumes personal experience gives a valid understanding. The forces and relationships you notice locally matter more than stock market changes, though bigger economic developments still affect daily life.

If personal experience shapes our understanding of the economy, what might we miss about how national or global policies affect our neighborhoods?

Most economists follow one type of economics called neoclassical economics. This kind of economics mixes beliefs with science. It was created in the late 1800s to support capitalism, not just explain it. (It favors the rich and powerful.) It still tries to “prove” ideas that aren’t true, like saying owning money is productive, everyone gets paid fairly, or unemployment isn’t real.

this shows how it makes economics sound like it can be explained just by looking at everyday things like work, shopping, tools, money, and the environment. While these are important parts of the economy, reducing economics to only what we see around us leaves out .

I fell like it is good to encourage people to learn about the economy, this view risks making economics seem too simple. Every single day people may have an interest in the economy, but that does not mean they can easily understand its complex systems. But you should always have the right knowledge because sometimes people will get the wrong idea.

I feel like everyone is within the economy but not everyone understands how things work. but the risk can be implied of the discipline of the economy. personal experiences may give people opinions about the economy, but that doesn’t mean they understand complex issues like government spending, trade, or banking. saying “everyone contributes” also ignores the fact that some people and groups have much more power to shape the economy than others

The Andromeda Galaxy seems so out of reach and disconnected, but it is interesting to see how physics is the one thing that connects us all as a part of the universe and nature. We think that everything is so out of reach, and we sometimes get so caught up in our own world, but our world is connected to others in ways we can't even see. However, we know that through these forces it is all related. We are all connected.

Weirdly enough, I always associate physics with much more physical (badum tss) reactions. In all reality, I think it stems solely from a vine that was popular before the death of the app that featured a high schooler showing her physics teacher shooting off what looked like a rocket in the classroom. That, followed by endless online compilations of teachers and professors showing a plethora of crazy things tends to give me a more exciting view of physics than what it actually is. There is a severe lack of Newton's cradles and electrical reactions, and instead, a lot more math and graphing. I understand that pulling the most interesting bits of a topic is what garners interest, but in alll my experience with physics, it more so feels like false advertising: like being told you're going to a waterpark and arriving to see a pitiful slip n' slide. I can't entirely say I'm sure what classifies as physics, as like I said before, my exposure has been so incredibly broad that it feels impossible to pinpoint the boxes that need to be checked for something to fall into that category. To make an overly long blog-post short, 'physics' makes me think of fun science-y stuff, but also leaves me both confused and disappointed by the reality of it.

I'm pretty sure Plato might say that physics gives us a glimpse of the shadows on the cave wall, the measurable patterns of motion, matter, and energy, but it also points us toward something deeper: the unchanging principles, or ‘Forms,’ that underlie the shifting physical world. In that sense, studying physics isn’t just about particles and forces; it’s about recognizing that everything is connected by a hidden order, a kind of mathematical harmony that reflects a greater reality beyond what we see. At least When I read this statement I felt as if the book hinted at something philosophical. but perhaps the book simply means "everything we can see and know" and not in a metaphysical or spiritual emotional "see and know."

I never knew that the ancient Greeks were among the first to try to understand nature. It was also very interesting to find out that physics stemmed from the prior knowledge of astronomy and mechanics.

throughout my college experience i have always had access to online textbook and the public library.

What are some extra steps in our list we can take for responding to off task behavior to avoid calling home and referral to the office? To give sympathy for a student struggling to behave.

I like the idea of involving students in setting classroom rules because it holds them accountable and gives them a sense of responsibility to uphold what they suggest.

annotations help me remember what i just read and how my brain reacts to the information.

I love canvas i use this for every project that i have it helps me organize my information.

I use my email everyday for my school work I think that everyone uses there email at least once in there lives.

I feel like social media has spread lots of fake news and so had the media in general.

I feel like this goes along with stereotypes about people that we tend to make up. We can even talk ourselves out of good ideas or plans because of conformational bias.

I always valued our technology that we were able to use throughout the semester. I could never afford a computer on my own and highschool always provided me with useful technology.

This is very true because many countries to not have access to the internet or technology.

Rehabilitation: The set of interventions undertaken to help a criminal with their criminal behavioral patterns.

Recidivism: The recurrent crimes committed by the same individuals after release from prison.

eLife Assessment

This study presents DeepTX, a valuable methodological tool that integrates mechanistic stochastic models with single-cell RNA sequencing data to infer transcriptional burst kinetics at genome scale. The approach is broadly applicable and of interest to subfields such as systems biology, bioinformatics, and gene regulation. The evidence supporting the findings is solid, with appropriate validation on synthetic data and thoughtful discussion of limitations related to identifiability and model assumptions.

Joint Public Review:

In this work, the authors present DeepTX, a computational tool for studying transcriptional bursting using single-cell RNA sequencing (scRNA-seq) data and deep learning. The method aims to infer transcriptional burst dynamics-including key model parameters and the associated steady-state distributions-directly from noisy single-cell data. The authors apply DeepTX to datasets from DNA damage experiments, revealing distinct regulatory patterns: IdU treatment in mouse stem cells increases burst size, promoting differentiation, while 5FU alters burst frequency in human cancer cells, driving apoptosis or survival depending on dose. These findings underscore the role of burst regulation in mediating cell fate responses to DNA damage.

The main strength of this study lies in its methodological contribution. DeepTX integrates a non-Markovian mechanistic model with deep learning to approximate steady-state mRNA distributions as mixtures of negative binomial distributions, enabling genome-scale parameter inference with reduced computational cost. The authors provide a clear discussion of the framework's assumptions, including reliance on steady-state data and the inherent unidentifiability of parameter sets, and they outline how the model could be extended to other regulatory processes.

The revised manuscript addresses many of the original concerns, particularly regarding sample size requirements, distributional assumptions, and the biological interpretation of inferred parameters. However, the framework remains limited by the constraints of snapshot data and cannot yet resolve dynamic heterogeneity or causality. The manuscript would also benefit from a broader contextualisation of DeepTX within the landscape of existing tools linking mechanistic modelling and single-cell transcriptomics. Finally, the interpretation of pathway enrichment analyses still warrants clarification.

Overall, this work represents a valuable contribution to the integration of mechanistic models with high-dimensional single-cell data. It will be of interest to researchers in systems biology, bioinformatics, and computational modelling.

Author response:

The following is the authors’ response to the original reviews

Joint Public Review:

In this work, the authors develop a new computational tool, DeepTX, for studying transcriptional bursting through the analysis of single-cell RNA sequencing (scRNA-seq) data using deep learning techniques. This tool aims to describe and predict the transcriptional bursting mechanism, including key model parameters and the steady-state distribution associated with the predicted parameters. By leveraging scRNA-seq data, DeepTX provides high-resolution transcriptional information at the single-cell level, despite the presence of noise that can cause gene expression variation. The authors apply DeepTX to DNA damage experiments, revealing distinct cellular responses based on transcriptional burst kinetics. Specifically, IdU treatment in mouse stem cells increases burst size, promoting differentiation, while 5FU affects burst frequency in human cancer cells, leading to apoptosis or, depending on the dose, to survival and potential drug resistance. These findings underscore the fundamental role of transcriptional burst regulation in cellular responses to DNA damage, including cell differentiation, apoptosis, and survival. Although the insights provided by this tool are mostly well supported by the authors' methods, certain aspects would benefit from further clarification.

The strengths of this paper lie in its methodological advancements and potential broad applicability. By employing the DeepTXSolver neural network, the authors efficiently approximate stationary distributions of mRNA counts through a mixture of negative binomial distributions, establishing a simple yet accurate mapping between the kinetic parameters of the mechanistic model and the resulting steady-state distributions. This innovative use of neural networks allows for efficient inference of kinetic parameters with DeepTXInferrer, reducing computational costs significantly for complex, multi-gene models. The approach advances parameter estimation for high-dimensional datasets, leveraging the power of deep learning to overcome the computational expense typically associated with stochastic mechanistic models. Beyond its current application to DNA damage responses, the tool can be adapted to explore transcriptional changes due to various biological factors, making it valuable to the systems biology, bioinformatics, and mechanistic modelling communities. Additionally, this work contributes to the integration of mechanistic modelling and -omics data, a vital area in achieving deeper insights into biological systems at the cellular and molecular levels.  

We thank the reviewers for their positive opinion on our manuscript. As reflected in our detailed responses to the reviewers’ comments, we will make significant changes to address their concerns comprehensively.

This work also presents some weaknesses, particularly concerning specific technical aspects. The tool was validated using synthetic data, and while it can predict parameters and steady-state distributions that explain gene expression behaviour across many genes, it requires substantial data for training. The authors account for measurement noise in the parameter inference process, which is commendable, yet they do not specify the exact number of samples required to achieve reliable predictions. Moreover, the tool has limitations arising from assumptions made in its design, such as assuming that gene expression counts for the same cell type follow a consistent distribution. This assumption may not hold in cases where RNA measurement timing introduces variability in expression profiles.

Thank reviewers for detailed and constructive feedback on our work. We will address the key concerns raised from the following points:

(1) Clarification on the required sample size: We tested the robustness of our inference method on simulated datasets by varying the number of single-cell samples. Our results indicated that the predictions of burst kinetics parameters become accurate when the number of cells reaches 500 (Supplementary Figure S3d, e). This sample size is smaller than the data typically obtained with current single-cell RNA sequencing (scRNA-seq) technologies, such as 10x Genomics and Smart-seq3 (Zheng GX et al., 2017; Hagemann-Jensen M et al., 2020). Therefore, we believed that our algorithm is well-suited for inferring burst kinetics from existing scRNA-seq datasets, where the sample size is sufficient for reliable predictions. We will clarify this point in the main text to make it easier for readers to use the tool.

(2) Assumption-related limitations: One of the fundamental assumptions in our study is that the expression counts of each gene are independently and identically distributed (i.i.d.) among cells, which is a commonly adopted assumption in many related works (Larsson AJM et al., 2019; Ochiai H et al., 2020; Luo S et al., 2023). However, we acknowledged the limitations of this assumption. The expression counts of the same gene in each cell may follow distinct distributions even from the same cell type, and dependencies between genes could exist in realistic biological processes. We recognized this and will deeply discuss these limitations from assumptions and prospect as an important direction for future research.  

The authors present a deep learning pipeline to predict the steady-state distribution, model parameters, and statistical measures solely from scRNA-seq data. Results across three datasets appear robust, indicating that the tool successfully identifies genes associated with expression variability and generates consistent distributions based on its parameters. However, it remains unclear whether these results are sufficient to fully characterise the transcriptional bursting parameter space. The parameters identified by the tool pertain only to the steady-state distribution of the observed data, without ensuring that this distribution specifically originates from transcriptional bursting dynamics.

We appreciate reviewers’ comments and the opportunity to clarify our study’s contributions and limitations. Although we agree that assessing whether the results from these three realistic datasets can represent the characterize transcriptional burst parameter space is challenging, as it depends on data property and conditions in biology, we firmly believe that DeepTX has the capacity to characterize the full parameter space. This believes stems from the extensive parameters and samples we input during model training and inference across a sufficiently large parameter range (Method 1.3). Furthermore, the training of the model is both flexible and scalable, allowing for the expansion of the transcriptional burst parameter space as needed. We will clarify this in the text to enable readers to use DeepTX more flexibly.

On the other hand, we agree that parameter identification is based on the steady-state distribution of the observed data (static data), which loses information about the fine dynamic process of the burst kinetics. In principle, tracking the gene expression of living cells can provide the most complete information about real-time transcriptional dynamics across various timescales (Rodriguez J et al., 2019).

However, it is typically limited to only a small number of genes and cells, which could not investigate general principles of transcriptional burst kinetics on a genome-wide scale. Therefore, leveraging the both steady-state distribution of scRNA-seq data and mathematical dynamic modelling to infer genome-wide transcriptional bursting dynamics represents a critical and emerging frontier in this field. For example, the statistical inference framework based on the Markovian telegraph model, as demonstrated in (Larsson AJM et al., 2019), offers a valuable paradigm for understanding underlying transcriptional bursting mechanisms. Building on this, our study considered a more generalized non-Mordovian model that better captures transcriptional kinetics by employing deep learning method under conditions such as DNA damage. This provided a powerful framework for comparative analyses of how DNA damage induces alterations in transcriptional bursting kinetics across the genome. We will highlight the limitations of current inference using steady-state distributions in the text and look ahead to future research directions for inference using time series data across the genome.

A primary concern with the TXmodel is its reliance on four independent parameters to describe gene state-switching dynamics. Although this general model can capture specific cases, such as the refractory and telegraph models, accurately estimating the parameters of the refractory model using only steadystate distributions and typical cell counts proves challenging in the absence of time-dependent data.

We thank reviewers for highlighting this critical concern regarding the TXmodel's reliance on four independent parameters to describe gene state-switching dynamics. We acknowledge that estimating the parameters of the TXmodel using only steady-state distributions and typical single-cell RNA sequencing (scRNA-seq) data poses significant challenges, particularly in the absence of timeresolved measurements.

As described in the response of last point, while time-resolved data can provide richer information than static scRNA-seq data, it is currently limited to a small number of genes and cells, whereas static scRNA-seq data typically capture genome-wide expression. Our framework leverages deep learning methods to link mechanistic models with static scRNA-seq data, enabling the inference of genome-wide dynamic behaviors of genes. This provides a potential pathway for comparative analyses of transcriptional bursting kinetics across the entire genome.

Nonetheless, the refractory model and telegraphic model are important models for studying transcription bursts. We will discuss and compare them in terms of the accuracy of inferred parameters.

Certainly, we agree that inferring the molecular mechanisms underlying transcriptional burst kinetics using time-resolved data remains a critical future direction. We will include a brief discussion on the role and importance of time-resolved data in addressing these challenges in the discussion section of the revised manuscript.

The claim that the GO analysis pertains specifically to DNA damage response signal transduction and cell cycle G2/M phase transition is not fully accurate. In reality, the GO analysis yielded stronger p-values for pathways related to the mitotic cell cycle checkpoint signalling. As presented, the GO analysis serves more as a preliminary starting point for further bioinformatics investigation that could substantiate these conclusions. Additionally, while GSEA analysis was performed following the GO analysis, the involvement of the cardiac muscle cell differentiation pathway remains unclear, as it was not among the GO terms identified in the initial GO analysis.

We thank the reviewer for this valuable feedback and for pointing out the need for clarification regarding the GO and GSEA analyses. We agree that the connection between the cardiac muscle cell differentiation pathway identified in the GSEA analysis and the GO terms from the initial analysis requires further clarification. This discrepancy arises because GSEA examines broader sets of pathways and may capture biological processes not highlighted by GO analysis due to differences in the statistical methods and pathway definitions used. We will revise the manuscript to address this point, explicitly discussing the distinct yet complementary nature of GO and GSEA analyses and providing a clearer interpretation of the results.

As the advancement is primarily methodological, it lacks a comprehensive comparison with traditional methods that serve similar functions. Consequently, the overall evaluation of the method, including aspects such as inference accuracy, computational efficiency, and memory cost, remains unclear. The paper would benefit from being contextualised alongside other computational tools aimed at integrating mechanistic modelling with single-cell RNA sequencing data. Additional context regarding the advantages of deep learning methods, the challenges of analysing large, high-dimensional datasets, and the complexities of parameter estimation for intricate models would strengthen the work.

We greatly appreciate your insightful feedback, which highlights important considerations for evaluating and contextualizing our methodological advancements. Below, we emphasize our advantages from both the modeling perspective and the inference perspective compared with previous model. As our work is rooted in a model-based approach to describe the transcriptional bursting process underlying gene expression, the classic telegraph model (Markovian) and non-Markovian models which are commonly employed are suitable for this purpose:

Classic telegraph model: The classic telegraph model allows for the derivation of approximate analytical solutions through numerical integration, enabling efficient parameter point estimation via maximum likelihood methods, e.g., as explored in (Larsson AJM et al., 2019). Although exact analytical solutions for the telegraph model are not available, certain moments of its distribution can be explicitly derived. This allows for an alternative approach to parameter inference using moment-based estimation methods, e.g., as explored in (Ochiai H et al., 2020). However, it is important to note that higher-order sample moments can be unstable, potentially leading to significant estimation bias. 

Non-Markovian Models: For non-Markovian models, analytical or approximate analytical solutions remain elusive. Previous work has employed pseudo-likelihood approaches, leveraging statistical properties of the model’s solutions to estimate parameters ,e.g., as explored in (Luo S et al., 2023).

However, the method may suffer from low inference efficiency. 

In our current work, we leverage deep learning to estimate parameters of TXmodel, which is nonMarkovian model. First, we represent the model's solution as a mixture of negative binomial distributions, which is obtained by the deep learning method. Second, through integration with the deep learning architecture, the model parameters can be optimized using automatic differentiation, significantly improving inference efficiency. Furthermore, by employing a Bayesian framework, our method provides posterior distributions for the estimated dynamic parameters, offering a comprehensive characterization of uncertainty. Compared to traditional methods such as moment-based estimation or pseudo-likelihood approaches, we believe our approach not only achieves higher inference efficiency but also delivers posterior distributions for kinetics parameters, enhancing the interpretability and robustness of the results. We will present and emphasize the computational efficiency and memory cost of our methods the revised version.

Recommendations for the authors:

There are various noise sources in biological progress. How transcriptional bursting fits within those as well as the reasons to focus only on this source needs to be clearly discussed in the introduction of the manuscript. Related to this last point, transcriptional bursting might not be the only mechanism to take advantage of the stochastic nature of biomolecular processes to make decisions. Once again, what are the implications of assuming this as the underlying mechanism?

Thank the reviewer for this valuable comment. We fully agree that biological systems are subject to multiple stochastic sources, which arise from both intrinsic and extrinsic noise (Eling N et al., 2019). Intrinsic noise is primarily driven by the stochastic biochemical effects that directly influence mRNA and protein expression in a gene-specific manner, such as DNA, epigenetic, transcription, and translation levels. Extrinsic noise arises from fluctuations in cell-specific manners, such as changes in cell size, cell cycle, or cell signaling. Given that DNA damage most directly perturbs transcription and translation processes, focusing on intrinsic noise sources is appropriate for mechanistically modeling gene-specific expression variability, particularly since this variability can be captured at the genome-wide scale by scRNA-seq data.

Among various intrinsic noise sources, transcriptional bursting offers a mechanistically wellcharacterized and quantifiable representation of gene expression variability (Tunnacliffe E & Chubb JR, 2020). It reflects the dynamic switching between active and inactive gene states and has been observed consistently across prokaryotic and eukaryotic cells (Eling N et al., 2019). Moreover, transcriptional bursting kinetics, defined by burst size and frequency, can be inferred from scRNA-seq data at the singlegene level using steady-state assumptions, making it an analytically tractable and biologically meaningful feature for large-scale inference (Rodriguez J & Larson DR, 2020).

We acknowledge that transcriptional bursting is not the only mechanism through which cells can utilize stochasticity for fate decisions. Other processes, such as translational noise and chromatin accessibility, may also contribute. However, given the data modality (static scRNA-seq) and the established theoretical framework for bursting, we assume transcriptional bursting as a representative and interpretable proxy of stochastic regulation. This assumption enables us to extract meaningful insights while remaining open to future model extensions, incorporating additional regulatory layers as more data types become available.

In this version of the manuscript, we have revised the introduction section to better clarify the rationale of this assumption and to more explicitly emphasize the important role of transcriptional bursting within stochastic noise.

More careful discussion of how the proposed method differentiates from previous work that employs scRNA-seq to elucidate the diverse sources of noise (pp.3).

Thank the reviewer for this suggestion. Our proposed method differs significantly from previous work that utilizes scRNA-seq data to study diverse noise sources from several aspects (Ochiai H et al., 2020; Eling N et al., 2019; Morgan MD & Marioni JC, 2018). Specifically, DeepTX infers genomewide burst kinetics by directly matching the full steady-state distribution of a mechanistic stochastic model to the observed scRNA-seq data, rather than relying solely on low-order statistics such as mean and variance. Moreover, by adopting a non-Markovian process that allows multi-step promoter switching, DeepTX extends beyond the classic telegraph model to better capture the complex molecular events underlying transcriptional activation and repression. Crucially, we used a deep-learning–based solver to obtain these intractable steady-state distributions rapidly and accurately. This combination of richer data usage, more realistic mechanistic assumptions, and scalable neural-network–accelerated computation lays the groundwork for incorporating additional noise sources into a unified inference framework in future work. 

In this version of the manuscript, we have revised the discussion section to highlight the difference with previous works.

The paper could benefit from being contextualised alongside other computational tools that aim to integrate mechanistic modelling with single-cell RNA sequencing data. This is an active area of research, and works such as Sukys and Grima (bioRxiv, 2024), Garrido-Rodriguez et al. (PLOS Computational Biology, 2021), Maizels (2024), and others could provide valuable context.

Thank the reviewer for suggesting these relevant works. Garrido-Rodriguez et al. (PLOS Comput. Biol., 2021) integrated single-cell and bulk transcriptomic data into mechanistic pathway models to infer signaling dynamics, an approach complementary to our mapping of burst kinetic parameters onto pathway enrichment for linking transcriptional bursting to functional outcomes. Sukys and Grima et al. (bioRxiv, 2024; Now in Nucleic Acids Res., 2025) demonstrated that cell-cycle stage and cellular age significantly modulate burst frequency and size, highlighting the potential to enhance DeepTX by incorporating cell-cycle–dependent variability into genome-wide burst inference. Maizels et al. (Philos. Trans. R. Soc. Lond. B. Biol. Sci., 2024) reviewed methods for capturing single-cell temporal dynamics across multi-omic modalities, underscoring how higher time-resolved data could refine and validate steady-state burst inference frameworks to better resolve causal gene-expression mechanisms.

We have cited these studies on the contextual relevance to DeepTX in the discussion sections.

As the advancement is primarily methodological, it lacks a comprehensive comparison with traditional methods that serve similar functions. Consequently, the overall evaluation of the method, including aspects such as inference accuracy, computational efficiency, and memory cost, remains unclear. We suggest incorporating these experiments to provide readers with a more complete understanding of the proposed method's performance.

Thank the reviewer for constructive suggestion regarding a comprehensive comparison with other previous methods. To address this problem, in this version, we compared DeepTX with our previous work, txABC, that utilized approximate Bayesian computation to infer parameters from the generalized telegraph model (Luo S et al., 2023). As a result, DeepTX achieved improvements in inference accuracy and computational efficiency (Supplementary Figure S4.). For memory cost during single-gene inference, DeepTX requires an average memory usage of approximately 70 MB, whose memory consumption accounts for only a small fraction of the total available memory on standard computing devices (typically exceeding 10 GB), while exhibiting superior inference efficiency compared to txABC. We have mentioned in the third result section.

Discuss the validity of the assumption of the static snapshot provided by the scRNA-seq data as in steadystate (i.e., stationary distribution), and the implications of this assumption being untrue (for the proposed method).

We thank the reviewer for the comment regarding the stationary assumption. We assume that each scRNA-seq snapshot approximates the steady-state (stationary) distribution of transcript counts because (i) typical single-cell experiments sample large, asynchronously dividing populations that collectively traverse many transcriptional burst cycles, and (ii) in the absence of a synchronized perturbation, mRNA production and degradation reach a dynamic balance on timescales much shorter than overall cell-type changes. Under these conditions, the empirical count distribution closely mirrors the model’s stationary solution, justifying steady-state inference of burst size and frequency from a single time point. This assumption is commonly adopted in probabilistic models of transcriptional bursting (Larsson AJM et al., 2019; Raj A & van Oudenaarden A, 2008).

However, this steady-state assumption has some limitations. First, in some scenarios, the cell system may exhibit highly transient transcriptional programs that do not satisfy stationarity, leading to biased or misleading parameter estimates. For example, immediately following a synchronized developmental stimulus—such as serum shock–induced activation of immediate-early genes. Second, because DeepTX infers the mean burst frequency and size across the population, it cannot recover the underlying time-resolved dynamics or distinguish heterogeneous kinetic subpopulations. 

We have added a statement in the discussion to acknowledge these limitations and suggest future extensions—such as incorporating time-series measurements or latent pseudo time covariates—to address non-stationarity and recover temporal burst dynamics.

On page 3, "traditional telegraph model" is mentioned without any context. This model, and particularly the implications for the current work, might not be obvious to the reader. Take one or two sentences to give the reader context.

Thank the reviewer for this helpful comment. We acknowledge that the mention of the "traditional telegraph model" on page 3 may not be immediately clear to all readers. The traditional telegraph model is a mathematical framework commonly used to describe gene expression burst dynamics, in which genes stochastically switch between active (ON) and inactive (OFF) states, with exponentially distributed waiting times for state transitions. To provide the necessary context, we added a brief introduction to the traditional telegraph model and its relevance to our work in the revised manuscript.

A primary concern with the model used in Figure 2a (TXmodel) is its reliance on four independent parameters to describe gene state switching dynamics. While this general model can encompass specific cases such as the refractory model (Science 332, 472 (2011)) and the telegraph model, accurately estimating the parameters of the refractory model using only steady-state distributions and typical cell numbers (10³-10⁴) is challenging without time-dependent data. To address this, we suggest that the authors provide parameter inference results for each individual parameter, rather than only for burst size and burst frequency, based on synthetic data. This would help clarify the model's effectiveness and improve understanding of its estimation precision.

Thank the reviewer for highlighting this important concern. We agree that the lack of timeresolved measurements may affect the accuracy of inferences about dynamic parameters, especially the unidentifiability of parameters inferred from steady-state distributions, i.e., multiple parameters leading to the same steady-state distribution. The unidentifiability of individual parameters is a common and critical problem in systems biology studies. To address this issue, for example, Trzaskoma et al. developed StochasticGene, a computationally efficient software suite that uses Bayesian inference to analyze arbitrary gene regulatory models and quantify parameter uncertainty across diverse data types (Trzaskoma P et al., 2024). Alexander et al. adopt a Bayesian approach to parameter estimation by incorporating prior knowledge through a prior distribution and classify a parameter as practically nonidentifiable if it cannot be uniquely determined beyond the confidence already provided by the prior (Browning AP et al., 2020). Hence, in DeepTX, we employed a Bayesian approach based on loss potential to infer the posterior distributions of the parameters (Figure 3E). 

Although DeepTX also encounters the issue of unidentifiability for individual parameters (Supplementary Figure S11), the multimodal nature of the posterior distribution suggests that multiple distinct parameter sets can produce similarly good fits to the observed data, highlighting the inherent non-identifiability of the model. Nevertheless, in the multimodal posterior distribution, at least one of the posterior peaks aligns closely with the ground truth, thereby demonstrating the validity of the inferred result. Moreover, inference results on synthetic data confirm that the BS and BF can be accurately estimated (Supplementary Figure S3b and S3c). We also performed robustness analyses on synthetic datasets. As shown in Supplementary Figure S3d and S3e, our model reliably recovers the ground-truth burst kinetics of models when the number of cells reaches ~1000, which is within the range of typical single-cell RNA-seq experiments. 

We have explicitly pointed out the potential issue of unidentifiability due to the lack of temporal resolution information in the discussion section. 

Noteworthy, transcriptional is always a multi-step process (depending on the granularity with which the process is described). What do the authors mean by saying that "DNA damage turns transcription into a multi-step process rather than a single-step process"?

Thank the reviewer for pointing out the lack of precision in our original statement. We agree that the phrasing could be misleading. Transcription is inherently a multi-step process, but most mechanistic studies simplify it to a single-step “telegraph” model for tractability. In the context of DNA damage, however, damage-induced pausing and repair-mediated delays introduce additional intermediary states in the transcription cycle that cannot be approximated by a single step. To capture these damage-specific interruptions, DeepTX explicitly consider a multi-step promoter switching framework rather than combining all transitions into one. What we originally wanted to express was the necessity of multi-step process modeling. We have replaced the original sentence in introduction with: “However, the presence of DNA damage necessitates modeling the transcriptional process as a multistep process, rather than a single-step process, to capture the additional complexity introduced by the damage”.

It is unclear why the authors have chosen a different definition in Equation (2) rather than the commonly used burst frequency, 1/(k_deg * tau_off), as reported in the literature. Unlike the traditional definition, which is unit-free, the definition in Eq. (2) includes units, raising questions about its interpretability and consistency with established conventions. Clarifying this choice would improve the understanding and consistency of the methodology.

Thank the reviewer for raising this important point. We acknowledge that there are multiple definitions of burst frequency (BF) in the literature. Here, we provide a detailed explanation, clarifying the differences between these definitions, including the one used and the traditional definition .

First, the definition of burst frequency we adopt has been widely used in recent literatures, such as Benjamin Zoller et al. (Zoller B et al., 2018), Caroline Hoppe et al. (Hoppe C et al., 2020) and Daniel Ramsköld (Ramsköld D et al., 2024). And its quantity represents the average time it takes for the promoter to complete one full stochastic cycle between its active and inactive states . Secondly, the traditional definition can be regarded as a simplified version of our definition, under the assumptions that τ_on is negligible and k_deg =1 (i.e., rate parameters are normalized to be unit-free). Although it is reasonable to neglecting activate time τ_on, as it is typically much shorter than inactive time under some conditions, we chose a more complete way to define the burst frequency so that it is applicable to more general situations. In addition, by defining the burst frequency as , the mean transcription level can be analytically represented as the product of burst size and burst frequency.

This explanation has been clarified in the methods 1.2 section.

The authors mention the need to model "more realistic gene expression processes". How is this exactly being incorporated into the model?

Thank the reviewer for raising this important question. To incorporate "more realistic gene expression processes" into our model, we considered two critical aspects into DeepTX that are often oversimplified in traditional approaches:

(1) Integration of gene expression and sequencing processes: Observations from scRNA-seq data are influenced by both the intrinsic gene expression processes and the subsequent sequencing procedure. Traditional models often focus solely on gene expression, neglecting the stochastic effects introduced by the sequencing process. Our model explicitly incorporates both the gene expression and sequencing processes, providing a more comprehensive and realistic representation of the observed data.

(2) Modeling gene expression as a multi-step process: Gene expression is inherently a multi-step process. However, traditional telegraph models typically simplify gene state switching as a single-step process for tractable analysis, often assuming Markovian dynamics where transition waiting times follow exponential distributions. In contrast, our model accounts for the multi-step nature of gene state transitions by allowing the waiting times to follow non-exponential (non-Markovian) distributions. This model is more suitable for gene expression dynamics that cannot be simplified to a single-step process, such as DNA damage, which may introduce an intermediate state to represent pausing and repair in the transcription process.

By addressing these factors, our model better reflects the complexity and stochastic nature of gene expression processes, aligning more closely with the data generated from biological systems. We have added detailed explanations after this sentence for clarification in the first result section.

Better explanation of the previously developed TXmodel, and the assumption of a non-Markovian system. In particular, it isn't clear how using arbitrary distributions for the waiting times implies a non-Markovian process (as the previous state(s) of the system is not used to inform the transition probability, at least as explained in pp. 4). Without a clear discussion of the so-called arbitrary waiting time distribution, it isn't clear how these represent a mechanistic model. In general, a more careful discussion of the "mechanistic" model is needed.

Thank the reviewer for this thoughtful comment. In this revised version, we provided a more detailed explanation of the relationship between the TXmodel and the non-Markovian system in the revised manuscript. Specifically, we will clarify the following points:

(1) Why non-Markovian system: In a Markovian system, the waiting times for events are exponentially distributed, meaning that the state transitions depend solely on the current state and are memoryless (Van Kampen NG, 1992). However, when the waiting times follow non-exponential distributions, such as Gamma or Weibull distributions, the state transitions are no longer independent of the system's previous states. This introduces memory into the system, making it non-Markovian.

(2) Why mechanistic model: First, it is important to clarify that regardless of whether the waiting time is arbitrary or exponential (corresponding to non-Markovian and Markovian systems), our TXmodel is a mechanistic model because it models the dynamic process of transcription bursts with interpretable kinetic parameters. Second, although we introduced arbitrarily distributed waiting times, reasonable selection of waiting time distributions can still make the distribution parameters mechanistically interpretable. For example, in the context of modeling ON and OFF state switching times using a Gamma distribution, the two parameters have clear interpretations: the shape parameter represents the number of sequential exponential (memoryless) steps required for the transition to occur, capturing the complexity or multi-step nature of the switching process, while the scale parameter denotes the average duration of each of these steps. We have added the explanation in methods 1.2 section.

Include a brief discussion about the metric used to compare distributions (and introduce KL abbreviation).

Thank the reviewer for this suggestion. In the second result and methods 1.3 section of revised manuscript, we have included a brief discussion to introduce and clarify the metric used to compare distributions. Specifically, we have given more explanation for the Kullback-Leibler (KL) divergence, which is a widely used metric for quantifying the difference between two probability distributions. We also ensured that the abbreviation "KL" is properly introduced when it first appears in the text, along with a concise description of its mathematical definition and interpretation within the context of our analysis. 

What does the "CTM" model stand for (in supplementary information)? And "TX" model?

Thank the reviewer for highlighting this point. We revised the supplementary information to explicitly define the "CTM" and "TX" models and clarify their distinctions.

CTM model: The "CTM" model refers to the classic telegraph model, a widely used model for capturing Markovian gene expression burst kinetics. The CTM describes stochastic gene expression as a sequence of four biochemical reactions involving two gene states (ON and OFF), mRNA transcription and degradation:

k_off as the rate at which the gene switches from OFF to ON, k_on  as the rate at which the gene switches from ON to OFF, k_syn  as the rate of mRNA synthesis and k_deg  as the rate of mRNA degradation. In this model, gene switching between active and inactive states is governed by a memoryless Markovian process, where the waiting times for transitions follow exponential distributions (Van Kampen NG, 1992).

TX model: In contrast, the "TX" model is a more generalized telegraph model for transcriptional processes.

Different from the CTM, the waiting times for state transitions between ON and OFF in the TX model follow arbitrary waiting time distributions. This implies that the future state of the system depends not only on the current state but may also be influenced by its historical trajectories. Consequently, the TX model exhibits non-Markovian behavior. We have added more detailed description on these two models in section 1.1 of supplementary text.

Leaky transcription (in the OFF promoter state) is not considered. What would be the implications of its presence in the data?

Thank the reviewer for pointing out the potential role of leaky transcription in our analysis. We acknowledge that leaky transcription, occurring in the promoter OFF state, was not explicitly considered in our current model. Our decision to exclude it assumed that the leaky transcription rate is relatively small and its impact on the observed data is negligible. This assumption is consistent with previous studies that similarly disregard leaky transcription in gene expression modeling due to its minimal contribution to the overall dynamics (Larsson AJM et al., 2019).

However, we recognize that the leaky transcription should be considered, particularly in systems where the leaky rate is significant relative to the active transcription rate. In such cases, it may introduce additional variability to the observed expression levels or obscure the distinction between ON and OFF states. We have added relevant statements in the discussion section.

In the main text, the waiting time for state transitions is described by two parameters, while in the methods/supplementary information only one parameter is considered per distribution (without a clear discussion of the so-called "dwell time distributions").

Thank the reviewer for this comment. We recognize the need to clarify the discrepancy between the descriptions of waiting times in the main text and supplementary materials.

Dwell time distribution refers to the probability distribution of the time in which a gene remains in a particular transcriptional state (ON or OFF) before transitioning to the other state. While in Markovian models the dwell time follows an exponential distribution, more complex or non-Markovian regulatory mechanisms may give rise to Gamma, Weibull, or other non-exponential dwell time distributions.

In our model, we denote the dwell time distributions in the OFF and ON states by and , respectively, where w represents a vector of parameters characterizing the distribution, the dimensionality of which depends on the specific form of the distribution. For example, when an exponential distribution is assumed, w consists of a single rate parameter; in contrast, for distributions such as the Gamma or Weibull, w includes two parameters. In the main text, both and are modeled using Gamma distributions, whereas in the Supplementary Materials, we assume exponential distributions for both, resulting in a single-parameter representation. We have added relevant statements in the methods 1.2 section.

Related, but more general, across the manuscript there are problems with the consistency in terminology. This is especially problematic with the figures. It makes it incredibly hard to follow the work. Better integration of the information, and consistency with the terminology, would improve the understanding for the reader.

Thank the reviewer for the valuable feedback. To enhance clarity and readability, we have carefully revised the manuscript to ensure consistent terminology throughout the text and figures e.g., unifying terms such as "untreatment" and "control" under the consistent label "control"—across both the text and figures.

One of the four main assumptions behind the model is that "the solution of the model can be explained by a mixed negative binomial distribution". The logic and implications of this assumption need to be discussed in the paper. (Methods, pp.13.) All four assumptions need to be carefully argued in the paper. 

We appreciate the reviewer’s comment regarding the assumptions underlying our model. Here, we would like to clarify the rationale and implications of each assumption. 

Assumption 1 (The gene expression of cells was in a stationary distribution during sequencing.) has been extensively used in previous studies for the inference and modeling of scRNA-seq data, demonstrating effectiveness in capturing mRNA expression distributions and inferring underlying dynamic parameters (Larsson AJM et al., 2019; Luo S et al., 2023; Ramsköld D et al., 2024; Gupta A et al., 2022).

For Assumption 2 (Gene expression counts of the same cell type follow the same distribution.) is as follows: cell types are typically defined based on gene expression profiles or functional characteristics. Cells with similar functions often exhibit consistent transcriptional programs, leading to approximately identical gene expression distributions. This assumption has been widely adopted in previous research (Larsson AJM et al., 2019; Gupta A et al., 2022).

Regarding Assumption 3 (The solution of the model can be approximated by a mixed negative binomial distribution.), in the most general formulation, a chemical master equation (CME) model of biological systems converges to a stationary distribution P(n;θ) over n∈ℕ. And P(n;θ) afford a real Poisson representation (Gardiner CW & Chaturvedi S, 1977): where F is a mixing cumulative distribution function (CDF). If such a Poisson representation exists, we can always write down a finite approximation over K Poisson kernels: , where w_k are weights on a K-dimensional simplex. Further, as k →∞,Q→P . More problematically, convergence in the number of kernels in K is typically slow. Negative binomial kernels P_Poisson (n m_k,l_k), which are continuous Poisson mixtures with a gamma mixing density can accelerate convergence in K (Gorin G et al., 2024). Hence, the solution of the TX model can be approximated by a mixed negative binomial distribution. 

For Assumption 4 (The state space sampled from a sufficiently long single simulation is statistically equivalent to that obtained from multiple simulations at steady state in gene expression models.), when a sample trajectory of the model is simulated for a sufficiently long period, it is assumed to have traversed the entire stationary state space (Kuntz J et al., 2021). Therefore, by performing truncated statistical analysis on the trajectory, the corresponding stationary distribution of the model can be obtained. We have added the explanation in methods 1.1 section.

The authors propose that the waiting times between promoter states follow a non-exponential distribution, but the choice of gamma distribution and the implications for the method and the biological conclusions need to be discussed.

We thank the reviewer for this comment. To account for the impact of DNA damage on the transcription process, our model assumes that both the "ON" and "OFF" states of the promoter consist of multiple underlying sub-states. When a promoter switches from the "ON" state to the "OFF" state, the transition is governed by multiple distinct waiting time distributions that follow exponential distributions. Similarly, when a promoter switches from the "OFF" state to the "ON" state, there may be multiple transitions from different "OFF" sub-states. Consequently, the waiting times for the transitions from the "OFF" state to the "ON" state, and vice versa, must account for multiple exponential waiting time distributions associated with each "ON" state transition. We can map a multiple exponential-waiting-times reaction process to a single-step reaction process with a non-exponential waiting time distribution. Therefore, we use a Gamma distribution for dwell time of promoter switching, which can be expressed as the convolution of multiple exponential distributions (corresponding to a sum of multiple exponential variables). Additionally, other non-exponential distributions, such as those discussed in our previous studies (Zhang J & Zhou T, 2019), may also be considered, and we recognize that alternative choices could be made depending on the specific characteristics of the system. We have added the explanation in methods 1.2 section.

BF - burst frequency; BS - burst size. These terms represent the main data output, but they are only mathematically defined in the methods, and never the intuition of the specific expression explained (e.g., why not using tON/(tON+tOFF) as BF instead of 1/(tON+tOFF), and why not kSYN*tON as BS instead of kSYN*tON).

We appreciate the reviewer’s comment and agree that clarifying the biological intuition behind the mathematical definitions of burst frequency (BF) and burst size (BS) is important. Below, we provide a more detailed explanation of these definitions.

BF: The definition of burst frequency we adopt has been widely used in previous literature, such as Benjamin Zoller et al (Zoller B et al., 2018), Caroline Hoppe et al (Hoppe C et al., 2020) and Daniel Ramsköld (Ramsköld D et al., 2024). And its quantity represents the average time it takes for the promoter to complete one full stochastic cycle between its active and inactive states.

BS: The definition of burst size BS = we adopt is consistent with the definition proposed by the reviewer. Burst size refers to the average number of mRNA transcripts produced during a single transcriptional activation event of a gene. It reflects the quantity of gene product synthesized per activation and is influenced by the rate of transcription and the duration of the active state of the gene. Our definition aligns with this biological interpretation and is mathematically formulated as BS = , where k_syn is the transcription rate and is the average duration of the active state.

In addition, the mean transcription level can be analytically represented as the product of burst size and burst frequency. This analytical result has been included in the methods 1.2 section of revised manuscript.

One can assume from the methods that omegaON and omegaOFF are the vector of (2) parameters describing the distribution, but the reader would benefit from some clarity here. The authors claim that they proved that the distribution moments can be obtained through an iterative process. How much does this rely on the assumption of an underlying binomial distribution?

Thank the reviewer for this helpful suggestion. To clarify, the vectors omegaON and omegaOFF represent the parameters characterizing the waiting time distributions of the promoter's active and inactive states, respectively. The exact form and interpretation of these vectors depend on the specific distributional choice for the waiting times. For instance, when the waiting time distribution follows a Gamma distribution with shape parameter α>0  and scale parameter β>0 , denoted as , then w_on = (α,β) . Conversely, when the waiting time distribution follows a Weibull distribution, denoted as , with shape parameter k >0 and scale parameter l>0, then w_on = (l,k) . We have clarified it in the Methods 1.2 section of the revised manuscript.

For the question about the binomial distribution, in our work, we use the binomial moment method to compute distributional statistics of chemical master equation (Zhang J et al., 2016). Binomial moments of the mRNA stationary distribution P(m) are defined as , where the symbol represents the combinatorial number. This technique refers to a mathematical tool for moment calculation and is not based on the assumption that the underlying distribution is binomial distribution (Luo S et al., 2023). Hence, our approach is general and does not require the distribution itself to follow a binomial form.

More details about the parameter sampling are required. For instance, why are the specific ranges chosen and their implications? And is the space explored in logarithmic scale?  

Thank the reviewer for the insightful comment regarding parameter sampling. In our study, we considered five parameters: . The parameters k_off  and k_on represent the number of intermediate reaction steps involved in transcriptional state transitions. These values were sampled uniformly from the range 1 to 15, which aligns with biological evidence indicating that most genes undergo either direct (single-step) transitions or a small number of intermediate steps, typically fewer than ten (Tunnacliffe E & Chubb JR, 2020). This range is sufficient to capture both widely used singlestep models and more detailed multi-step mechanisms without introducing biologically implausible complexity. 

Among these parameters, r_off and r_on denote the rate constants governing stochastic transitions between the OFF and ON transcriptional states, respectively. The mean duration of the OFF state, which corresponds to the time between transcriptional bursts, is given by = k_off / r_off , and falls within the range ∈(0.1,150).Experimental measurements report a median value of approximately 3.7 (Gupta A et al., 2022), which is well contained within this range. Similarly, the mean duration of the ON state, referred to as the burst duration, is defined by = k_on / r_on , and spans the interval ∈(0.1,1500). The experimentally observed median value of 0.12 (Gupta A et al., 2022) confirms that the parameter range adequately captures biologically realistic dynamics.

The parameter k_syn  represents the normalized synthesis rate after accounting for molecular degradation. Its range was chosen based on empirical observations of transcriptional burst sizes, which typically vary from single molecules to several dozen (Gupta A et al., 2022). Considering the relationship BS = k_syn * , the selected range of k_syn ensures that the experimentally observed burst sizes are well represented within the defined parameter space. We have added the explanation in methods 1.2 section and supplementary text 4.

We fully recognize the advantages of logarithmic sampling, particularly when parameters span several orders of magnitude. Logarithmic scaling ensures balanced exploration across wide ranges and prevents sampling bias towards larger values. However, in our work, we applied Sobol sampling directly within the original (linear) parameter space. Although we did not explicitly transform parameters into logarithmic scale, Sobol sequences provide low-discrepancy, quasi-random coverage, which promotes uniform sampling across bounded domains (Sobol IM, 1967). Further, if necessary, we can increase the parameter range adaptively, and perform simulation algorithm to obtain sample and train a new model to solve a larger parameter range. 

On page 15, the rationale for selecting the parameter space is unclear. This is crucial, as fully connected neural networks typically exhibit poor extrapolation beyond their training parameter space. If the parameter space of an experimental dataset significantly differs from the training range, the inference results may become unreliable. We suggest further clarification on how the alignment between the parameter spaces of the experimental data and the training dataset can be ensured to maintain inference accuracy.

We appreciate the reviewer’s insightful comment regarding the extrapolation limitations of fully connected neural networks. To address this concern, we have implemented a truncation strategy during inference, which constrains the inferred parameters to remain within the bounds of the training parameter space. This ensures that the neural network operates within a regime where its predictive accuracy has been validated, thereby enhancing the robustness of our results. Additionally, we have carefully selected the training parameter space to be reasonable, based on the characteristics of the experimental data. These ranges have been validated through domain knowledge and data analysis, ensuring that even when the experimental data approaches the boundaries of the training range, the inference results remain reliable and accurate.

On page 16, it is unclear why the authors chose to incorporate the Fano factor instead of using the coefficient of variation or variance. Clarifying the reasoning behind the selection of the Fano factor over these other statistical measures would provide better insight into its relevance for their analysis.  

We thank the reviewer for raising this point. Although the loss term is described using the Fano factor, its formulation actually involves both the variance and the mean. Specifically, the loss we use is: . We chose to use the Fano factor because it is particularly well-suited for quantifying transcriptional noise in systems where the mean expression level varies across conditions or parameters. Unlike variance, the Fano factor normalizes variability by the mean, making it more robust for comparing noise levels across genes or regulatory regimes with different expression levels. Compared to the coefficient of variation (CV), which normalizes by the square of the mean, the Fano factor tends to be less sensitive to low expression regimes and is commonly used in stochastic gene expression studies, especially when the distribution is skewed or over dispersed (i.e., variance exceeds the mean). This makes it a more appropriate metric in our context, where transcriptional bursting often leads to over dispersed expression distributions. We have added an explanation in the methods 1.3 of revised manuscript to explain this choice.

On page 17, the definition of "sample" is unclear. Does it refer to the number of parameters sets or to the simulated trajectories generated by stochastic simulation algorithms?

Thank reviewers for your valuable feedback. The term "sample" in this context refers to the data points used in the neural network training set. To eliminate any ambiguity, we included a precise mathematical definition of "sample" (θ_i,P_simulation,i ) in the methods 1.3 section of revised manuscript.

Additionally, it is unclear how the authors determined the number of simulated trajectories per parameter set to ensure training accuracy. Furthermore, it would be relevant to address whether including moments during neural network training is beneficial.

We appreciate the reviewer’s insightful questions regarding the simulation and training process. To clarify, for each parameter set, we did not simulate multiple trajectories to obtain the corresponding distribution. Instead, we simulated the system for a sufficiently long period to ensure that the system reached a steady-state distribution. From this steady-state data, we then used interpolation methods to derive the corresponding distribution for each parameter set.

On the other hand, the moments were calculated theoretically without any approximations, providing higher accuracy. By incorporating the moments into the training process, we can effectively mitigate potential biases arising from insufficient sampling of the simulated data. Moreover, our experiments on the synthetic dataset demonstrate that introducing the moments as a loss function significantly enhances the model's performance on the test set (Figure 2E).

What is the intuition behind the choice of alpha_cg? On page 18, the rationale for setting the sampling probability to 0.5 is unclear. Could this parameter be inferred rather than being preset?  

We thank the reviewer for the insightful comment regarding the choice of α_cg. We acknowledge that the typical values of this parameter in related literature often fall within a narrower range (e.g., 0.06–0.32) (Zheng GX et al., 2017; Macosko EZ et al., 2015). However, our decision to set α_cg was based on a trade-off between sampling efficiency and computational tractability in our specific application context. While it is indeed possible to infer α_cg as a learnable parameter, we opted for a fixed value in this work to reduce model complexity and avoid unidentifiability issues. In addition, we conducted inference under different capture efficiencies (0.5, 0.3, and 0.2), and found that the inferred burst size (BS) and burst frequency (BF) remained strongly correlated across these conditions (Supplementary Figure S12). This indicates that variations in capture efficiency do not significantly impact the outcomes of downstream enrichment analyses. Nevertheless, we agree that adaptively learning α_cg could be a promising direction, and we plan to explore this in future work. We have added the explanation in methods 1.4 section.

On page 19, the authors employed gradient descent for parameter inference. However, as this method is sensitive to initial values, it is unclear how the starting points were selected.

We sincerely thank the reviewer for highlighting the sensitivity of gradient-based optimization methods to initial values. To address this concern, we adopted a black-box optimization strategy in the form of the adaptive differential evolution (DE) algorithm (Das S & Suganthan PN, 2010) to derive robust initial parameters for the parameter inference. The adaptive DE algorithm enables global exploration across a broad parameter space, thereby reducing the risk of convergence to suboptimal local minima. This yielded reasonably good initial estimates, which were subsequently refined using gradient-based optimization to identify high-quality solutions characterized by a vanishing gradient norm. This hybrid strategy, which combines global and local search, is widely adopted in optimization literature to alleviate the risk of entrapment in local optima (Ahandani MA et al., 2014). We have clarified this detail in the third result of the revised manuscript.

Furthermore, clarification on how the gradients were computed - whether through finite difference approximation or other methods - would offer additional insight into the robustness and accuracy of their approach.

Thank reviewers for valuable feedback. Regarding the computation of gradients, we use the chain rule in neural networks, and the gradients are computed through backpropagation. Specifically, we rely on automatic differentiation to efficiently calculate the gradients. Unlike finite difference approximation, automatic differentiation directly computes the derivative of the loss function with respect to each parameter, ensuring accurate gradient calculations (Baydin AG et al., 2018). We have clarified this detail in the discussion section of the revised manuscript.

The paper presents several comparisons between continuous and discrete distributions in Figure 2B and Supplementary Figures S4, S6, and S8, described as a "comparison between mRNA distribution and inferred distribution by DeepTX for scRNA-seq data" or a "comparison between SSA results and DeepTX prediction results." This may lead to confusion for the reader, as the paper focuses on transcriptional bursting, a process where we would typically expect the distributions to be discrete. Clarifying this point would help align the figures with the main topic and enhance the reader's understanding.

We sincerely thank the reviewer for this insightful comment. We understand the concern that the distributions shown in Figure 2B and Supplementary Figures S4, S6, and S8 may appear to be continuous, which could be confusing given that transcriptional bursting naturally results in discrete mRNA count distributions.

We have clarified that in all these figures, both the empirical mRNA distributions derived from scRNAseq data and the model-predicted distributions from DeepTX are inherently discrete. To visualize the empirical distributions, we used histograms where the x-axis corresponds to discrete mRNA copy numbers and the y-axis represents the normalized frequency (density). To illustrate the DeepTX-inferred probability mass function, we plotted the predicted probabilities at each integer count as points and connected them with lines for clarity. While the connecting lines give the appearance of continuity, this is a standard graphical convention used to better show trends and model fit in discrete distributions.

We suggest that Figure 3E could present the error as a percentage of the parameter value, as this would provide a more equitable comparison and better illustrate the relative accuracy of the parameter estimation.

Thank reviewers for suggestion regarding Figure 3E. We agree that presenting the error as a percentage of the parameter value would offer a more equitable basis for comparison and better highlight the relative accuracy of our parameter estimation. Accordingly, we have revised Figure 3E to include the relative percentage error for each parameter.

Figure 4A could be improved for better legibility. The contour plots are somewhat confusing, and the light blue points are difficult to distinguish. Additionally, the x-axis label "Untreatment" appears throughout the manuscript-could this term be referring to the control experiment?

Thank reviewers for constructive feedback. We have revised Figure 4A to improve its clarity and legibility. Specifically, we adjusted the display style of the contour plots and enhanced the visibility of the light green points to make them more distinguishable.

Additionally, we recognize the potential confusion caused by the term "Untreatment" and have replaced it with "Control" throughout the revised manuscript to ensure consistency and accuracy in terminology.

Figure 4B was unclear, and further explanation would be helpful for understanding its purpose.

Thank reviewers for feedback. The purpose of Figure 4B is to illustrate the relationship between bursting kinetics and the mean and variance of the model. In the revised manuscript, we will provide a more detailed explanation of how the figure captures these relationships, highlighting the key insights it offers into the underlying dynamics.

Figure 4B illustrates the quantitative relationships among BS, BF, and gene expression noise within the framework of the transcriptional model. In this log-log-log 3D space, the mean expression level is constrained on a blue plane defined by the equation log(BS)+log(BF) = log(Mean), highlighting that the product of burst size and burst frequency determines the mean expression level. The orange plane represents a scaling relationship between expression noise and burst kinetics, expressed as log(BS)+log(BF) = klog(Noise), where k is a constant indicating how the burst kinetics co-vary with noise. Notably, the trajectory of the green sphere demonstrates that, under a fixed mean expression level (i.e., remaining on the blue plane), an increase in gene expression noise arises primarily from an increase in burst size. We have revised the caption of Figure 4B.

In Figure 4D, two of the GO analysis terms are highlighted in red, but the meaning behind this emphasis is not clear. The same question applies to Figure 5E, where the green dots are missing from the plot.

Clarification on these points would enhance the overall clarity.  

We appreciate the reviewer’s thoughtful comments. We have added further clarification regarding the enrichment analysis results presented in Figure 4D. Specifically, we highlighted the "cell cycle G2/M phase transition" pathway because a delay in the G2/M phase transition has been shown to increase the probability of cell differentiation, which is a key aspect of our study. In addition, since IdU treatment is known to induce DNA damage, we emphasized the DNA damage-related pathway to support the biological relevance and consistency of our enrichment results. Similarly, in Figure 5E, we highlighted the apoptosis-related pathway. Apoptosis in this context is closely associated with cellular responses to toxic substances and mitochondrial dynamics. The enrichment of pathways related to these processes enables us to hypothesize the underlying mechanisms driving apoptosis in our system. Further, the absence of green dots in Figure 5E was due to an error in the figure caption. We have revised the figure caption accordingly to accurately describe all elements presented in the figure.

Clarify axis labels in figures, particularly the y-axis in Figure 5A and the x-axis in Figure 6G. In the first case, it isn't clear what this "value" represents. In the second case, the x-label is very confusing. As I understand the figure description, in these plots you are always comparing the G0 arrested genes between control and treated cells. But the x-label says "G0 (0 D)", "Cycle (50 D)".

Thank reviewers for pointing out the issues with the axis labels. We have made the necessary revisions to eliminate any confusion. In Figure 5A, the label for the y-axis has been changed from "value" to "log2 (value)" for clarity. The “value” in y-axis represents the value of statistical measure indicated at top of each panel. In Figure 6G, the x-axis label "Cycle (50 D)" has been updated to "G0 (50 D)" to accurately reflect the comparison between the G0-arrested genes in control and treated cells. We have revised the text of Figure 5A and Figure 6G.

Figure 6 uses a QS metric (quality score), but the definition of this metric is not provided. Including a brief explanation of its meaning would be helpful for clarity.  

Thank reviewers for feedback. In this version, we provided explanation of the QS (Quality Score) metric in the supplementary text 3 for better clarity. The QS is calculated based on the difference in z-scores derived from GSVA (Gene set variation analysis) of gene sets upregulated and downregulated during the quiescent phase, and is defined as QS = z(up genes)− z(down genes) , as described in the literature (Wiecek AJ et al., 2023). z(up genes) represents the standardized enrichment score of the gene set upregulated during quiescence in each sample. A higher value indicates that the quiescenceassociated upregulated genes are actively expressed, suggesting that the sample is more likely to be in a quiescent (G0) state. z(down genes)  corresponds to the standardized enrichment score of genes downregulated during quiescence. A lower value implies effective suppression of these genes, which is also consistent with quiescence. The difference score QS serves as an integrated indicator of the quiescent state: A higher value reflects simultaneous activation of quiescence-associated upregulated genes and repression of downregulated genes, indicating a gene expression profile that strongly aligns with the G0/quiescent state. A lower or negative value suggests a deviation from the quiescent signature, potentially reflecting a proliferative state or failure to enter quiescence. 

In Figure 6G, light grey lines are shown, but their significance is unclear. It would be useful to specify what these lines represent.

Thank reviewers for observation. In Figure 6G, each point represents a single gene, and the light grey lines indicate the trend of changes in the corresponding bursting kinetics values, mean and variance for genes. We have added the explanation in the caption of Figure 6G.

Additionally, the manuscript should include references to the specific pathways used in the GO analysis to provide more context for the reader.

Thank reviewers for the suggestion. We have included references to the specific pathways used in the GO analysis in the revised manuscript to provide additional context for the readers.

In the discussion, sentences like "IdU drug treatment-induced BS enhancement delays the cell mitosis phase transition, impacting cell reprogramming and differentiation" are problematic as they imply causality, which I believe cannot be determined through the present analysis. The strength of the conclusions needs to be better argued (or toned down).

We acknowledge that the original sentence lacked precision and may have conveyed a misleading implication of causality not fully supported by our current analysis. In the discussion section of revised manuscript, we have rephrased the statement to present a more nuanced interpretation: IdU drug treatment-induced BS enhancement of genes may be associated with a delayed transition in the cell mitosis phase, which could potentially influence cell reprogramming and differentiation.  

Other (minor) comments:

On pp. 10, "the BS down-regulates differential genes were mainly enriched..." appears to have a grammatical error/typo, "down-regulated"?

We have made correction. We have revised “down-regulates” to “down-regulated” for grammatical consistency.

Equation 2 doesn't match Figure 1A.

We have made correction. The definition of BF = in Equation 2 is correct. We have revised the definition of BF in Figure 1A to ensure consistency with Equation 2.

Reference

Zheng, G.X., Terry, J.M., Belgrader, P., Ryvkin, P., Bent, Z.W., Wilson, R., Ziraldo, S.B., Wheeler, T.D., McDermott, G.P., Zhu, J., Gregory, M.T., Shuga, J., Montesclaros, L., Underwood, J.G., Masquelier, D.A., Nishimura, S.Y., Schnall-Levin, M., Wyatt, P.W., Hindson, C.M., Bharadwaj, R., Wong, A., Ness, K.D., Beppu, L.W., Deeg, H.J., McFarland, C., Loeb, K.R., Valente, W.J., Ericson, N.G., Stevens, E.A., Radich, J.P., Mikkelsen, T.S., Hindson, B.J., Bielas, J.H. 2017. Massively parallel digital transcriptional profiling of single cells. Nature Communications 8: 14049. DOI: https://dx.doi.org/10.1038/ncomms14049, PMID: 28091601

Hagemann-Jensen, M., Ziegenhain, C., Chen, P., Ramsköld, D., Hendriks, G.J., Larsson, A.J.M., Faridani, O.R., Sandberg, R. 2020. Single-cell RNA counting at allele and isoform resolution using Smart-seq3. Nature Biotechnology 38: 708714. DOI: https://dx.doi.org/10.1038/s41587-020-0497-0, PMID: 32518404

Larsson, A.J.M., Johnsson, P., Hagemann-Jensen, M., Hartmanis, L., Faridani, O.R., Reinius, B., Segerstolpe, A., Rivera, C.M., Ren, B., Sandberg, R. 2019. Genomic encoding of transcriptional burst kinetics. Nature 565: 251-254. DOI: https://dx.doi.org/10.1038/s41586-018-0836-1, PMID: 30602787

Ochiai, H., Hayashi, T., Umeda, M., Yoshimura, M., Harada, A., Shimizu, Y., Nakano, K., Saitoh, N., Liu, Z., Yamamoto, T., Okamura, T., Ohkawa, Y., Kimura, H., Nikaido, I. 2020. Genome-wide kinetic properties of transcriptional bursting in mouse embryonic stem cells. Science Advances 6: eaaz6699. DOI: https://dx.doi.org/10.1126/sciadv.aaz6699, PMID: 32596448

Luo, S., Wang, Z., Zhang, Z., Zhou, T., Zhang, J. 2023. Genome-wide inference reveals that feedback regulations constrain promoter-dependent transcriptional burst kinetics. Nucleic Acids Research 51: 68-83. DOI: https://dx.doi.org/10.1093/nar/gkac1204, PMID: 36583343

Rodriguez, J., Ren, G., Day, C.R., Zhao, K., Chow, C.C., Larson, D.R. 2019. Intrinsic dynamics of a human gene reveal the basis of expression heterogeneity. Cell 176: 213-226.e218. DOI: https://dx.doi.org/10.1016/j.cell.2018.11.026, PMID: 30554876

Luo, S., Zhang, Z., Wang, Z., Yang, X., Chen, X., Zhou, T., Zhang, J. 2023. Inferring transcriptional bursting kinetics from single-cell snapshot data using a generalized telegraph model. Royal Society Open Science 10: 221057. DOI: https://dx.doi.org/10.1098/rsos.221057, PMID: 37035293

Eling, N., Morgan, M.D., Marioni, J.C. 2019. Challenges in measuring and understanding biological noise. Nature Reviews Genetics 20: 536-548. DOI: https://dx.doi.org/10.1038/s41576-019-0130-6, PMID: 31114032

Tunnacliffe, E., Chubb, J.R. 2020. What is a transcriptional burst? Trends in Genetics 36: 288-297. DOI: https://dx.doi.org/10.1016/j.tig.2020.01.003, PMID: 32035656

Rodriguez, J., Larson, D.R. 2020. Transcription in living Cells: molecular mechanisms of bursting. Annual Review of Biochemistry 89: 189-212. DOI: https://dx.doi.org/10.1146/annurev-biochem-011520-105250, PMID: 32208766

Morgan, M.D., Marioni, J.C. 2018. CpG island composition differences are a source of gene expression noise indicative of promoter responsiveness. Genome Biology 19: 81. DOI: https://dx.doi.org/10.1186/s13059-018-1461-x, PMID: 29945659

Raj, A., van Oudenaarden, A. 2008. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell 135: 216-226. DOI: https://dx.doi.org/10.1016/j.cell.2008.09.050, PMID: 18957198

Trzaskoma, P., Jung, S., Pękowska, A., Bohrer, C.H., Wang, X., Naz, F., Dell’Orso, S., Dubois, W.D., Olivera, A., Vartak, S.V. 2024. 3D chromatin architecture, BRD4, and Mediator have distinct roles in regulating genome-wide transcriptional bursting and gene network. Science Advances 10: eadl4893. DOI: https://dx.doi.org/https://www.science.org/doi/10.1126/sciadv.adl4893, PMID: 

Browning, A.P., Warne, D.J., Burrage, K., Baker, R.E., Simpson, M.J. 2020. Identifiability analysis for stochastic differential equation models in systems biology. Journal of the Royal Society Interface 17: 20200652. DOI: https://dx.doi.org/10.1098/rsif.2020.0652, PMID: 33323054

Zoller, B., Little, S.C., Gregor, T. 2018. Diverse spatial expression patterns emerge from unified kinetics of transcriptional bursting. Cell 175: 835-847.e825. DOI: https://dx.doi.org/10.1016/j.cell.2018.09.056, PMID: 30340044

Hoppe, C., Bowles, J.R., Minchington, T.G., Sutcliffe, C., Upadhyai, P., Rattray, M., Ashe, H.L. 2020. Modulation of the promoter activation rate dictates the transcriptional response to graded BMP signaling levels in the drosophila embryo. Dev Cell 54: 727-741.e727. DOI: https://dx.doi.org/10.1016/j.devcel.2020.07.007, PMID: 32758422

Ramsköld, D., Hendriks, G.J., Larsson, A.J.M., Mayr, J.V., Ziegenhain, C., Hagemann-Jensen, M., Hartmanis, L., Sandberg, R. 2024. Single-cell new RNA sequencing reveals principles of transcription at the resolution of individual bursts. Nature Cell Biology 26: 1725-1733. DOI: https://dx.doi.org/10.1038/s41556-024-01486-9, PMID: 39198695 Van Kampen, N.G. 1992. Stochastic Processes in Physics and Chemistry. Elsevier.

Gupta, A., Martin-Rufino, J.D., Jones, T.R., Subramanian, V., Qiu, X., Grody, E.I., Bloemendal, A., Weng, C., Niu, S.Y., Min, K.H., Mehta, A., Zhang, K., Siraj, L., Al' Khafaji, A., Sankaran, V.G., Raychaudhuri, S., Cleary, B., Grossman, S., Lander, E.S. 2022. Inferring gene regulation from stochastic transcriptional variation across single cells at steady state. Proceedings of the National Academy of Sciences 119: e2207392119. DOI: https://dx.doi.org/10.1073/pnas.2207392119, PMID: 35969771

Gardiner, C.W., Chaturvedi, S. 1977. The Poisson representation. I. A new technique for chemical master equations. Journal of Statistical Physics 17: 429-468. DOI: https://dx.doi.org/https://doi.org/10.1007/BF01014349, PMID: 

Gorin, G., Carilli, M., Chari, T., Pachter, L. 2024. Spectral neural approximations for models of transcriptional dynamics. Biophysical Journal 123: 2892-2901. DOI: https://dx.doi.org/10.1016/j.bpj.2024.04.034, PMID: 38715358

Kuntz, J., Thomas, P., Stan, G.-B., Barahona, M. 2021. Stationary distributions of continuous-time Markov chains: a review of theory and truncation-based approximations. SIAM Review 63: 3-64. DOI: 

Zhang, J., Zhou, T. 2019. Computation of stationary distributions in stochastic models of cellular processes with molecular memory. bioRxiv: 521575. DOI: https://dx.doi.org/https://doi.org/10.1101/521575, PMID: 

Zhang, J., Nie, Q., Zhou, T. 2016. A moment-convergence method for stochastic analysis of biochemical reaction networks. The Journal of chemical physics 144. DOI: 

Sobol, I.M. 1967. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput. Math. Math. Phys. 7: 784-802. DOI: https://dx.doi.org/10.1016/0041-5553(67)90144-9, PMID: 

Macosko, E.Z., Basu, A., Satija, R., Nemesh, J., Shekhar, K., Goldman, M., Tirosh, I., Bialas, A.R., Kamitaki, N., Martersteck, E.M., Trombetta, J.J., Weitz, D.A., Sanes, J.R., Shalek, A.K., Regev, A., McCarroll, S.A. 2015. Highly parallel genome-wide expression profiling of individual cells using nanoliter dsroplets. Cell 161: 1202-1214. DOI: https://dx.doi.org/10.1016/j.cell.2015.05.002, PMID: 26000488

Das, S., Suganthan, P.N. 2010. Differential evolution: A survey of the state-of-the-art. IEEE transactions on evolutionary computation 15: 4-31. DOI: https://dx.doi.org/10.1109/TEVC.2010.2059031, PMID: 

Ahandani, M.A., Vakil-Baghmisheh, M.-T., Talebi, M. 2014. Hybridizing local search algorithms for global optimization. Computational Optimization and Applications 59: 725-748. DOI: https://dx.doi.org/https://doi.org/10.1007/s10589014-9652-1, PMID: 

Baydin, A.G., Pearlmutter, B.A., Radul, A.A., Siskind, J.M. 2018. Automatic differentiation in machine learning: a survey. Journal of machine learning research 18: 1-43. DOI: https://dx.doi.org/https://dl.acm.org/doi/abs/10.5555/3122009.3242010, PMID: 

Wiecek, A.J., Cutty, S.J., Kornai, D., Parreno-Centeno, M., Gourmet, L.E., Tagliazucchi, G.M., Jacobson, D.H., Zhang, P., Xiong, L., Bond, G.L., Barr, A.R., Secrier, M. 2023. Genomic hallmarks and therapeutic implications of G0 cell cycle arrest in cancer. Genome Biology 24: 128. DOI: https://dx.doi.org/10.1186/s13059-023-02963-4, PMID: 37221612

To me this is what people use to visualize things that’s sometimes better left unsaid but rather be viewed upon.

To some, art might even seem as if it’s simple but in some cases art makes u wonder/curious.

mucky purple and bets the cloth to clean him.

7<br /> 00:04:43

The excerpt "mucky purple and bets the cloth to clean him" likely describes a messy situation involving a purple substance and an attempt to tidy up. The use of "mucky" implies dirtiness or uncleanliness, with "purple" possibly indicating the color of the mess. The phrase "bets the cloth to clean him" suggests an action taken with a cloth, perhaps involving a gamble or risk in effectively cleaning up.

在这段摘录中，“mucky purple”和“bets the cloth to clean him”可能描述了一种令人烦恼的情况，涉及一种紫色的污垢和清理的尝试。“mucky”意味着肮脏，而“purple”可能暗示了污物的颜色。“bets the cloth to clean him”暗示用布进行某种清理的行动，似乎是在冒险能否有效清理干净。

Camron Newcomb

CC BY-NC-SA 4.0

The figure of Phineus, a blind prophet saved by male heroes (the Boreades) and rewarded with the masculine coded power of foresight reveals the way heroism is constructed through patriarchal intervention. In both Apollonius and the dramatized version by Aeschylus, Phineus is repositioned from a victim to a hero via male deliverance. Notably, the Harpies female monsters, represent chaos and disruption that must be tamed by male force, reinforcing traditional gender binaries where femininity is aligned with disorder, and masculine action with order and divine favor.

In this context, The Persians uses Phineus as a mythic prologue to set up Xerxes's downfall as a failure to embody the virtues of Hellenic masculinity: discipline, moderation, and obedience to divine will. Aeschylus contrasts the heroic male ideal of prophecy (Phineus, Darius) with the failed heroism of Xerxes, whose excessive ambition marks a deviation from the masculine ideal and leads to ruin.

Comparing this version of The Persians (as interpreted through 19th century scholarship like Plumptre’s) with Robert Auletta’s modern adaptation (1993) shows how gender is reframed over time. While Plumptre’s translation emphasizes stoic, hierarchical masculinity in line with Victorian values, Auletta’s contemporary version inserts more emotional vulnerability into Xerxes, complicating the classical heroic ideal. This shows how gender expectations shift with culture and time, revealing translation as an act of ideological transmission, not just linguistic rendering.

Camron Newcomb CC BY-NC-SA 4.0

This lyrical passage from The Persians offers a striking contrast between heroic masculinity and feminine suffering, deeply encoded in the gender politics of ancient Greek tragedy. The women are defined not by their own actions but by the absence of their men, reinforcing a binary where male heroism exists on the battlefield while female identity is rooted in passive emotional endurance.

The imagery “widow’d couch,” “pensive breast,” and “love lorn woe” frames these women as emotional vessels, symbolically tethered to the physical and martial exertions of men. Their suffering is romanticized and gendered, grief is feminized, domestic, and private, while heroism is masculinized, public, and glorified. This pattern reveals how female subjectivity is subordinated to the narrative arc of the male hero, echoing patriarchal ideologies.

From a linguistic standpoint, the poetic diction emphasizes emotional melodrama and uses bodily metaphors ("throb," "streaming tear") to anchor femininity in physical vulnerability. In contrast, men are described earlier in the text through martial ornamentation: “blazing with gold,” “proud steeds,” “massy spears,” etc. The translation here (Robert Potter’s 1777 version) clearly reflects the 18th century lens, romanticizing grief in highly gendered Victorian prose, potentially amplifying the patriarchal dimensions more than Aeschylus himself might have done in the original Greek.

Comparatively, this portrayal of women mirrors Sita’s position in The Ramayana and Soudabeh’s emotional manipulation in Shahnameh. In both cases, women are symbols of honor, temptation, or mourning, rather than autonomous actors. Meanwhile, male heroes like Rama, Siavash, and Beowulf embody courage through sacrifice and public duty, reaffirming a cultural pattern that links masculinity to action and femininity to reaction.

This annotation demonstrates how epic and dramatic literature across cultures constructs gendered heroism by emotionally loading female grief and idealizing male war-making. The juxtaposition deepens our understanding of how literary canon preserves, and sometimes critiques patriarchal hero myths.

PMID: 37288251

Gene: SCN5A

HGNC ID: 10593

Case#: 55-year-old man

DiseaseAssertion: single coronary artery (SCA) and presented with dilated cardiomyopathy (DCM)

FamilyInfo: Unremarkable

ParentalTesting: NR

CasePresentingHPOs: HP:0002094, HP:0031352, HP:0001638, HP:0001644, HP:0010741

CaseHPOFreeText: chest tightness and dyspnoea after activity lasting for 2 months. CTCA showed congenital absence of the right coronary artery. TTE revealed enlargement of the left heart and cardiomyopathy. CMR revealed DCM. oedema of both lower limbs. Laboratory data in Table 1.

CaseNotHPOs: NR

CaseNotHPOFreeText: Stenosis

CasePreviousTesting: See NGS results in Supplementary Table 1

Genotyping Method: Genetic screening (NGS results in Supplementary Table 1) with confirmation by Sanger

FunctionalAnalysis: NR

Variant: c.1858C>T (p.Arg620Cys)

ClinVar: 67694

CAID: CA015449

gnomAD: v4.1.0 GrpMax FAF: 0.00002033 (European non-Finnish)

AdditionalInfo: The patient also has APOA5:c.990_993delAACA (p. Asp332Valfs*5) (P/LP in ClinVar with 2 stars)

Camron Newcomb

CC BY-NC-SA 4.0

In this scene, the Danish coast guard stops Beowulf’s ship and immediately identifies one of the Geats (Beowulf himself) as an extraordinary figure: “Never a greater one / Of earls o’er the earth have I had a sight of.” This response not only reflects the cultural idealization of the hero’s physical appearance, but also reinforces how masculinity is visually constructed and recognized in warrior societies. The coast guard reads Beowulf’s armor, stature, and composure as clear signs of high status and heroic capability, connecting external form with internal worth, a hallmark of the gendered construction of the hero in epic literature.

The linguistic emphasis on “hero in armor,” “low ranking fellow,” and “striking appearance” shows that visual markers of masculinity, armor, weapons, height, posture are treated as symbolic credentials, establishing heroic identity before action even begins. This reveals a form of performative masculinity, where being seen as a man and a hero is almost as important as actually acting like one. The narrative rewards the ability to appear heroic even before deeds confirm it.

From a gender politics standpoint, this reinforces a patriarchal worldview where male bodies are not only expected to perform heroism but also to embody it visually a privilege and pressure that aligns with martial and aristocratic ideals of masculinity. Female figures, by contrast, are often rendered invisible or are not physically described unless tied to beauty or emotional traits, emphasizing how gender roles are linguistically and culturally encoded in unequal ways.

Comparatively, this construction mirrors figures like Rama in the Ramayana, whose beauty and bearing identify him as dharmic, or Siavash in Shahnameh, whose dignity and divine aura precede his moral trials. Each reinforces how masculinity and the heroic ideal are visually coded across traditions, reflecting the shared patriarchal values of ancient epic literature.

Camron Newcomb CC BY-NC-SA 4.0

This early commentary on Beowulf draws attention to the composite nature of its protagonist, a "mixed Nature" figure blending myth and historical legend. Importantly, it hints at the constructed ideal of the heroic masculine figure, shaped by cultural memory, mythic exaggeration, and evolving political ideologies. Beowulf is positioned here as a man whose identity is not only historical or literary, but mythological, molded to meet the gender expectations of the societies that retold his story.

From a gendered lens, this portrayal reinforces a masculine heroic archetype rooted in supernatural achievement. By emphasizing “supernatural Character” over mortal vulnerability, Beowulf is gendered as more than man, he is mythologized masculinity, capable of performing feats that symbolize the ultimate virtues of patriarchal societies: strength, courage, conquest, and leadership.

Linguistically, the passage’s formal register (“futile Speculation,” “mixed Personage”) reflects the Victorian scholarly tone, but also subtly upholds patriarchal values by assuming the centrality of male heroism as the proper subject of epic literature. The gender invisibility of women in both the narrative and its analysis further underlines the male-dominated interpretive tradition in early philology and mythology.

Comparatively, Beowulf as a mythologized male hero aligns with Feridoun in Shahnameh, Rama in Ramayana, and Siavash in Persian mythology. All are men elevated to semi-divine status, reinforcing a cultural preference for men as saviors, kings, and spiritual ideals. This canonization of male figures obscures feminine agency and often reinterprets communal or spiritual archetypes through a gendered lens of masculine dominance.

Its crazy to think how many kids and people died on my favorite date.On this date is my birthday something that I enjoy celebrating to others is the worst day.![] (https://giphy.com/gifs/BismarckStateCollege-happy-birthday-bismarck-state-college-73CsKNxdfQPc9gVOJX)

Camron Newcomb CC BY-NC-SA 4.0

This passage highlights Siavash’s mythological role as a “vegetation god” figure, symbolizing death and rebirth, a motif deeply tied to cycles of nature and agricultural fertility. In terms of gender politics, Siavash embodies a masculine hero archetype, whose sacrificial martyrdom and regeneration reflect culturally constructed ideals of male heroism, where strength is paired with self sacrifice for communal renewal.

Linguistically, the text frames Siavash as a divine masculine figure “drawn to earth to carry out their duty,” emphasizing active male agency and responsibility. The repeated focus on “martyrdom and regeneration” underscores a cultural valorization of male suffering as necessary for social and cosmic balance. This resonates with the patriarchal worldview of ancient Iranian and surrounding societies, where heroic masculinity is defined by endurance, sacrifice, and regeneration.

Moreover, the narrative subtly contrasts Siavash’s enduring symbolic vitality with the loss of “sanctity and old nature” following Aryan immigration, which may reflect the cultural and linguistic layers imposed by successive translators and editors, each influencing the gendered portrayal according to their historical context. For example, earlier texts may emphasize Siavash’s divine qualities, while later versions humanize him, aligning heroism with mortal virtues.

Comparatively, this construction of Siavash parallels other vegetation gods like Tammuz and Osiris, where masculine death and rebirth cycles serve as metaphors for heroic masculinity, blending divine and human traits. Unlike many epic female figures who embody passivity or relational roles, Siavash’s heroism is active and sacrificial, a key marker of masculine ideals across cultures.

Camron Newcomb CC BY-NC-SA 4.0

In this pivotal moment, Kaweh asserts his heroic masculinity through fearless defiance against the Shah’s corrupt authority. His refusal to “lend [his] hand unto this lie” marks a bold moral stance, positioning him as a masculine ideal rooted in honor, integrity, and resistance to tyranny. The physical act of tearing the declaration symbolizes the rejection of false authority and the destructive power of oppressive patriarchy embodied by Zohak’s regime.

Linguistically, the choice of words, “wicked,” “ignoble,” and “ally of Deevs”, reflects a clear moral binary tied to gendered power, Kaweh embodies righteous masculinity, while the Shah and his associates are cast as corrupt and weak. The narrative empowers Kaweh’s masculine agency, highlighting his ability to act alone and command respect (“none dared uplift a finger to restrain him”), underscoring patriarchal values that prioritize male leadership and courage.

Comparatively, other versions of this story and adaptations may soften or amplify Kaweh’s defiance depending on the translator’s cultural context and gender politics. For instance, a more modern feminist influenced version might explore Kaweh’s role in a communal or collaborative context, but this traditional version centers on individual male heroism as the driver of social justice. This focus parallels other epic heroes like Gilgamesh or Beowulf, where masculinity is inseparable from heroic authority and moral righteousness.

If every object that scientist measure has uncertainty how do they evaluate if it is "good enough" for a important decision on a experiment or project?

This is a major contract from America's 'individualistic' ideology. For example, American advertising and TV commercial typically focus on you as an individual and sell things that will make you better, stronger, smarter, make you stand out from others, etc. Filipino advertising is focused on encouraging you to do better for your family or to show how each person makes an impact on the community (example: https://youtu.be/K9vFWA1rnWc?si=CZwV3a1wRgag5flx)

This was my inspiration for researching values, culture and beliefs in contrast to U.S. culture. The Philippines has a much greater appreciation and honor for their elders, seek guidance and never disrespect their elders. For example, I learned when I lived there that I never initiate a conversation with an elder and I must show my respect by bowing down, reaching out for their hand and bringing their hand to my forehead to show my gratitude and seek their approval before moving forward.

Both the U.S. and the Philippines have a focus on marketing women empowerment

Additionally, The U.S. influence is prevalent to this day because the U.S. freed the Philippines from Japanese rule and the American influence is seen in Filipino culture such as basketball being the top sport in the Philippines (most NBA fans/viewership per capita) and the Philippines celebrates the 4th of July as independence day.

This is an extreme opposite of the United States, who avoid faith-driven businesses and focus on mission-driven businesses that do not isolate any other religions. For example, about 18% of Americans identify as Catholic. In the Philippines, nearly 82% identify as Catholic and this guides their business standards and branding.

Filipino advertising promotes advertising based on collective inclusivity and growth, stressing faith and women empowerment to pave a path to communal inclusivity. U.S. advertising really focuses on reaching a particular demographic to inspire individuals from that minority demographic to step up or that you are noticed, often separately from a society or communal setting.

They didn't even care about him, they didn't even see him

The imagery shown makes me want to annotate so there's can see

Camron Newcomb

CC BY-NC-SA 4.0

This passage presents Gilgamesh not just as a mythical hero, but as part of a historically rooted dynastic tradition, with his lineage traced through a priestly father. This connection between priesthood and kingship reflects how masculine authority in ancient Mesopotamian heroism is both divine and hereditary. Heroism is gendered male from its very origin, the right to rule and to be remembered is passed from man to man, sanctified by both blood and religion.

Interestingly, while Langdon's translation was significant in making this version accessible, later scholars (like Clay and Jastrow) criticized it for errors and misreadings, many of which reinforce patriarchal norms through selective emphasis. Langdon frequently positions Gilgamesh’s actions in a romanticized light, elevating masculine conquest and omitting or downplaying the influence of female characters like Shamhat or Ninsun.

Linguistically, this version contributes to the gendered image of the hero through epithets like “builder of walls” or “conqueror,” aligning Gilgamesh with male coded acts of power and civilization. The omission of female influence in Langdon’s translation suggests not only a flaw in scholarship but a reflection of early 20th century gender norms, where masculinity was seen as synonymous with leadership, and femininity as peripheral or subversive.

this is just for site and apps that offer porn

He acts as if the words the fans are saying are checkmarks exclamation and asterisks.

this just showed very vivid imagery

I chose this as I relate to this as thinking of their face and thinking who they are as a person and who wrote that.

Note: This response was posted by the corresponding author to Review Commons. The content has not been altered except for formatting.

Learn more at Review Commons

Reply to the reviewers

We thank the reviewers for their positive comments. Our manuscript is to our knowledge the first to investigate the role of VAIL (V-ATPase—ATG16L1 induced LC3 lipidation), a form of CASM (Conjugation of ATG8s to single membranes) in SARS-CoV-2 replication. We demonstrate that SARS-CoV-2 Envelope (E) induces VAIL and this contributes to viral replication, including by using a reverse genetics system to make an E mutant virus. There have been many high quality studies examining the role of canonical autophagy in SARS-CoV-2 replication and our manuscript does not argue that all or even most LC3 lipidation during infection is via VAIL. We will try to make this point more clearly in the text. We do not think this detracts from the novelty and importance of our manuscript.

*Reviewer #1 (Evidence, reproducibility and clarity (Required)): *

• Figueras-Novoa et al present a short report demonstrating the induction of LC3 lipidation on single membranes by SARS-CoV-2 through a noncanonical autophagy pathway referred to as VAIL. The authors utilize elegant genetic tools to show that the induction of LC3 lipidation upon viral infection is mainly due to VAIL rather than canonical autophagy. They demonstrate that the activity of the viral E protein that can cause neutralization of acidic vesicles leads to the activation of non-canonical LC3 lipidation on single membranes. Interestingly, the authors also conclude that the impairment of VAIL leads to a reduction of viral load as a result of a defect in later stages of viral infection, although the underlying mechanism was not further explored. *

• Overall, this is an elegant and well controlled study that provides a clear conclusion. I only have some minor comments.*

We thank the reviewer for their assessment of our manuscript.

In some experiments, LC3 lipidation does not appear to be fully disrupted upon VAIL inhibition (e.g. Fig.'s 1H, 3D, 4A). As other labs have shown that SARS-CoV2 blocks autophagic flux, this could be further clarified in this manuscript as both VAIL and autophagy may be co-induced upon viral infection.

We agree with the reviewer that there is a contribution of canonical macroautophagy to the LC3B lipidation observed in SARS-CoV-2. We will extend the discussion in the manuscript to clarify this point for the readers.

Can the authors test the induction of LC3 lipidation in cells expressing K490 mutant of ATG16L1 in ATG16L1 KO cells to compare them with ATG16L1-ATG13 double knockouts?

The western blot in figure 3F (quantified in Figure 3G) shows LC3B lipidation in response to E expression in ATG16L1-ATG13 double knock out cells reconstituted with wild type ATG16L1 but not in cells complimented with ATG16L1 K490A mutant. We agree that the referee’s suggestion to perform these experiments in the context of infection would be informative. However in spite of numerous attempts, we have so far been unable to generate a cell clone fully devoid of ATG16L1 in a cell line that can be productively infected with SARS-CoV-2. For reasons unclear to us there appears to be a very low level of residual ATG16L1 activity despite multiple different CRISPR/Cas9 targeting attempts. The suggested complementation experiments might still be informative in the context of low level ATG16L1 expression so we will pursue this. Alternatively, as a contingency we can try to produce SARS-CoV-2 infectable cells with mutations in ATG16L1’s binding partner V1H, this interaction is required for VAIL. A further contingency could be to assess LC3B lipidation during infection and treatment with a Vps34 inhibitor, which inhibits canonical autophagy.

Minor points: * * The difference between Fig. 1F&G is unclear and why the authors are including both analyses. Similarly figures 4G&H.

We included both metrics to show that the decrease in LC3B lipidation in cells expressing SopF during infection is robust and observed in two separate readouts. While spot area measures the area of infected cells covered by GFP-LC3B fluorescence, spot intensity is a reading of the intensity of the area defined in an infected cell as being LC3 positive. Theoretically, these measurements could change in different ways. For example, if the same amount of lipidated LC3 were to distribute over a larger area of the cell. We prefer to keep both measurements in the manuscript.

The authors should show boxed colocalisation of all images, including negative controls. For examples, the authors have shown boxed magnifications in only the lowest panel in Figure 2A but not the upper two panels. Figures 4E&F should include boxed examples. This serves to clarify both positive and negative colocalisation events.

Boxed magnifications will be added to all images.

• Reviewer #1 (Significance (Required)): *

• Overall an elegant and well controlled study demonstrating the induction of non-canonical LC3 conjugation on single membranes (VAIL) during SARS-CoV2 infection. A further exploration of canonical autophagy (as previously published by others) in addition to VAIL would enhance this study.*

As the reviewer noted, several excellent studies have explored canonical autophagy during SARS-CoV-2 infection, many of which we cite in our manuscript. Our focus, however, is to demonstrate that SARS-CoV-2 E induces LC3 lipidation via VAIL. We believe that exploring the diverse roles of canonical autophagy mechanisms in SARS-CoV-2 infection is beyond the scope of this study.

*This study is of interest to researchers studying autophagy, viruses, immunology, single membrane LC3 lipidation, and lysosomes as well as potentially clinicians treating SARS-CoV2 infecteted individuals. *

• This reviewer is experienced in autophagy research.*

We thank the reviewer for this assessment of our manuscript.

*Reviewer #2 (Evidence, reproducibility and clarity (Required)): *

• Major Comments *

• Figure 1D does not very clearly show an overlap between V1D and LC3B. Both proteins seem broadly present across the cell and there is no easily identifiable change in V1D distribution upon infection. As such the overlay may be purely stochastic. The authors should quantify the observed co-localization events across multiple cells and biological replicates and compare them to other protein(s) with a similar cellular distribution pattern.*

We agree there is no obvious change in V1D staining on infection. The images in Figure 1D are purely intended to illustrate that LC3 and the V-ATPase can colocalise, not to demonstrate a change in V-ATPase distribution or to suggest a direct interaction. We will make this point more clearly in the text. We will also carry out analyses of the kind (see also response to the first two Minor Comments). We would be happy to provide an alternative method of visualising the V-ATPase (we could use any suitable antibody to the V-ATPase, or the bacterial effector SidK) if required. In response to reviewer 3’s comments, we will carry out a pull-down experiment to test the association of the V-ATPase and ATG16L1 during E expression, as this is a key interaction during VAIL activation.

Based on Figure 2F the authors suggest that virus entry is unaffected by the inhibition of VAIL in early timepoints. However, according to the figure legend, the timepoint used is 7hpi, while 2D uses 24hpi. Some SARS-CoV-2 papers suggest 7-10 hours is sufficient time to release new virions (Ban-On et al., 2020). As such 7hpi can not necessarily be seen as an early time point. Did the authors test earlier ones? Also, based on this, would it be possible that the effects observed at 24hpi are actually secondary infections, meaning that the virus utilizes pathway components for virion production and a lack thereof reduces infectivity of newly formed virions? In this case it would be interesting to set up an assay that can distinguish between primary and secondary infection to study both individually more closely.

Whereas 7 hours may be sufficient to release new virions, it is not sufficient to establish infections in other cells – this is why we chose that time point. The observation that there is no difference in the percentage of infected cells at 7 h p.i. (figure 2F) led us to suggest that viral entry is unaffected . We then confirmed this through the pseudovirus assay in Figure 2G, where no difference is found between SopF and mCherry expressing cells. For this assay, GFP-expressing, replication incompetent, lentiviral particles pseudotyped with Spike from different SARS-CoV-2 lineages were used to transduce mCherry and SopF expressing cells. A change in the percentage of GFP-positive cells would indicate an effect on viral entry, but no such change was observed in SopF-expressing cells.

We agree with the reviewer that the effects observed at 24 hpi are likely due to a defect in subsequent rounds of infection, since no difference was observed at 7 hpi or with our pseudovirus assay. We will attempt to make this point in the text as clearly as possible.

The authors nicely show in their study an involvement of VAIL in SARS-CoV-2 mediated LC3 lipidation. However, the observed effects are relatively moderate in several experiments, indicating that there may be another contributor to the observed phenotype. It would be nice to highlight this in the discussion and debate potential mechanisms that are causing the observed effects during infection.

We agree with the reviewer’s analysis. We have discussed the contribution of canonical autophagy in the second paragraph of the discussion, but we will expand on this in a revised manuscript. E expression levels are moderate during infection, other structural proteins such as N and M are present in much higher amounts. Since E is the key protein in VAIL initiation, a moderate effect of VAIL inhibition in perhaps expected. Nonetheless this still plays a crucial role in the viral life cycle.

*Minor Comments *

• The re-localization events shown in Fig 3A should be quantified.*

This quantification of GFP-LC3 relocalisation will be carried out and included.

• The co-localization events displayed in Fig 4A should be quantified.*

The quantification of V1D, E and GFP-LC3 will be carried out and included.

For Figure 2H-K the authors perform KDs of ATG16L1 and ATG13. While the results for the two specific proteins are certainly convincing, the authors would strengthen their argument by testing additional proteins in the autophagy pathway to support their claim that VAIL but not autophagy affects protein abundance of N (OPTIONAL).

As discussed in response to reviewer 1, we will attempt to infect ATG16L1 KO cells reconstituted with a K490A ATG16L1 mutant, which is an established tool and has been validated to be deficient in VAIL but not canonical autophagy.

***Referee cross-commenting** *

• Overall I agree with the comments of my co-reviewers and I think the suggested experiments/comments are sensible. *
• I in part already eluted to it my analysis, but I tend to agree with reviewer 3 on the limited effect VAIL seems to have on LC3b lipidation.*

As outlined above in response to reviewer 1 and below to reviewer 3, we agree that there is a modest contribution of VAIL to overall LC3 lipidation, which correlates with a modest amount of E expression in SARS-CoV-2 infection. VAIL is clearly important for the viral life cycle, thus whatever the proportion of LC3 lipidation attributable to this pathway it must be biologically significant.

*Reviewer #2 (Significance (Required)): *

• While previous publications have shown interaction between SARS-CoV2 and autophagy, the authors of this manuscript demonstrate that V-ATPase-ATG16L1 induced LC3 lipidation (VAIL) is activated during infection and affects viral replication. *

• This study provides an interesting new aspect to host-SARS_CoV-2 interactions. *

• The manuscript is of interest for people studying virus-host cell interaction, as well as for researchers in the fields of infectious diseases, specifically SARS-CoV2, and autophagy/VAIL*.

We thank the reviewer for their assessment of our manuscript.

R*eviewer #3 (Evidence, reproducibility and clarity (Required)): *

• The interaction of SARS-CoV-2 with canonical autophagy has been well documented. However, whether SARS-CoV-2 infection induces and benefits from non-canonical autophagy is unclear. In this manuscript, the authors demonstrated that SARS-CoV-2 infection induces V-ATPase-ATG16L1-induced LC3 lipidation (VAIL), a form of non-canonical autophagy in which LC3 is conjugated to single membranes. The SARS-CoV-2 envelope protein, through its ion channel activity, triggers the V-ATPase proton pump and induces VAIL during SARS-CoV-2 infection. Inhibiting VAIL during SARS-CoV-2 infection with SopF, a Salmonella effector, attenuates SARS-CoV-2 egress. *

• While these findings are interesting and demonstrate that SARS-CoV-2 infection triggers VAIL for its own benefit, the mechanism by which VAIL promotes SARS-CoV-2 replication remains unclear. Moreover, the contribution of VAIL to LC3 lipidation during SARS-CoV-2 infection appears to be minimal, as blocking VAIL through SoPF expression only marginally reduced LC3B lipidation (Fig. 1H). Therefore, the contribution of VAIL to LC3 lipidation during SARS-CoV-2 infection is minimal.*

We thank the reviewer for their assessment of our manuscript. As we have already alluded to in our response, we agree that only part of the LC3 lipidation observed during infection can be attributed to VAIL. There is a reproducible effect on viral replication which we have demonstrated in multiple ways, therefore the contribution of VAIL is of biological importance.

*Comments: *

• The authors show that the ion channel activity of E is essential for VAIL induction during SARS-CoV-2 infection. Since V-ATPase recruits the ATG16L complex to induce VAIL, and to clarify how SARS-CoV-2 infection triggers VAIL, the authors should examine whether SARS-CoV-2 infection or the expression of E induces V-ATPase-ATG16L interaction and whether this interaction is disrupted when SopF is expressed.*

We agree with the reviewer that this would be an informative experiment. We can carry out this experiment in an E expression system, rather than infection. This is due to the difficulty of getting enough material to carry out this kind of pull-down experiment in infected cells (at the time of writing these experiments still have to be carried out in CL3).

• Since the authors suggest that expression of SopF attenuates viral exit, one would expect that the number of N-positive cells will increase in SopF-expressing cells compared to the mCherry control cells. However, as shown in Figure 2D, this is not the case. Could the authors discuss why N-positive cells will be reduced in SopF-expressing cells when viral egress is impeded in these cells*?

This is a reflection of multi-cycle kinetics. N is still very strongly expressed in infected cells, even after virions have egressed. SARS-CoV-2 can infect VAIL-deficient cells and expresses the same levels of N prior to subsequent rounds of infection (at 7 hours after infection for example). Egress in VAIL-deficient, SopF-expressing cells is defective. Therefore, fewer cells will be infected in subsequent rounds of infection in SopF expressing cells, resulting in fewer N-positive cells in the SopF expressing cell population (most obvious after 24 hours).

Figure 2H. The authors show that knockdown of ATG16L1 reduces the expression of N during SARS-CoV-2 infection compared to the controls. To confirm that knockdown of ATG16L1, which is required for both canonical autophagy and VAIL, reduces N staining via VAIL, the authors should examine the impact of SopF expression on N levels in ATG16L KD cells. This experiment will confirm if the reduction in N staining in ATG16L1 KD cells is due to VAIL.

As stated in the response to reviewer 1, we can attempt this experiment in an ATG16L1 KO system complemented with K490A ATG16L1, which is deficient in VAIL and not canonical autophagy.

• Figure 2J. The quality of the Western blot data is poor.*

In this western the exposure is deliberately turned up to show that minimal ATG13 was left after knock down. We will also show the full blot with less exposure – this will demonstrate high quality.

Also, N appears as a single band in Figure 2J, but appears as double bands in Figures 2A and H. Could the authors explain this?

An extra band can be seen in 2J for N. However, as the reviewer points out, the intensity of the lower band is fainter than in 2A or 2H. The biology of SARS-CoV-2 N is interesting and complicated, with different truncated isoforms and phosphorylation patterns observed (see for example Mears et al., 2025 PMID:39836705). We observed changes in abundance of the second band between experiments, but this did not obviously depend on VAIL. We therefore consider this to be beyond the scope of this investigation.

*Reviewer #3 (Significance (Required)): *

• This manuscript proposes a role for VAIL in LC3 lipidation during SARS-CoV-2 infection. While the findings are interesting, VAIL only marginally contributes to LC3 lipidation during SARS-CoV-2 infection. Therefore, the significance of VAIL to LC3B lipidation during SARS-CoV-2 infection is unclear.*

Our experiments show unambiguously that VAIL contributes to viral replication. Therefore even if As alluded to above, we do not think a further investigation of canonical macroautophagy and SARS-CoV-2 would enhance the quality of our manuscript. We will try to make our description of the contribution of macroautophagy clearer in the revised manuscript (without providing a full literature review). We also do not think that exploring the nature of the multiple N bands on western blot is within the scope of this paper.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #3

Evidence, reproducibility and clarity

The interaction of SARS-CoV-2 with canonical autophagy has been well documented. However, whether SARS-CoV-2 infection induces and benefits from non-canonical autophagy is unclear. In this manuscript, the authors demonstrated that SARS-CoV-2 infection induces V-ATPase-ATG16L1-induced LC3 lipidation (VAIL), a form of non-canonical autophagy in which LC3 is conjugated to single membranes. The SARS-CoV-2 envelope protein, through its ion channel activity, triggers the V-ATPase proton pump and induces VAIL during SARS-CoV-2 infection. Inhibiting VAIL during SARS-CoV-2 infection with SopF, a Salmonella effector, attenuates SARS-CoV-2 egress.

While these findings are interesting and demonstrate that SARS-CoV-2 infection triggers VAIL for its own benefit, the mechanism by which VAIL promotes SARS-CoV-2 replication remains unclear. Moreover, the contribution of VAIL to LC3 lipidation during SARS-CoV-2 infection appears to be minimal, as blocking VAIL through SoPF expression only marginally reduced LC3B lipidation (Fig. 1H). Therefore, the contribution of VAIL to LC3 lipidation during SARS-CoV-2 infection is minimal.

Comments:

The authors show that the ion channel activity of E is essential for VAIL induction during SARS-CoV-2 infection. Since V-ATPase recruits the ATG16L complex to induce VAIL, and to clarify how SARS-CoV-2 infection triggers VAIL, the authors should examine whether SARS-CoV-2 infection or the expression of E induces V-ATPase-ATG16L interaction and whether this interaction is disrupted when SopF is expressed.

Since the authors suggest that expression of SopF attenuates viral exit, one would expect that the number of N-positive cells will increase in SopF-expressing cells compared to the mCherry control cells. However, as shown in Figure 2D, this is not the case. Could the authors discuss why N-positive cells will be reduced in SopF-expressing cells when viral egress is impeded in these cells?

Figure 2H. The authors show that knockdown of ATG16L1 reduces the expression of N during SARS-CoV-2 infection compared to the controls. To confirm that knockdown of ATG16L1, which is required for both canonical autophagy and VAIL, reduces N staining via VAIL, the authors should examine the impact of SopF expression on N levels in ATG16L KD cells. This experiment will confirm if the reduction in N staining in ATG16L1 KD cells is due to VAIL.

Figure 2J. The quality of the Western blot data is poor. Also, N appears as a single band in Figure 2J, but appears as double bands in Figures 2A and H. Could the authors explain this?

Significance

This manuscript proposes a role for VAIL in LC3 lipidation during SARS-CoV-2 infection. While the findings are interesting, VAIL only marginally contributes to LC3 lipidation during SARS-CoV-2 infection. Therefore, the significance of VAIL to LC3B lipidation during SARS-CoV-2 infection is unclear.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #2

Evidence, reproducibility and clarity

Major Comments

Figure 1D does not very clearly show an overlap between V1D and LC3B. Both proteins seem broadly present across the cell and there is no easily identifiable change in V1D distribution upon infection. As such the overlay may be purely stochastic. The authors should quantify the observed co-localization events across multiple cells and biological replicates and compare them to other protein(s) with a similar cellular distribution pattern.

Based on Figure 2F the authors suggest that virus entry is unaffected by the inhibition of VAIL in early timepoints. However, according to the figure legend, the timepoint used is 7hpi, while 2D uses 24hpi. Some SARS-CoV-2 papers suggest 7-10 hours is sufficient time to release new virions (Ban-On et al., 2020). As such 7hpi can not necessarily be seen as an early time point. Did the authors test earlier ones? Also, based on this, would it be possible that the effects observed at 24hpi are actually secondary infections, meaning that the virus utilizes pathway components for virion production and a lack thereof reduces infectivity of newly formed virions? In this case it would be interesting to set up an assay that can distinguish between primary and secondary infection to study both individually more closely.

The authors nicely show in their study an involvement of VAIL in SARS-CoV-2 mediated LC3 lipidation. However, the observed effects are relatively moderate in several experiments, indicating that there may be another contributor to the observed phenotype. It would be nice to highlight this in the discussion and debate potential mechanisms that are causing the observed effects during infection.

Minor Comments

The re-localization events shown in Fig 3A should be quantified.

The co-localization events displayed in Fig 4A should be quantified.

For Figure 2H-K the authors perform KDs of ATG16L1 and ATG13. While the results for the two specific proteins are certainly convincing, the authors would strengthen their argument by testing additional proteins in the autophagy pathway to support their claim that VAIL but not autophagy affects protein abundance of N (OPTIONAL).

Referee cross-commenting

Overall I agree with the comments of my co-reviewers and I think the suggested experiments/comments are sensible. I in part already eluted to it my analysis, but I tend to agree with reviewer 3 on the limited effect VAIL seems to have on LC3b lipidation.

Significance

While previous publications have shown interaction between SARS-CoV2 and autophagy, the authors of this manuscript demonstrate that V-ATPase-ATG16L1 induced LC3 lipidation (VAIL) is activated during infection and affects viral replication.

This study provides an interesting new aspect to host-SARS_CoV-2 interactions.

The manuscript is of interest for people studying virus-host cell interaction, as well as for researchers in the fields of infectious diseases, specifically SARS-CoV2, and autophagy/VAIL.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #1

Evidence, reproducibility and clarity

Figueras-Novoa et al present a short report demonstrating the induction of LC3 lipidation on single membranes by SARS-CoV-2 through a noncanonical autophagy pathway referred to as VAIL. The authors utilize elegant genetic tools to show that the induction of LC3 lipidation upon viral infection is mainly due to VAIL rather than canonical autophagy. They demonstrate that the activity of the viral E protein that can cause neutralization of acidic vesicles leads to the activation of non-canonical LC3 lipidation on single membranes. Interestingly, the authors also conclude that the impairment of VAIL leads to a reduction of viral load as a result of a defect in later stages of viral infection, although the underlying mechanism was not further explored.

Overall, this is an elegant and well controlled study that provides a clear conclusion. I only have some minor comments.

In some experiments, LC3 lipidation does not appear to be fully disrupted upon VAIL inhibition (e.g. Fig.'s 1H, 3D, 4A). As other labs have shown that SARS-CoV2 blocks autophagic flux, this could be further clarified in this manuscript as both VAIL and autophagy may be co-induced upon viral infection. Can the authors test the induction of LC3 lipidation in cells expressing K490 mutant of ATG16L1 in ATG16L1 KO cells to compare them with ATG16L1-ATG13 double knockouts?

Minor points:

The difference between Fig. 1F&G is unclear and why the authors are including both analyses. Similarly figures 4G&H.

The authors should show boxed colocalisation of all images, including negative controls. For examples, the authors have shown boxed magnifications in only the lowest panel in Figure 2A but not the upper two panels. Figures 4E&F should include boxed examples. This serves to clarify both positive and negative colocalisation events.

Significance

Overall an elegant and well controlled study demonstrating the induction of non-canonical LC3 conjugation on single membranes (VAIL) during SARS-CoV2 infection. A further exploration of canonical autophagy (as previously published by others) in addition to VAIL would enhance this study.

This study is of interest to researchers studying autophagy, viruses, immunology, single membrane LC3 lipidation, and lysosomes as well as potentially clinicians treating SARS-CoV2 infecteted individuals.

This reviewer is experienced in autophagy research.

Students engage in online discussion forums as well as interaction with peers in synchronous remote class sessions. Kendra frames the discussions with questions for personal reflection, guiding student preparation. This supports both multiple means of expression and engagement.

Course content includes multiple means of representation. Kendra assigns documentaries and legal documents alongside textbook chapters.

Note that course content includes multiple means of representation. Denise assigns chapters from the equity-minded textbook, dynamic web content, and relevant media relating to the week's topic.

Denise planned activities to allow for multiple means of engagement. This includes inviting multiple guest speakers, coordinating class trips to museums and the campus library, incorporating small group work into class sessions, and facilitating all-class discussion.

A personal exploration of Nigerian cuisine through vibrant food photography. Featuring dishes like beans with fried potatoes, jollof rice with plantain, and palm oil sauce with yam. While emphasizing the cultural significance and unique flavors of these meals.

Make this a separate item - is the data dashboard the project dashboard?

May publish the project application as a marker paper in iPlaces?

I think this is all part of the published project application/plan

the author is saying that if you don't want to read this u have every right to choice to do what you want, but that doesn't mean you should push your opinion on others because if you say its inappropriate that doesn't mean people can't have their own opinion on it.

This text is a personal narrative from a writer who describes their experience and censorship and book banning. This all serves as a direct ,clear statement of the authors position on the issue.

Hopkins being uninvited from the school event because a student's parent pointed out that her books may be inappropriate, For middle school students. Having her books removed from all middle schools in the district and her being prohibited from speaking at one is unfair.

The author is explaining the tension between the positive impact that their book has on people (readers ) like healing, connection, and some sort of hope

This was a unique and powerful example of how the writer's honest, personal story can directly save or change a person's life. It was inspired by their daughter's addiction and how it resonates with the young women in the audience.

Her book seemed to have inspired the woman into try and change her life, get away from drugs and turn her life around.

The author must've spoked to the lovely young woman and touched her heart based off talking about her journey

eLife Assessment

This study reveals that female moths use ultrasonic sounds emitted by dehydrated plants to guide their oviposition decisions. It highlights sound as an additional sensory modality in host searching, adding an important piece to the puzzle of how insects and plants interact. Through convincing experimental approaches, the authors provide insights that advance our understanding of plant-insect interactions.

Reviewer #1 (Public review):

Thank you for your thorough reply to my review and for conducting the additional experiments. All points from my previous review have been addressed.

Reviewer #2 (Public review):

This paper presents interesting and fresh approach as it investigates whether female moths utilize plant-emitted ultrasounds, particularly those associated with dehydration stress, in their egg-laying decision-making process. It provides the first empirical evidence suggesting that acoustic information may contribute to insect-plant interactions.

The revised version is significantly strengthened by the addition of supplementary data and improved explanations. The authors present robust results across multiple experiments, enhancing the credibility of their conclusions.

Female moths showed a preference for moist, fresh plants over dehydrated ones in experiments using actual plants. Additionally, when both plants were fresh but ultrasonic sounds specific to dehydrated plants were presented from one side, the moths chose the silent plant. However, in experiments without plants, contrary to the hypothesis derived from the above results, the moths preferred to oviposit near ultrasonic playback mimicking the sounds of dehydrated plants.　

These results clearly indicate that moths can perceive plant presence through sound.　The findings also highlight the need for future investigation into the multi-modal nature of moth decision-making, as acoustic cues alone may not fully explain the behavioral choices observed across different contexts.

Overall, the results are intriguing, and I think the experiments are very well designed. The authors successfully demonstrate that plant-derived acoustic signals influence oviposition behavior in female moths, thereby achieving the study's objectives. The experimental design and analysis protocols are reproducible and well suited for adaptation to other species.

Author response:

The following is the authors’ response to the original reviews.

Reviewer #1 (Public review):

(1) The authors demonstrate that female Spodoptera littoralis moths prefer to oviposit on wellwatered tomato plants and avoid drought-stressed plants. The study then recorded the sounds produced by drought-stressed plants and found that they produce 30 ultrasonic clicks per minute. Thereafter, the authors tested the response of female S. littoralis moths to clicks with a frequency of 60 clicks per minute in an arena with and without plants and in an arena setting with two healthy plants of which one was associated with 60 clicks per minute. These experiments revealed that in the absence of a plant, the moths preferred to lay eggs on the side of the area in which the clicks could be heard, while in the presence of a plant the S. littoralis females preferred to oviposit on the plant where the clicks were not audible. In addition, the authors also tested the response of S. littoralis females in which the tympanic membrane had been pierced making the moths unable to detect the click sounds. As hypothesised, these females placed their eggs equally on both sites of the area.

Finally, the authors explored whether the female oviposition choice might be influenced by the courtship calls of S. littoralis males which emit clicks in a range similar to a drought-stressed tomato plant. However, no effect was found of the clicks from ten males on the oviposition behaviour of the female moths, indicating that the females can distinguish between the two types of clicks. Besides these different experiments, the authors also investigated the distribution of egg clusters within a longer arena without a plant, but with a sugar-water feeder. Here it was found that the egg clusters were mostly aggregated around the feeder and the speaker producing 60 clicks per minute. Lastly, video tracking was used to observe the behaviour of the area without a plant, which demonstrated

that the moths gradually spent more time at the arena side with the click sounds.

We thank the reviewers for their helpful comments. We agree with the summary, but would like to note that in the control experiment (Figure 2) we used a click rate of 30 clicks per minute—a design choice driven by the editor’s feedback. We have clarified this and, to further probe the system’s dynamics, added a second experiment employing the same click rate (30 clicks per minute) with a dehydrated plant (see details below). In both experiments, females again showed a clear tendency to oviposit nearer the speaker; these findings are described in the updated manuscript.

(2) The study addresses a very interesting question by asking whether female moths incorporate plant acoustic signals into their oviposition choice, unfortunately, I find it very difficult to judge how big the influence of the sound on the female choice really is as the manuscript does not provide any graphs showing the real numbers of eggs laid on the different plants, but instead only provides graphs with the Bayesian model fittings for each of the experiments. In addition, the numbers given in the text seem to be relatively similar with large variations e.g. Figure 1B3: 1.8 {plus minus} 1.6 vs. 1.1 {plus minus} 1.0. Furthermore, the authors do not provide access to any of the raw data or scripts of this study, which also makes it difficult to assess the potential impact of this study. Hence, I would very much like to encourage the authors to provide figures showing the measured values as boxplots including the individual data points, especially in Figure 1, and to provide access to all the raw data underlying the figures.

We acknowledge that there are researchers who favor Bayesian graphical representation versus raw data visualization. Therefore, we have added chartplots of the raw data from Figure 1 in the supplementary section. We are aware of the duplication in presentation and apologize for this redundancy.  

Regarding the variance and means we obtained in our experiment, we have analyzed all raw data using the statistical model presented, and if statistical significance was found despite a particular mean difference or variance, this is meaningful from a biological perspective. One can certainly discuss whether this difference has biological importance, but it should be remembered that in this experimental system, we are trying to isolate the acoustic signal from a complex system that includes multiple signals. Therefore, at no point we’ve suggested that this is a standalone factor, but rather proposed it as an informative and significant component. 

In addition to the experiments described above, we conducted an experiment in which we counted both eggs and clusters. The results indicate that cluster counts are a reliable proxy for reproductive investment at a given location. In this experiment, we present cluster numbers alongside egg counts (Figure 2).

Furthermore, we apologize for the technical error that prevented our uploaded data files from reaching the reviewers. We have also uploaded updated data and code.

(3) Regarding the analysis of the results, I am also not entirely convinced that each night can be taken as an independent egg-laying event, as the amount of eggs and the place were the eggs are laid by a female moth surely depends on the previous oviposition events. While I must admit that I am not a statistician, I would suggest, from a biological point of view, that each group of moths should be treated as a replicate and not each night. I would therefore also suggest to rather analyse the sum of eggs laid over the different consecutive nights than taking the eggs laid in each night as an independent data point.

We thank the reviewer for this question. This is a valid and point that we will address in three aspects: 

First, regarding our statistical approach, we used a model that takes into account the sequence of nights and examines whether there is an effect of the order of nights, i.e., we used GLMMs, with the night nested within the repetition. This is equivalent to addressing this as a repeated measure and is, to our best knowledge, the common way to treat such data. 

Second, following the reviewer's comment, we also reran the statistics of the third experiment (i.e., “sound gradient experiments”, Figure 2 and Supplementary figure 4) when only taking the first night when the female/s laid eggs to avoid the concern of dependency. This analysis revealed the same result – i.e., a significant preference for the sound stimulus. We have now updated our methods and results section to clarify this point.  

Third, an important detail that may not have been clearly specified in the methods: at the end of each night, we cleaned the arena of counted egg clusters using a cloth with ethanol, so that on the subsequent night, we would not expect there to be evidence of previous oviposition but thus would not exclude some sort of physiological or cognitive memories. We have now updated our methods section to clarify this important procedural point. 

(4) Furthermore, it did not become entirely clear to me why a click frequency of 60 clicks per minute was used for most experiments, while the plants only produce clicks at a range of 30 clicks per minute. Independent of the ecological relevance of these sound signals, it would be nice if the authors could provide a reason for using this frequency range. Besides this, I was also wondering about the argument that groups of plants might still produce clicks in the range of 60 clicks per minute and that the authors' tests might therefore still be reasonable. I would agree with this, but only in the case that a group of plants with these sounds would be tested. Offering the choice between two single plants while providing the sound from a group of plants is in my view not the most ecologically reasonable choice. It would be great if the authors could modify the argument in the discussion section accordingly and further explore the relevance of different frequencies and dBlevels.

This is an excellent point. We originally increased the click rate generate a strong signal. However, it was important for us to verify that there was ecological relevance in the stimulus we implemented in the system. For this purpose, we recorded a group of dehydrated plants at a distance of ~20cm and we measured a click rate of 20 clicks per minute (i.e., 0.33 Hz) (see Methods section). Therefore, as mentioned at the beginning of this letter, in the additional experiment described in Figure 2, we reduced the click frequency to 30 clicks per minute, and at this lower rate, the effect was maintained. Increasing plant density would probably lead to a higher rate of 30 clicks per minute. 

(5) Finally, I was wondering how transferable the findings are towards insects and Lepidopterans in general. Not all insects possess a tympanic organ and might therefore not be able to detect the plant clicks that were recorded. Moreover, I would imagine that generalist herbivorous like Spodoptera might be more inclined to use these clicks than specialists, which very much rely on certain chemical cues to find their host plants. It would be great if the authors would point more to the fact that your study only investigated a single moth species and that the results might therefore only hold true for S. littoralis and closely related species, but not necessary for other moth species such as Sphingidae or even butterflies.

Good point. Our research uses a specific model system of one moth species and one plant species in a particular plant-insect interaction where females select host plants for their offspring. As with any model-based research that attempts to draw broader conclusions, we've taken care to distinguish between our direct findings and potential wider implications. We believe our system may represent mechanisms relevant to a wider group of herbivorous insects with hearing capabilities, particularly considering that several moth families and other insect orders can detect ultrasound. However, additional research examining more moth and plant species is necessary to determine how broadly applicable these findings are. We have made these clarifications in the text.

Reviewer #2 (Public review):

(6) The results are intriguing, and I think the experiments are very well designed. However, if female moths use the sounds emitted by dehydrated plants as cues to decide where to oviposit, the hypothesis would predict that they would avoid such sounds. The discussion mentions the possibility of a multi-modal moth decision-making process to explain these contradictory results, and I also believe this is a strong possibility. However, since this remains speculative, careful consideration is needed regarding how to interpret the findings based solely on the direct results presented in the results section.  

Thank you for this insightful observation. We agree that the apparent attraction of females to dehydrated-plant sounds contradicts our initial prediction. Having observed this pattern consistently across multiple setups, we have now added a targeted choice experiment to the revised manuscript: here female moths were offered a choice between dehydrated plants broadcasting their natural ultrasonic emissions and a control. These results—detailed in the Discussion and presented in full in the Supplementary Materials (Supplementary Figure 4)—show that when only a dehydrated plant is available, moths would prefer it for oviposition, supporting our hypothesis that in the absence of a real plant, the plant’s sounds might represent a plant..

(7) Additionally, the final results describing differences in olfactory responses to drying and hydrated plants are included, but the corresponding figures are placed in the supplementary materials. Given this, I would suggest reconsidering how to best present the hypotheses and clarify the overarching message of the results. This might involve reordering the results or re-evaluating which data should appear in the main text versus the supplementary materials

Thank you for this suggestion. We have reorganized the manuscript and removed the olfactory response data from the current version to maintain a focused narrative on acoustic cues. We agree that a detailed investigation of multimodal interactions deserves a separate study, which we plan to pursue in future work. 

(8) There were also areas where more detailed explanations of the experimental methods would be beneficial.

Thank you for highlighting this point. We have expanded and clarified the Methods section to provide comprehensive detail on our experimental procedures.

Reviewer #1 (Recommendations for the authors):

(9) Line 1: Please include the name of the species you tested also in the title as your results might not hold true for all moth species.

We do not fully agree with this comment. Please see comment 5.

(10) Line 19-20: Please rephrase the sentence so that it becomes clear that the "dehydration stress" refers to the plant and not to the moths.

Thank you for the suggestion; we have clarified the text accordingly

(11) Line 31: Male moths might provide many different signals to the females, maybe better "male sound signals" or similar.

Thank you for the suggestion; we have clarified the text accordingly.

(12) Line 52-53: Maybe mention here that not all moth species have evolved these abilities.

Thank you for the suggestion; we have clarified the text accordingly.

(13) Line 77: add a space after 38.

Thank you for the suggestion; we have clarified the text accordingly.

(14) Line 88: Maybe change "secondary predators" to "natural enemies".

Thank you for the suggestion; we have clarified the text accordingly.

(15) Line 134: Why is "notably" in italics? I would suggest using normal spelling/formatting rules here.

Thank you for the suggestion; we have clarified the text accordingly.

(16) Line 140-144: If you did perform the experiment also with the more ecological relevant playback rate, why not present these findings as your main results and use the data with the higher playback frequency as additional support?

Thank you for this suggestion. We agree that the ecologically relevant playback data are important; as described in detail at the beginning of this letter and also in comment 4, however, to preserve a clear and cohesive narrative, we have maintained the original ordering of this section. Nevertheless, the various experiments conducted in Figure 1 differ in several components from Figure 2 and the work that examined sounds in plant groups in the appendices. Therefore, we find it more appropriate to use them as supporting evidence for the main findings rather than creating a comparison between different experimental systems. For this reason, we chose to keep them as a separate description in "The ecological playback findings (Lines 140–144) remain fully described in the Results and serve to reinforce the main observations without interrupting the manuscript's flow.

(17) Line 146: Please explain already here how you deafened the moths.

Thank you for the suggestion; we have clarified the text accordingly.

(18) Line 181: should it be "male moths' " ?

Thank you for the suggestion; we have clarified the text accordingly.

(19) Line 215: Why is "without a plant" in italics? I would suggest using normal spelling/formatting rules here.

Thank you for the suggestion; we have clarified the text accordingly.

(20) Line 234: I do not understand why this type of statistic was used to analyse the electroantennogram (EAG) results. Would a rather simple Student's t-test or a Wilcon rank sum test not have been sufficient? I would also like to caution you not to overinterpret the data derived from the EAG, as you combined the entire headspace into one mixture it is no longer possible to derive information on the different volatiles in the blends. The differences you observe might therefore mostly be due to the amount of emitted volatiles.

We have reorganized the manuscript and removed the olfactory response data from the current version to maintain a focused narrative on acoustic cues (See comment 7). 

(21) Line 268: It might be nice to add an additional reference here referring to the multimodal oviposition behaviour of the moths.

Thank you for the suggestion; we have clarified the text accordingly.

(22) Line 284: If possible, please add another reference here referring to the different cues used by moths during oviposition.

Thank you for the suggestion; we have clarified the text accordingly.

(23) Line 336: What do you mean by "closed together"?

Thank you for the suggestion; we have clarified the text accordingly.

(24) Line 434-436: Please see my overall comments. I do not think that you can call it ecologically relevant if the signal emitted by multiple plants is played in the context of just a single plant.

Please see comments 1 and 4.

(25) Line 496: Please change "stats" to statistics.

Thank you for the suggestion; we have clarified the text accordingly.

(26) Line 522-524: I am not sure whether simply listing their names does give full credit to the work these people did for your study. Maybe also explain how they contributed to your work.

Thank you for the suggestion; we have clarified the text accordingly.

Reviewer #2 (Recommendations for the authors):

(27) L54 20-60kHz --> 20Hz-60kHz or 20kHz - 60kHz?

OK. We have replaced it.

(28) L124 Are the results for the condition where nothing was placed and the condition where a decoy silent resistor was placed combined in the analysis? If so, were there no significant differences between the two conditions? Comparing these with a condition presenting band-limited noise in the same frequency range as the drought-stressed sounds might also have been an effective approach to further isolate the specific role of the ultrasonic emissions.

We have used both conditions due to technical constrains and pooled them tougher for analysis— statistical tests confirmed no significant differences between them—and this clarification has now been added to the Methods section including the results of the statistical test.

(29) L125 (Fig. 1A), see Exp. 1 in the Methods). -> (Fig.1B. See Exp.1 in the Methods).

Thank you for the suggestion; we have clarified the text accordingly.

(30) L132 "The opposite choice to what was seen in the initial experiment (Fig.1B)"

Thank you for the suggestion; we have clarified the text accordingly.

(31) L137-143 If you are writing about results, why not describe them with figures and statistics? The current description reads like a discussion.

These findings were not among our primary research questions; however, we believe that including them in the Results section underscores the experimental differences. In our opinion, introducing an additional figure or expanding the statistical analysis at this point would disrupt the narrative flow and risk confusing the reader.

(32) L141 "This is higher than the rate reported for a single young plant" Are you referring to the tomato plants used in the experiments? It might be helpful to include in the main text the natural click rate emitted by tomato plants, as this information is currently only mentioned in the Methods section.

See comment 4.  

(33) L191 Is the main point here to convey that the plant playback effect remained significant even when the sound presentation frequency was reduced to 30 clicks per minute? The inclusion of the feeder element, however, seems to complicate the message. To simplify the results, moving the content from lines 185-202 to the supplementary materials might be a better approach. Additionally, what is the rationale for placing the sugar solution in the arena? Is it to maintain the moths' vitality during the experiment? Clarifying this in the methods section would help provide context for this experimental detail.

In this series of experiments, we manipulated four variables—single moths, ultrasonic click rate, arena configuration (from a two-choice design to an elongated enclosure), and the response metric (total egg counts rather than cluster counts)—to evaluate moth oviposition under more ecologically realistic conditions. We demonstrate the system’s robustness and validity in a more realistic setting (by tracking individual moths, counting single eggs, etc.).  

As noted in the text, feeders were included to preserve the moths’ natural behavior and vitality. We have further clarified this in the revised manuscript.

(34) L215 Is the click presentation frequency 30 or 60 per minute? Since Figure 3 illustrates examples of moth movement from the experiment described in Figure 1, it might be more effective to present Figure 3 when discussing the results of Figure 1 or to include it in the supplementary materials for better clarity and organization.

See comments 1 and 4. As mentioned in the above 

(35) L291 Please provide a detailed explanation of the experiments and measurements for the results shown in Figure S3 (and Figure S2). If the multi-modal hypothesis discussed in the study is a key focus, it might be better to include these results in the main results section rather than in the supplementary materials.

Thank you for this suggestion. Figure S2 was removed, see comments above. We’ve added now the context to figure S3.

(36) L303 It might be helpful to include information about the relationship between the moth species used in this study and tomato plants somewhere in the text. This would provide an important context for understanding the ecological relevance of the experiments.

Thank you for the suggestion; we have clarified the text accordingly.

(37) Table 1 The significant figures in the numbers presented in the tables should be consistent.

Thank you for the suggestion; we have clarified the text accordingly.

(38) L341 The text mentions that experiments were conducted in a greenhouse, but does this mean the arena was placed inside the greenhouse? Also, the term "arena" is used - does this refer to a sealed rectangular case or something similar? For the sound presentation experiments, it seems that the arena cage was placed inside a soundproof room. If the arena is indeed a case-like structure, were there any specific measures taken to prevent sound scattering within the case, such as the choice of materials or structural modifications?

Here, “arena” refers to the plastic boxes used throughout this study. In this particular experiment, we presented plants alone—reflecting ongoing debate in the literature—and used these trials as a baseline for our subsequent sound-presentation experiments, during which we measured sound intensity as described in the Methods section. All sound-playback experiments were conducted in sound-proof rooms, and acoustic levels were measured beforehand—sound on the control side fell below our system’s detection threshold. 

(39) L373 "resister similar to the speaker" Could you explain it in more detail? I think this would depend on the type of speaker used-particularly whether it includes magnets. From an experimental perspective, presenting different sounds such as white noise from the speaker might have been a better control. Was there a specific reason for not doing so? Additionally, the study does not clearly demonstrate whether the electric and magnetic field environments on both sides of the arena were appropriately controlled. Without this information, it is difficult to evaluate whether using a resistor as a substitute was adequate.

Thank you for this comment. We have now addressed this point in the Discussion. We acknowledge that we did not account for the magnetic field, which might have differed between the speaker and the resistor. We agree that using an alternative control, such as white noise, could have been informative, and we now mention this as a limitation in the revised Methods.

(40) L435 60Hz? The representation of frequencies in the text is inconsistent, with some values expressed in Hz and others as "clicks per second." It would be better to standardize these units for clarity, such as using Hz throughout the manuscript.

We agree that this is confusing. We reviewed the text and made sure that when we addressed click per second, we meant how many clicks were produced and when we addressed Hz units it was in the context of sound frequencies.  

(41) L484 "we quantified how many times each individual crossed the center of the arena" Is this data being used in the results?

Yes. Mentioned in the text just before Figure 3. L220

eLife Assessment

IL-10 balances protective and deleterious immune functions in mice and humans, but if IL-10 also controls avian intestinal homeostasis remains unclear. Generating genetic knockouts, Meunier et al. established that a complete lack of IL-10 strengthened immunity against enteric bacteria in chickens, while also aggravating infection-inflicted inflammatory tissue damage and dysbiosis upon parasite infection, but unlike mouse models, IL-10 deficiency did not provoke spontaneous colitis in chickens. The findings presented are valuable, and the strength of evidence is convincing. The observation may have implications for the livestock industry and additional studies involving genetic knockouts may further unravel conserved and distinct avian IL-10 controls.

Reviewer #1 (Public review):

Summary:

In this study, Meunier et al. investigated the functional role of IL-10 in avian mucosal immunity. While the anti-inflammatory role of IL-10 is well established in mammals, and several confirmatory Knock-out models available in mice, IL-10's role in avian mucosal immunity is so far correlative. In this study the authors generated two different models of IL-10 ablation in Chickens. A whole body knock-out model, and an enhancer KO model leading to reduced IL10 expression. The authors first performed in vitro LPS stimulation based experiments, and then in vivo two different infection models employing C. jejuni and E. tenella, to demonstrate that complete ablation of IL10 leads to enhanced inflammation related pathology and gene expression, and enhanced pathogen clearance. At a steady-state level, however, IL-10 ablation did not lead to spontaneous colitis.

Strengths:

Overall the study is well executed and establishes an anti-inflammatory role of IL-10 in birds. While the results are expected, and not surprising, this appears to be the first report to conclusively demonstrate IL-10's anti-inflammatory role upon its genetic ablation in avian model. Provided the applicability of this information in combating pathogen infection in livestock species in sustainable industries like poultry, the study is worth publishing.

Weaknesses:

The study is primarily a confirmation of the already established anti-inflammatory role of IL-10.

Comments on revised version:

The authors have incorporated most of the points raised, and provided a reasonable argument for not considering DSS mediated colitis as an additional model.

Reviewer #2 (Public review):

Summary:

The authors were to investigate functional role of IL10 on mucosal immunity in chickens. CRISPR technology was employed to generate IL10 knock out chickens in both exon and putative enhancer regions. IL10 expressions were either abolished (knockout in exon) or reduced (enhancer knock-out). IL-10 play an important role in the composition of the caecal microbiome. Through various enteric pathogens challenge, deficient IL10 expression was associated with enhanced pathogen clearance, but with more severe lesion score and body weight loss.

Strengths:

Both in vitro and in vivo knock-out in abolished and reduced IL10 expression and broad enteric pathogens were challenged in vivo and various parameters were examined to evaluate the functional role of IL10 on mucosal immunity.

Weaknesses:

Overexpression of IL10 either in vitro or in vivo may further support the findings from this study.

Comments on revised version:

The authors' response and justifications are appropriate.

Author response:

The following is the authors’ response to the original reviews.

Reviewer #1 (Public review): 

Summary: 

In this study, Meunier et al. investigated the functional role of IL-10 in avian mucosal immunity. While the anti-inflammatory role of IL-10 is well established in mammals, and several confirmatory knockout models are available in mice, IL-10's role in avian mucosal immunity is so far correlative. In this study, the authors generated two different models of IL-10 ablation in Chickens. A whole body knock-out model and an enhancer KO model leading to reduced IL10 expression. The authors first performed in vitro LPS stimulation-based experiments, and then in vivo two different infection models employing C. jejuni and E. tenella, to demonstrate that complete ablation of IL10 leads to enhanced inflammation-related pathology and gene expression, and enhanced pathogen clearance. At a steady-state level, however, IL-10 ablation did not lead to spontaneous colitis. 

Strengths: 

Overall, the study is well executed and establishes an anti-inflammatory role of IL-10 in birds. While the results are expected and not surprising, this appears to be the first report to conclusively demonstrate IL-10's anti-inflammatory role upon its genetic ablation in the avian model. Provided this information is applicable in combating pathogen infection in livestock species in sustainable industries like poultry, the study will be of interest to the field. 

Weaknesses: 

The study is primarily a confirmation of the already established anti-inflammatory role of IL-10. 

We do not agree that this work is primarily confirmatory. The anti-inflammatory role of IL10 was indeed known previously from studies in mammals. The much more general insight from the current study is our demonstration of the intrinsic trade-off between inflammation and tolerance in the response to both the microbiome (which was significantly altered in the IL10 knockout birds) and mucosal pathogens. The study of Eimeria challenge in particular highlights the fact that it may be better for the host to tolerate a potential pathogen than to take on the cost of elimination.

Reviewer #2 (Public review): 

Summary: 

The authors were to investigate the functional role of IL10 on mucosal immunity in chickens. CRISPR technology was employed to generate IL10 knock-out chickens in both exon and putative enhancer regions. IL10 expressions were either abolished (knockout in exon) or reduced (enhancer knock-out). IL-10 plays an important role in the composition of the caecal microbiome. Through various enteric pathogen challenges, deficient IL10 expression was associated with enhanced pathogen clearance, but with more severe lesion scores and body weight loss. 

Strengths: 

Both in vitro and in vivo knock-out abolished and reduced IL10 expression, and broad enteric pathogens were challenged in vivo, and various parameters were examined to evaluate the functional role of IL10 on mucosal immunity. 

Weaknesses: 

Overexpression of IL-10 either in vitro or in vivo may further support the findings from this study. 

An overexpression experiment, regardless of outcome, would not necessarily support or invalidate the findings of the current study. It would address the question of whether the absolute concentration of IL10 produced alters the outcome of an infection.

Reviewer #1 (Recommendations for the authors): 

The following are the recommendations that, in my opinion, will be helpful to enhance the quality of the study. 

Major point: 

The authors at a steady state did not observe any sign of spontaneous colitis. Since IL-10 KO in mice leads to enhanced pathological score upon DSS-mediated induction of colitis, and several colitis models are well established in birds, it will be worthwhile to test the consequence of experimentally inducing colitis in this context. 

One of the novel features of this study is the observation that the microbiome is modified in the IL10KO HOM chicks, which may serve to mitigate potential spontaneous pathology; we now mention this in the discussion. We agree that it could be worthwhile in the future to look at additional challenge models. However, we would argue that the Eimeria challenge is a sufficiently adequate experimentally-induced model of colitis to demonstrate the increased inflammation that occurs in an IL10-deficient bird. This is further supported by evidence of enhanced inflammatory responses in the caeca of IL10KO HOM birds challenged with Campylobacter or Salmonella relative to WT controls. See in the revised manuscript (pages 12-13).

Minor points: 

(1) In Figure 2B, the authors should confirm whether the ROS-AV163 groups also have LPS treatment. 

The legend for Figure 2B already states that neutralizing anti-IL10 antibody was added to LPS-stimulated BMDMs: “Nitric oxide production was assessed by measuring nitrite levels using Griess assay for LPS-stimulated BMDMs […] in the absence or presence of neutralizing anti-IL10 antibody ROS-AV163”. However, for added clarity we have now modified the x-axis label for Figure 2B (“+ROS-AV163” replaced by “+LPS +anti-IL10”) and we have also made minor changes to the figure legend. See in the revised manuscript (page 33).

(2) In Figure 3F, the authors should discuss why the duodenum of KO birds has enhanced infiltration compared to WT? 

We are not sure what the reviewer is referring to here. Although not specifically mentioned in Figure 3F, there is no statistically significant difference in cellular infiltration in the duodenum of IL10KO WT and HOM birds raised in our specified pathogen-free (SPF) facility, nor in the duodenum of IL10KO WT and HOM birds raised in our conventional facility (Mann-Whitney U tests, p>0.1 in both cases); this can be seen in the sums of histopathological scores shown in Figures 3C (SPF facility) and 3E (conventional facility). Figure 3F shows that there is a statistically significant difference in cellular infiltration scores in the duodenum and proximal colon of both IL10KO WT and HOM birds based on the environment they are raised in (SPF vs conventional). We have made minor changes to the text to clarify this. See in the revised manuscript (page 7).

(3) The authors should discuss the observed differences in the C. jejuni colonization results among the two cohorts at week 1 and week 2 post-infection. 

Numbers of C. jejuni in the caeca of IL10KO HOM birds were markedly lower than for WT controls at 1-week post-infection in cohort 1, and at both time intervals post-infection in cohort 2 (Figure 4A). This reached statistical significance at 1-week post-infection in cohort 1 and at 2-weeks post-infection in cohort 2. It is evident from Figure 4A that considerable inter-animal variance existed in each group, and in the IL10KO HOM birds in particular. This is typical of C. jejuni colonisation in chickens, where bacterial population structures have been reported to be variable and unpredictable (Coward et al., Appl Environ Microbiol 2008, PMID: 18424530). Similar variation between time intervals, birds and repeated experiments has been reported when evaluating vaccines against C. jejuni colonisation (e.g. Buckley et al., Vaccine 2010, PMID: 19853682; Nothaft et al., Front Microbiol 2021, PMID: 34867850). We performed two independent studies for this reason. Taken together, we consider that our data provide convincing evidence of elevated pro-inflammatory responses upon C. jejuni infection in IL10KO HOM birds relative to WT controls that associates with reduced bacterial burden. Our data is also consistent with a published observation that a commercial broiler line with low IL10 expression had correspondingly elevated expression of CXCLi-1, CXCLi-2 and IL-1b (Humphrey et al., mBio 2014, reference 33 in our original submission). We have added text to the discussion to capture the points above.  See in the revised manuscript (page 13).

Reviewer #2 (Recommendations for the authors): 

For the animal challenging experiments, both IL10KO HOM and IL10EnKO HOM chickens were used for Eimeria challenge, but not for Salmonella and Campylobacter. Could the authors justify why? 

The Eimeria challenge produced a much higher and more reproducible level of inflammation than either of the bacterial challenge models. Within the parasite challenge cohorts, IL10KO HET and IL10EnKO HOM birds were only marginally different from WT controls (e.g. parasite replication: Figures 5A and B; lesion scores: Figures 5E and F; body weight gain: Figures 5G and H). Given the more limited response and the inter-individual variation in the bacterial challenge models, we felt that analysis of a sufficiently large cohort of the IL10KO HOM was appropriate, while additional cohorts of IL10KO HET and IL10EnKO HOM birds large enough to detect statistically significant differences could not be justified.

In the M&M, there was no mention of # of birds generated for IL10EnKO HOM, HET, etc. 

Full details of bird numbers can be found in SI Appendix Table S1 “Number of IL10KO and IL10EnKO WT, HET and HOM chicks hatched in the NARF SPF chicken facility in the first (G1) and second (G2) generations”. Table S1 is already referred to in the Results section “Generation of IL10-deficient chickens”; we have now also clearly referred to it in the “Animals” and “Generation of surrogate host chickens and establishment of the IL10KO and IL10EnKO lines under SPF conditions” sections of the Materials and Methods. In all three sections we have also added some text to clarify that the table details G1 and G2 bird numbers. See in the revised manuscript (pages 5, 15, 17).

From the results of Campylobacter challenge, the results from the cohort 1 and cohort 2 were not consistent at both 1 and 2 weeks of post-infection. There is not much discussion on this inconsistency. What is the final conclusion: significant difference in week 1 or week 2, OR none of them, OR both of them. What would happen if an additional cohort were conducted for Salmonella and Eimeria? 

As noted in response to Reviewer 1 (minor point 3), we have now added text to the discussion on the partial inconsistency between independent C. jejuni challenge studies. We do not feel that additional experiments to address this comment are required. Highly significant increases in the infiltration of lymphoplasmacytic cells and heterophils were detected in IL10KO HOM chickens relative to WT controls in the caeca, a key site of Campylobacter colonisation. This was consistently observed in two independent cohorts at both 1- and 2-weeks post-infection (SI Appendix Figures S7 and S8) and was reflected in similar patterns of expression of pro-inflammatory genes at these intervals in both cohorts (Figure 4B). As our laboratory has observed substantially less variation between repeated Salmonella challenges, a single study was performed, but with adequate power to detect statistical differences.  The effects of E. tenella infection in IL10KO WT and HOM birds were replicated (compare Figure 4 with data from day 6 in Figure 5).

eLife Assessment

The authors present a software (TEKRABber) to analyze how expression of transposable elements (TEs) and TE silencing factors KRAB zinc finger (KRAB-ZNF) genes are correlated in experimentally validated datasets. TEKRABber is used to reconstruct regulatory networks of KRAB-ZNFs and TEs during human brain evolution and in Alzheimer's disease. The direction of the work is important, with potentially significant interest from others looking for a tool for correlative gene expression analysis across individual genomes and species. However, the reviews identified biases and shortcomings in the pipeline that could lead to an unacceptable number of false positive and negative signals and thus impact the conclusions, leaving the work in its current form incomplete.

Reviewer #1 (Public review):

The authors present their new bioinformatic tool called TEKRABber, and use it to correlate expression between KRAB ZNFs and TEs across different brain tissues, and across species. While the aims of the authors are clear and there would be significant interest from other researchers in the field for a program that can do such correlative gene expression analysis across individual genomes and species, the presented approach and work display significant shortcomings. In the current state of the analysis pipeline, the biases and shortcomings mentioned below, for which I have seen no proof of that they are accounted for by the authors, are severely impacting the presented results and conclusions. It is therefore essential that the points below are addressed, involving significant changes in the TEKRABber progamm as well as the analysis pipeline, to prevent the identification of false positive and negative signals, that would severely affect the conclusions one can raise about the analysis.

My main concerns are provided below:

One important shortcoming of the biocomputational approach is that most TEs are not actually expressed, and others (Alus) are not a proxy of the activity of the TE class at all. I will explain: While specific TE classes can act as (species-specific) promoters for genes (such as LTRs) or are expressed as TE derived transcripts (LINEs, SVAs), the majority of other older TE classes do not have such behavior and are either neutral to the genome or may have some enhancer activity (as mapped in the program they refer to 'TEffectR'. A big focus is on Alus, but Alus contribute to a transcriptome in a different way too: They often become part of transcripts due to alternative splicing. As such, the presence of Alu derived transcripts is not a proxy for the expression/activity of the Alu class, but rather a result of some Alus being part of gene transcripts (see also next point). Bottom line is that the TEKRABber software/approach is heavily prone to picking up both false positives (TEs being part of transcribed loci) and false negatives (TEs not producing any transcripts at all) , which has a big implication for how reads from TEs as done in this study should be interpreted: The TE expression used to correlate the KRAB ZNF expression is simply not representing the species-specific influences of TEs where the authors are after.

With the strategy as described, a lot of TE expression is misinterpreted: TEs can be part of gene-derived transcripts due to alternative splicing (often happens for Alus) or as a result of the TE being present in an inefficiently spliced out intron (happens a lot) which leads to TE-derived reads as a result of that TE being part of that intron, rather than that TE being actively expressed. As a result, the data as analysed is not reliably indicating the expression of TEs (as the authors intend too) and should be filtered for any reads that are coming from the above scenarios: These reads have nothing to do with KRAB ZNF control, and are not representing actively expressed TEs and therefore should be removed. Given that from my lab's experience in brain (and other) tissues, the proportion of RNA sequencing reads that are actually derived from active TEs is a stark minority compared to reads derived from TEs that happen to be in any of the many transcribed loci, applying this filtering is expected to have a huge impact on the results and conclusions of this study.

Another potential problem that I don't see addressed is that due to the high level of similarity of the many hundreds of KRAB ZNF genes in primates and the reads derived from them, and the inaccurate annotations of many KZNFs in non-human genomes, the expression data derived from RNA-seq datasets cannot be simply used to plot KZNF expression values, without significant work and manual curation to safeguard proper cross species ortholog-annotation: The work of Thomas and Schneider (2011) has studied this in great detail but genome-assemblies of non-human primates tend to be highly inaccurate in appointing the right ortholog of human ZNF genes. The problem becomes even bigger when RNA-sequencing reads are analyzed: RNA-sequencing reads from a human ZNF that emerged in great apes by duplication from an older parental gene (we have a decent number of those in the human genome) may be mapped to that older parental gene in Macaque genome: So, the expression of human-specific ZNF-B, that derived from the parental ZNF-A, is likely to be compared in their DESeq to the expression of ZNF-A in Macaque RNA-seq data. In other words, without a significant amount of manual curation, the DE-seq analysis is prone to lead to false comparisons which make the stategy and KRABber software approach described highly biased and unreliable.

There is no doubt that there are differences in expression and activity of KRAB-ZNFs and TEs repspectively that may have had important evolutionary consequences. However, because all of the network analyses in this paper rely on the analyses of RNA-seq data and the processing through the TE-KRABber software with the shortcomings and potential biases that I mentioned above, I need to emphasize that the results and conclusions are likely to be significantly different if the appropriate measures are taken to get more accurate and curated TE and KRAB ZNF expression data.

Finally, there are some minor but important notes I want to share:

The association with certain variations in ZNF genes with neurological disorders such as AD, as reported in the introduction is not entirely convincing without further functional support. Such associations could be merely happen by chance, given the high number of ZNF genes in the human genome and the high chance that variations in these loci happen associate with certatin disease associated traits. So using these associations as an argument that changes in TEs and KRAB ZNF networks are important for diseases like AD should be used with much more caution.

There is a number of papers where KRAB ZNF and TE expression are analysed in parallel in human brain tissues. So the novelty of that aspect of the presented study may be limited.

Additional note after reviewing the revised version of the manuscript:

After reviewing the revised version of the manuscript, my criticism and concerns with this study are still evenly high and unchanged. To clarify, the revised version did not differ in essence from the original version; it seems that unfortunately, no efforts were taken to address the concerns raised on the original version of the manuscript, the results section as well as the discussion section are virtually unchanged.

Author response:

The following is the authors’ response to the current reviews.

Reviewer #1 (Public review): 

The authors present their new bioinformatic tool called TEKRABber, and use it to correlate expression between KRAB ZNFs and TEs across different brain tissues, and across species. While the aims of the authors are clear and there would be significant interest from other researchers in the field for a program that can do such correlative gene expression analysis across individual genomes and species, the presented approach and work display significant shortcomings. In the current state of the analysis pipeline, the biases and shortcomings mentioned below, for which I have seen no proof of that they are accounted for by the authors, are severely impacting the presented results and conclusions. It is therefore essential that the points below are addressed, involving significant changes in the TEKRABber progamm as well as the analysis pipeline, to prevent the identification of false positive and negative signals, that would severely affect the conclusions one can raise about the analysis. 

Thank you very much for the insightful review of our manuscript. Since most of the comments on our revised version are not different from the comments on our first version, we repeated our previous answer, but wrote a new reply to the new concerns (please see the last two paragraphs). 

We would also like to reiterate here that most of the critique of the reviewer concerns the performance of other tools and not TEKRABber presented in our manuscript. We consider it out of scope for this manuscript to improve other tools.

My main concerns are provided below: 

One important shortcoming of the biocomputational approach is that most TEs are not actually expressed, and others (Alus) are not a proxy of the activity of the TE class at all. I will explain: While specific TE classes can act as (species-specific) promoters for genes (such as LTRs) or are expressed as TE derived transcripts (LINEs, SVAs), the majority of other older TE classes do not have such behavior and are either neutral to the genome or may have some enhancer activity (as mapped in the program they refer to 'TEffectR'. A big focus is on Alus, but Alus contribute to a transcriptome in a different way too: They often become part of transcripts due to alternative splicing. As such, the presence of Alu derived transcripts is not a proxy for the expression/activity of the Alu class, but rather a result of some Alus being part of gene transcripts (see also next point). Bottom line is that the TEKRABber software/approach is heavily prone to picking up both false positives (TEs being part of transcribed loci) and false negatives (TEs not producing any transcripts at all) , which has a big implication for how reads from TEs as done in this study should be interpreted: The TE expression used to correlate the KRAB ZNF expression is simply not representing the species-specific influences of TEs where the authors are after. 

With the strategy as described, a lot of TE expression is misinterpreted: TEs can be part of gene-derived transcripts due to alternative splicing (often happens for Alus) or as a result of the TE being present in an inefficiently spliced out intron (happens a lot) which leads to TE-derived reads as a result of that TE being part of that intron, rather than that TE being actively expressed. As a result, the data as analysed is not reliably indicating the expression of TEs (as the authors intend too) and should be filtered for any reads that are coming from the above scenarios: These reads have nothing to do with KRAB ZNF control, and are not representing actively expressed TEs and therefore should be removed. Given that from my lab's experience in brain (and other) tissues, the proportion of RNA sequencing reads that are actually derived from active TEs is a stark minority compared to reads derived from TEs that happen to be in any of the many transcribed loci, applying this filtering is expected to have a huge impact on the results and conclusions of this study. 

We sincerely thank the reviewer for highlighting the potential issues of false positives and negatives in TE quantification. The reviewer provided valuable examples of how different TE classes, such as Alus, LTRs, LINEs, and SVAs, exhibit distinct behaviors in the genome. To our knowledge, specific tools like ERVmap (Tokuyama et al., 2018), which annotates ERVs, and LtrDetector (Joseph et al., 2019), which uses k-mer distributions to quantify LTRs, could indeed enhance precision by treating specific TE classes individually. We acknowledge that such approaches may yield more accurate results and appreciate the suggestion. 

In our study, we used TEtranscripts (Jin et al., 2015) prior to TEKRABber. TEtranscripts applies the Expectation Maximization (EM) algorithm to assign ambiguous reads as the following steps. Uniquely mapped reads are first assigned to genes, and  reads overlapping genes and TEs are assigned to TEs only if they do not uniquely match an annotated gene. The remaining ambiguous reads are distributed based on EM iterations. While this approach may not be as specialized as the latest tools for specific TE classes, it provides a general overview of TE activity. TEtranscripts outputs subfamily-level TE expression data, which we used as input for TEKRABber to perform downstream analyses such as differential expression and correlation studies.

We understand the importance of adapting tools to specific research objectives, including focusing on particular TE classes. TEKRABber is designed not to refine TE quantification at the mapping stage but to flexibly handle outputs from various TE quantification tools. It accepts raw TE counts as input in the form of dataframes, enabling diverse analytical pipelines. We would also like to clarify that, since the input data is transcriptomic, our primary focus is on expressed TEs, rather than the effects of non-expressed TEs in the genome. In the revised version of our manuscript, we emphasize this distinction in the discussion and provide examples of how TEKRABber can integrate with other tools to enhance specificity and accuracy.

Another potential problem that I don't see addressed is that due to the high level of similarity of the many hundreds of KRAB ZNF genes in primates and the reads derived from them, and the inaccurate annotations of many KZNFs in non-human genomes, the expression data derived from RNA-seq datasets cannot be simply used to plot KZNF expression values, without significant work and manual curation to safeguard proper cross species ortholog-annotation: The work of Thomas and Schneider (2011) has studied this in great detail but genome-assemblies of non-human primates tend to be highly inaccurate in appointing the right ortholog of human ZNF genes. The problem becomes even bigger when RNA-sequencing reads are analyzed: RNA-sequencing reads from a human ZNF that emerged in great apes by duplication from an older parental gene (we have a decent number of those in the human genome) may be mapped to that older parental gene in Macaque genome: So, the expression of human-specific ZNF-B, that derived from the parental ZNF-A, is likely to be compared in their DESeq to the expression of ZNF-A in Macaque RNA-seq data. In other words, without a significant amount of manual curation, the DE-seq analysis is prone to lead to false comparisons which make the stategy and KRABber software approach described highly biased and unreliable. 

There is no doubt that there are differences in expression and activity of KRAB-ZNFs and TEs repspectively that may have had important evolutionary consequences. However, because all of the network analyses in this paper rely on the analyses of RNA-seq data and the processing through the TE-KRABber software with the shortcomings and potential biases that I mentioned above, I need to emphasize that the results and conclusions are likely to be significantly different if the appropriate measures are taken to get more accurate and curated TE and KRAB ZNF expression data. 

We thank the reviewer for raising the important issue of accurately annotating the expanded repertoire of KRAB-ZNFs in primates, particularly the challenges of cross-species orthology and potential biases in RNA-seq data analysis. Indeed, we have also addressed this challenge in some of our previous papers (Nowick et al., 2010, Nowick et al., 2011 and Jovanovic et al., 2021).

In the revised manuscript, we include more details about our two-step strategy to ensure accurate KRAB-ZNF ortholog assignments. First, we employed the Gene Order Conservation (GOC) score from Ensembl BioMart as a primary filter, selecting only one-to-one orthologs with a GOC score above 75% across primates. This threshold, recommended in Ensembl’s ortholog quality control guidelines, ensures high-confidence orthology relationships.(http://www.ensembl.org/info/genome/compara/Ortholog_qc_manual.html#goc).

Second, we incorporated data from Jovanovic et al. (2021), which independently validated KRAB-ZNF orthologs across 27 primate genomes. This additional layer of validation allowed us to refine our dataset, resulting in the identification of 337 orthologous KRAB-ZNFs for differential expression analysis (Figure S2).

We acknowledge that different annotation methods or criteria may for some genes yield variations in the identified orthologs. However, we believe that this combination provides a robust starting point for addressing the challenges raised, while we remain open to additional refinements in future analyses.

Finally, there are some minor but important notes I want to share:

The association with certain variations in ZNF genes with neurological disorders such as AD, as reported in the introduction is not entirely convincing without further functional support. Such associations could be merely happen by chance, given the high number of ZNF genes in the human genome and the high chance that variations in these loci happen associate with certatin disease associated traits. So using these associations as an argument that changes in TEs and KRAB ZNF networks are important for diseases like AD should be used with much more caution. 

We fully acknowledge the concern that, given the large number of KRAB-ZNFs and their inherent variability, some associations with AD or other neurological disorders could occur by chance. This highlights the importance of additional functional studies to validate the causal role of KRAB-ZNF and TE interactions in disease contexts. While previous studies have indeed analyzed KRAB-ZNF and TE expression in human brain tissues, our study seeks to expand on this foundation by incorporating interspecies comparisons across primates. This approach enabled us to identify TE:KRAB-ZNF pairs that are uniquely present in healthy human brains, which may provide insights into their potential evolutionary significance and relevance to diseases like AD.

In addition to analyzing RNA-seq data (GSE127898 and syn5550404), we have cross-validated our findings using ChIP-exo data for 159 KRAB-ZNF proteins and their TE binding regions in humans (Imbeault et al., 2017). This allowed us to identify specific binding events between KRAB-ZNF and TE pairs, providing further support for the observed associations. We agree with the reviewer that additional experimental validations, such as functional studies, are critical to further establish the role of KRAB-ZNF and TE networks in AD. We hope that future research can build upon our findings to explore these associations in greater detail.

There is a number of papers where KRAB ZNF and TE expression are analysed in parallel in human brain tissues. So the novelty of that aspect of the presented study may be limited. 

We agree with the reviewer that many studies have examined the expression levels of KRAB-ZNFs and TEs in developing human brain tissues (Farmiloe et al., 2020; Turelli et al., 2020; Playfoot et al., 2021, among others). However, the novelty of our study lies in comparing KRAB-ZNF and TE expression across primate species, as well as in adult human brain tissues from both control individuals and those with Alzheimer’s disease. To our knowledge, no previous study has analyzed these data in this context. We therefore believe that our results will be of interest to evolutionary biologists and neurobiologists focusing on Alzheimer’s disease.

Additional note after reviewing the revised version of the manuscript: 

After reviewing the revised version of the manuscript, my criticism and concerns with this study are still evenly high and unchanged. To clarify, the revised version did not differ in essence from the original version; it seems that unfortunately, no efforts were taken to address the concerns raised on the original version of the manuscript, the results section as well as the discussion section are virtually unchanged.

We regret that this reviewer was not satisfied with our changes. In fact, many of the points raised by this reviewer are important, but concern weaknesses of other tools. In our opinion, validating other tools would be out of scope for this paper. We want to emphasize that TEKRABber is not a quantification tool for sequencing data, but a software for comparative analysis across species. We provided a detailed answer to the reviewer and readers can refer to that answer in the public review above for further information.

The following is the authors’ response to the original reviews.

Reviewer #1 (Public review):

The authors present their new bioinformatic tool called TEKRABber, and use it to correlate expression between KRAB ZNFs and TEs across different brain tissues, and across species. While the aims of the authors are clear and there would be significant interest from other researchers in the field for a program that can do such correlative gene expression analysis across individual genomes and species, the presented approach and work display significant shortcomings. In the current state of the analysis pipeline, the biases and shortcomings mentioned below, for which I have seen no proof that they are accounted for by the authors, are severely impacting the presented results and conclusions. It is therefore essential that the points below are addressed, involving significant changes in the TEKRABber program as well as the analysis pipeline, to prevent the identification of false positive and negative signals, that would severely affect the conclusions one can raise about the analysis.

Thank you very much for the insightful review of our manuscript.

My main concerns are provided below:

(1) One important shortcoming of the biocomputational approach is that most TEs are not actually expressed, and others (Alus) are not a proxy of the activity of the TE class at all. I will explain: While specific TE classes can act as (species-specific) promoters for genes (such as LTRs) or are expressed as TE derived transcripts (LINEs, SVAs), the majority of other older TE classes do not have such behavior and are either neutral to the genome or may have some enhancer activity (as mapped in the program they refer to 'TEffectR'. A big focus is on Alus, but Alus contribute to a transcriptome in a different way too: They often become part of transcripts due to alternative splicing. As such, the presence of Alu derived transcripts is not a proxy for the expression/activity of the Alu class, but rather a result of some Alus being part of gene transcripts (see also next point). The bottom line is that the TEKRABber software/approach is heavily prone to picking up both false positives (TEs being part of transcribed loci) and false negatives (TEs not producing any transcripts at all), which has a big implication for how reads from TEs as done in this study should be interpreted: The TE expression used to correlate the KRAB ZNF expression is simply not representing the species-specific influences of TEs where the authors are after.

With the strategy as described, a lot of TE expression is misinterpreted: TEs can be part of gene-derived transcripts due to alternative splicing (often happens for Alus) or as a result of the TE being present in an inefficiently spliced out intron (happens a lot) which leads to TE-derived reads as a result of that TE being part of that intron, rather than that TE being actively expressed. As a result, the data as analysed is not reliably indicating the expression of TEs (as the authors intend to) and should be filtered for any reads that are coming from the above scenarios: These reads have nothing to do with KRAB ZNF control, and are not representing actively expressed TEs and therefore should be removed. Given that from my lab's experience in the brain (and other) tissues, the proportion of RNA sequencing reads that are actually derived from active TEs is a stark minority compared to reads derived from TEs that happen to be in any of the many transcribed loci, applying this filtering is expected to have a huge impact on the results and conclusions of this study.

We sincerely thank the reviewer for highlighting the potential issues of false positives and negatives in TE quantification. The reviewer provided valuable examples of how different TE classes, such as Alus, LTRs, LINEs, and SVAs, exhibit distinct behaviors in the genome. To our knowledge, specific tools like ERVmap (Tokuyama et al., 2018), which annotates ERVs, and LtrDetector (Joseph et al., 2019), which uses k-mer distributions to quantify LTRs, could indeed enhance precision by treating specific TE classes individually. We acknowledge that such approaches may yield more accurate results and appreciate the suggestion. 

In our study, we used TEtranscripts (Jin et al., 2015) prior to TEKRABber. TEtranscripts applies the Expectation Maximization (EM) algorithm to assign ambiguous reads as the following steps. Uniquely mapped reads are first assigned to genes, and  reads overlapping genes and TEs are assigned to TEs only if they do not uniquely match an annotated gene. The remaining ambiguous reads are distributed based on EM iterations. While this approach may not be as specialized as the latest tools for specific TE classes, it provides a general overview of TE activity. TEtranscripts outputs subfamily-level TE expression data, which we used as input for TEKRABber to perform downstream analyses such as differential expression and correlation studies.

We understand the importance of adapting tools to specific research objectives, including focusing on particular TE classes. TEKRABber is designed not to refine TE quantification at the mapping stage but to flexibly handle outputs from various TE quantification tools. It accepts raw TE counts as input in the form of dataframes, enabling diverse analytical pipelines. We would also like to clarify that, since the input data is transcriptiomic, our primary focus is on expressed TEs, rather than the effects of non-expressed TEs in the genome. In the revised version of our manuscript, we emphasize this distinction in the discussion and provide examples of how TEKRABber can integrate with other tools to enhance specificity and accuracy.

(2) Another potential problem that I don't see addressed is that due to the high level of similarity of the many hundreds of KRAB ZNF genes in primates and the reads derived from them, and the inaccurate annotations of many KZNFs in non-human genomes, the expression data derived from RNA-seq datasets cannot be simply used to plot KZNF expression values, without significant work and manual curation to safeguard proper cross species ortholog-annotation: The work of Thomas and Schneider (2011) has studied this in great detail but genome-assemblies of non-human primates tend to be highly inaccurate in appointing the right ortholog of human ZNF genes. The problem becomes even bigger when RNA-sequencing reads are analyzed: RNA-sequencing reads from a human ZNF that emerged in great apes by duplication from an older parental gene (we have a decent number of those in the human genome) may be mapped to that older parental gene in Macaque genome: So, the expression of human-specific ZNF-B, that derived from the parental ZNF-A, is likely to be compared in their DESeq to the expression of ZNF-A in Macaque RNA-seq data. In other words, without a significant amount of manual curation, the DE-seq analysis is prone to lead to false comparisons which make the strategy and KRABber software approach described highly biased and unreliable.

There is no doubt that there are differences in expression and activity of KRAB-ZNFs and TEs respectively that may have had important evolutionary consequences. However, because all of the network analyses in this paper rely on the analyses of RNA-seq data and the processing through the TE-KRABber software with the shortcomings and potential biases that I mentioned above, I need to emphasize that the results and conclusions are likely to be significantly different if the appropriate measures are taken to get more accurate and curated TE and KRAB ZNF expression data.

We thank the reviewer for raising the important issue of accurately annotating the expanded repertoire of KRAB-ZNFs in primates, particularly the challenges of cross-species orthology and potential biases in RNA-seq data analysis. Indeed, we have also addressed this challenge in some of our previous papers (Nowick et al., 2010, Nowick et al., 2011 and Jovanovic et al., 2021).

In the revised manuscript, we include more details about our two-step strategy to ensure accurate KRAB-ZNF ortholog assignments. First, we employed the Gene Order Conservation (GOC) score from Ensembl BioMart as a primary filter, selecting only one-to-one orthologs with a GOC score above 75% across primates. This threshold, recommended in Ensembl’s ortholog quality control guidelines, ensures high-confidence orthology relationships. (http://www.ensembl.org/info/genome/compara/Ortholog_qc_manual.html#goc).

Second, we incorporated data from Jovanovic et al. (2021), which independently validated KRAB-ZNF orthologs across 27 primate genomes. This additional layer of validation allowed us to refine our dataset, resulting in the identification of 337 orthologous KRAB-ZNFs for differential expression analysis (Figure S2).

We acknowledge that different annotation methods or criteria may for some genes yield variations in the identified orthologs. However, we believe that this combination provides a robust starting point for addressing the challenges raised, while we remain open to additional refinements in future analyses.

(3) The association with certain variations in ZNF genes with neurological disorders such as AD, as reported in the introduction is not entirely convincing without further functional support. Such associations could merely happen by chance, given the high number of ZNF genes in the human genome and the high chance that variations in these loci happen to associate with certain disease-associated traits. So using these associations as an argument that changes in TEs and KRAB ZNF networks are important for diseases like AD should be used with much more caution.

There are a number of papers where KRAB ZNF and TE expression are analysed in parallel in human brain tissues. So the novelty of that aspect of the presented study may be limited.

We fully acknowledge the concern that, given the large number of KRAB-ZNFs and their inherent variability, some associations with AD or other neurological disorders could occur by chance. This highlights the importance of additional functional studies to validate the causal role of KRAB-ZNF and TE interactions in disease contexts. While previous studies have indeed analyzed KRAB-ZNF and TE expression in human brain tissues, our study seeks to expand on this foundation by incorporating interspecies comparisons across primates. This approach enabled us to identify TE:KRAB-ZNF pairs that are uniquely present in healthy human brains, which may provide insights into their potential evolutionary significance and relevance to diseases like AD.

In addition to analyzing RNA-seq data (GSE127898 and syn5550404), we have cross-validated our findings using ChIP-exo data for 159 KRAB-ZNF proteins and their TE binding regions in humans (Imbeault et al., 2017). This allowed us to identify specific binding events between KRAB-ZNF and TE pairs, providing further support for the observed associations. We agree with the reviewer that additional experimental validations, such as functional studies, are critical to further establish the role of KRAB-ZNF and TE networks in AD. We hope that future research can build upon our findings to explore these associations in greater detail.

Reviewer #1 (Recommendations for the authors):

It is essential before this work can be considered for publication, that the points above are addressed, involving significant changes in the TEKRABber program as well as the analysis pipeline, to prevent the identification of false positive and negative signals, that would severely affect the conclusions one can raise about the analysis.

We sincerely appreciate the reviewer’s insightful recommendations and constructive feedback. Each specific point has been carefully addressed in detail in the public reviews section above.

Reviewer #2 (Public review)

Summary:

The aim was to decipher the regulatory networks of KRAB-ZNFs and TEs that have changed during human brain evolution and in Alzheimer's disease.

Strengths:

This solid study presents a valuable analysis and successfully confirms previous assumptions, but also goes beyond the current state of the art.

Weaknesses:

The design of the analysis needs to be slightly modified and a more in-depth analysis of the positive correlation cases would be beneficial. Some of the conclusions need to be reinterpreted.

We sincerely thank the reviewer for the thoughtful summary, positive evaluation of our study, and constructive feedback. We appreciate the recognition of the strengths in our analysis and the valuable suggestions for improving its design and interpretation. 

We would like to briefly comment on the suggested modifications to the design here and will provide a detailed point-by-point review later with our revised manuscript. 

The reviewer recommended considering a more recent timepoint, such as less than 25 million years ago (mya), to define the "evolutionary young group" of KRAB-ZNF genes and TEs when discussing the arms-race theory. This is indeed a valuable perspective, as the TE repressing functions by KRAB-ZNF proteins  may have evolved more recently than the split between Old World Monkeys (OWM) and New World Monkeys (NWM) at 44.2 mya we used. 

Our rationale for selecting 44.2 mya is based on certain primate-specific TEs such as the Alu subfamilies, which emerged after the rise of Simiiformes and have been used in phylogenetic studies (Xing et al., 2007 and Williams et al., 2010). This timeframe allowed us to investigate the potential co-evolution of KRAB-ZNFs and TEs in species that emerged after the OWM-NWM split (e.g., humans, chimpanzees, bonobos, and macaques used for this study). However, focusing only on KRAB-ZNFs and TEs younger than 25 million years would limit the analysis to just 9 KRAB-ZNFs and 92 TEs expressed in our datasets. While we will not conduct a reanalysis using this more recent timepoint, we will integrate the recommendation into the discussion section of the revised manuscript. 

Furthermore, we greatly appreciate the reviewer's detailed insights and suggestions for refining specific descriptions and interpretations in our manuscript. We will address these points in the revised version to ensure the content is presented with greater precision and clarity.

Once again, we thank both reviewers for their valuable feedback, which provides significant input for strengthening our study.

Reviewer #2 (Recommendations for the authors):

We thank the reviewer for the very insightful comments, which helped a lot in our interpretation and discussion of our results and in improving some of our statements.

The present study seeks to uncover how the repression of transposable elements (TEs) by rapidly evolving KRAB-ZNF genes, which are known for their role in TE suppression, may influence human brain evolution and contribute to Alzheimer's disease (AD). Utilizing their previously developed tool, TEKRABber, the researchers analyze transcriptome datasets from the brains of four species of Old World Monkeys (OWM) alongside samples from healthy human individuals and AD patients.

Through bipartite network analysis, they identify KRAB-ZNF/Alu-TE interactions as the most negatively correlated in the network, highlighting the repression of Alu elements by KRAB-ZNF proteins. In AD patient samples, they observe a reduction in a subnetwork comprising 21 interactions within an Alu TE module. These findings are consistent with earlier evidence that: (1) KRAB-ZNFs are involved in suppressing evolutionarily young Alu TEs; and (2) specific Alu elements have been reported to be deregulated in AD. The study also validates previous experimental ChIP-exo data on KRAB-ZNF proteins obtained in a different cell type (Imbeault et al., 2017).

As a novely, the study identifies a human-specific amino acid variation in ZNF528, which directly contacts DNA nucleotides, showing signs of positive selection in humans and several human-specific TE interactions.

Interestingly, in addition to the negative links, the researchers observed predominantly positive connections with other TEs, suggesting that while their approach is consistent with some previous observations, the authors conclude that it provides limited support for the 'genetic arms race' hypothesis.

The reviewer is a specialist in TE and evolutionary research.

Major issues:

The study demonstrates the usefulness of the TEKRABber tool, which can support and successfully validate previous observations. However, there are several misconceptions and problems with the interpretation of the results.

KRAB-ZNF proteins in repressing TEs in vertebrates  In the Abstract: "In vertebrates, some KRAB-ZNF proteins repress TEs, offering genomic protection."

Although some KRAB-ZNF proteins exist in vertebrates, their TE-suppression role is not as prominent or specialized as it is in mammals, where it serves as a key defense mechanism against the mobilization of TEs.

We appreciate the reviewer’s clarification regarding the role of KRAB-ZNF proteins in vertebrates. To improve accuracy and precision, we have revised the wording to specify that this mechanism is primarily observed in mammals rather than vertebrates.

The definition of young and old

The study considers the evolutionary age of young ({less than or equal to} 44.2 mya) and old(> 44.2 mya). This is the time of the Old World Monkey (OWM) and New World Monkey (NWM) split. Importantly, however, the KRAB-ZNF / KAP1 suppression system primarily suppresses evolutionarily younger TEs (< 25 MY old). These TEs are relatively new additions to the genome, i.e. they are specific to certain lineages (such as primates or hominins) and are more likely to be actively transcribed (and recognized as foreign by innate immunity) or have residual activity upon transposition. Examples include certain subfamilies of LINE-1, Alu (Y, S, less effective for J), SVA and younger human endogenous retroviruses (HERVs) such as HERV-K. The KRAB-ZNF / KAP1 system therefore focuses primarily on TEs that have evolved more recently in primates, in the last few million years (within the last 25 million years). Older TEs are controlled by broader epigenetic mechanisms such as DNA methylation, histone modifications, etc. Therefore, the age ({less than or equal to} 44.2 mya) is not suitable to define it as young.

In this context, the specific TEs of the Simiiformes cannot be considered as 'recently evolved' (in the Abstract). The Simiiformes contain both OWM and NWM. Notably, the study includes four species, all of which belong to the OWMs.

The 'genetic arms race' theory

Unfortunately, the problematic definition of young and old could also explain why the authors conclude that their data only weakly support the 'genetic arms race' hypothesis.

The KRAB-ZNF proteins evolve rapidly, similar to TEs, which raises the 'genetic arms race' hypothesis. This hypothesis refers to the constant evolutionary struggle between organisms and TEs. TEs constantly evolve to overcome host defences, while host genomes develop mechanisms to suppress these potentially harmful elements. Indeed, in mammals, an important example is the KRAB-ZNF/TE interaction. The KRAB-ZNF proteins rapidly evolve to target specific TEs, creating a 'genetic arms race' in which each side - TEs and the KRAB-ZNF/KAP1 (alias TRIM28) repressor complex - drives the evolution of the other in response to adaptive pressure. Importantly, the 'genetic arms race' hypothesis describes the evolutionary process that occurs between TE and host when the TE is deleterious. Again, this includes the young TEs (< 25 MY old) with residual transposition activity or those that actively transcribed and exacerbate cellular stress and inflammatory responses. Approximately 25 million years ago, the superfamilies Hominoidea (apes) and Cercopithecoidea (Old World monkeys, I.e. macaque) split.

Just to clarify, our initial study aim was to examine whether TEs exhibit any evolutionary relationships with KRAB-ZNFs across the four studied species (human, chimpanzee, bonobo, and macaque). For investigating the arms-race hypothesis, we really appreciate the reviewer suggesting a more recent time point, such as less than 25 million years ago (mya), to define the "evolutionary young group" of TEs and KRAB-ZNF genes. This is indeed a valuable recommendation, as 25 mya marks the emergence of Hominoidea (Figure 2C in the manuscript), making it a meaningful reference point for studying recently evolved KRAB-ZNFs and TEs. However, restricting the analysis to elements younger than 25 mya would reduce the dataset to only 9 KRAB-ZNFs and 92 TEs. Nevertheless, we provide here our results for those elements in Table S7:

We observed that among the correlations in the < 25 mya subset, negative correlations (7) outnumbered positive ones (2). However, these correlations were derived from only 3 out of 9 KRAB-ZNFs and 9 out of 92 TE subfamilies. Therefore, based on our data, while the < 25 mya group shows a higher proportion of negative correlations, the sample size is too limited to derive networks or draw robust conclusions in our analysis, especially when compared to our original evolutionary age threshold of 44.2 mya. For this reason, we chose not to reanalyze the data but rather to acknowledge that our current definition of “young” may not be optimal for testing the arms-race model in humans. While previous studies (Jacobs et al., 2014; Bruno et al., 2019; Zuo et al., 2023) have explored relevant KRAB-ZNF and TE interactions, our review of the KRAB-ZNFs and TEs highlighted in those works suggests that a specific focus on elements <25 mya has not been a primary emphasis. 

"our findings only weakly support the arms-race hypothesis. Firstly, we noted that young TEs exhibit lower expression levels than old TEs (Figure 2D and 5B), which might not be expected if they had recently escaped repression". - This is a misinterpretation. These old TEs are no longer harmful. This is not the case of the 'genetic arms race'.

We sincerely appreciate the reviewer’s comments, which have helped us refine our interpretation to prevent potential misunderstandings. Our initial expectation, based on the arms-race hypothesis, was that young TEs would exhibit higher expression levels due to a recent escape from repression, while young KRAB-ZNFs would show increased expression as a counter-adaptive response. However, our findings indicate that both young TEs and young KRAB-ZNFs exhibit lower expression levels. This observation does not align with the classical arms-race model, which typically predicts an ongoing cycle of adaptive upregulation. We rephrase the sentences in our discussion to hopefully make our idea more clear. In addition, we added the notion that older TEs might not be harmful anymore, which we agree with.

"Additionally, some young TEs were also negatively correlated with old KRAB-ZNF genes, leading to weak assortativity regarding age inference, which would also not be in line with the arms-race idea."

This is not a contradiction, as an old KRAB-ZNF gene could be 'reactivated' to protect against young TEs. (It might be cheaper for the host than developing a brand new KRAB-ZNF gene.

We agree with the reviewer's point that older KRAB-ZNFs may be reactivated to suppress young TEs, potentially as a more cost-effective evolutionary strategy than the emergence of entirely new KRAB-ZNFs. We have incorporated this perspective into the revised manuscript to provide a more detailed discussion of our findings.

TEs remain active

In the abstract: "Notably, KRAB-ZNF genes evolve rapidly and exhibit diverse expression patterns in primate brains, where TEs remain active."

This is not precise. TEs are not generally remain active in the brain. It is only the autonomous LINE-1 (young) and non-autonomous Alu (young) and SVA (young) elements that can be mobilized by LINE-1. In addition, the evolutionary young HERV-K is recognized as foreign and alerts the innate immune system (DOI: 10.1172/jci.insight.131093 ) and is a target of the KRAB-ZNF/KAP1 suppression system.

In the abstract: "Evidence indicates that transposable elements (TEs) can contribute to the evolution of new traits, despite often being considered deleterious."

Oversimplification: The harmful and repurposed TEs are washed together.

We appreciate the reviewer’s detailed suggestions for improving the precision of our abstract. While we previously mentioned LINE-1 and Alu elements in the introduction, we now explicitly specify in the abstract that only certain TE subfamilies, such as autonomous LINE-1 and non-autonomous Alu and SVA elements, remain active in the primate brain. Additionally, we have refined the phrasing regarding the role of TEs in evolution to clearly distinguish between their deleterious effects and their potential for functional repurposing. These clarifications have been incorporated into the revised abstract to ensure greater accuracy and nuance.

Positive links

"The high number of positive correlations might be surprising, given that KRAB-ZNFs are considered to repress TEs."

Based on the above, it is not surprising that negative associations are only found with young (< 25 my) TEs. In fact, the relationship between old KRAB-ZNF proteins and old (non-damaging) TEs could be neutral/positive. The case of ZNF528 could be a valuable example of this.

We thank the reviewer for providing this plausible interpretation and added it to the manuscript.

"276 TE:KRAB-ZNF with positive correlations in humans were negatively correlated in bonobos"  It would be important to characterise the positive correlations in more detail. Could it be that the old KRAB-ZNF proteins lost their ability to recruit KAP1/TRIM28? Demonstrate it.

The strategy of developing sequence-specific DNA recognition domains that can specifically recognise TEs is expensive for the host. Recent studies suggest that when the TE is no longer harmful, these proteins/connections can be occasionally repurposed. The repurposed function would probably differ from the original suppressive function.

In my opinion, the TEKRABber tool could be useful in identifying co-option events:

We appreciate the reviewer’s suggestion regarding the characterization of positive correlations. While it is possible that some old KRAB-ZNF proteins have lost their ability to recruit KAP1/TRIM28, we cannot conclude this definitively for all cases. To address this, we examined ChIP-exo data from Imbeault et al. (2017) (Accession: GSE78099) and analyzed the overlap of binding sites between KRAB-ZNFs, KAP1/TRIM28, and RepeatMasker-annotated TEs. Our results indicate that some old KRAB-ZNFs still exhibit binding overlap with KAP1 at TE regions, suggesting that their repressive function may be at least partially retained (Author response image 1).

Author response image 1.<br /> Overlap of KAP1, Zinc finger proteins, and RepeatMasker annotation. Here we detect the overlap of ChIP-exo binding events using KAP1/TRIM28, with KRAB-ZNF genes (one at a time) and RepeatMasker annotation. (115 old and 58 young KRAB-ZNFs, Mann-Whitney, p<0.01).<br />

Minor

"Lead poisoning causes lead ions to compete with zinc ions in zinc finger proteins, affecting proteins such as DNMT1, which are related to the progression of AD (Ordemann and Austin 2016)."

Not precise: While DNMT1 does contain zinc-binding domains, it is not categorized as a zinc finger protein.

We appreciate the reviewer’s insight regarding the classification of DNMT1. After careful consideration, we have removed this sentence from the introduction to maintain focus on KRAB zinc finger proteins.

Definition of TEs

"There were 324 KRAB-ZNFs and 895 TEs expressed in Primate Brain Data." Define it more precisely. It is not clear, what the authors mean by TEs: Are these TE families, subfamilies? Provide information on copy numbers of each in the analysed four species.

We appreciate the reviewer’s suggestion to clarify our definition of TEs. To improve precision, we have specified that the analysis was conducted at the subfamily level. Additionally, we have provided the copy numbers of TEs for the four analyzed species in Table S4.

Occupancy of TEs in the genome

"TEs comprise (i) one third to one half of the mammalian genome and are (ii) not randomly distributed..."

(i) The most accepted number is 45%. However, some more recent reports estimate over 50%, thus the one third is an underestimation.

(ii) Not randomly distributed among the mammalian species?

(i) We thank the reviewer for pointing out that our statement about the abundance of TEs was outdated. We have updated the estimate to reflect that TEs can occupy more than half of the genome, based on recent publications.

(ii) We acknowledge the reviewer’s concern regarding the distribution of TEs. Although TEs are interspersed throughout the genome, their insertion sites are not entirely random, as they tend to exhibit preferences for certain genomic regions. To clarify this, we have revised the wording in the paragraph accordingly.

We would like to express our sincere gratitude to both reviewers for their insightful feedback, which has been instrumental in enhancing the quality of our study.

This was the original docgraph retraction notice.

We did not know it at the time, but the "Update Oct 5 2012" message that we took a screen-shot of would be the only acknowledgement that CMS ever made about the problem.

This is the second post on the DocGraph retraction. It documents the fact that CMS would not collaborate with us to ensure that the data they were releasing was correct, but the data was released and "looked right".

Note: This response was posted by the corresponding author to Review Commons. The content has not been altered except for formatting.

Learn more at Review Commons

Reply to the reviewers

We appreciate the constructive and supportive feedback on our manuscript. All three reviewers acknowledged the significance and novelty of our work on bacterial telomere protection. In response to their suggestions, we have conducted the requested experiments and revised the manuscript accordingly. These changes have enhanced the rigor of our study and clarified our interpretations and explanations.

Moreover, we characterized an additional truncation mutant of TelN (TelN Δ445–631), which lacks the two C-terminal domains. Despite this deletion, the mutant retained protection activity (Supplementary Figure S4B), indicating that the remaining regions of the protein are sufficient to confer efficient protection in this assay.

Finally, we removed three sequence alignments (previously Supplementary Figures S6A and S7), as we recognized that the high degree of sequence divergence could hinder proper alignment and potentially lead to misinterpretation.

Reviewer #1 (Evidence, reproducibility and clarity (Required)):

This study addresses how the bacterial telomere protein TelN protects telomere ends against the action of the Mre11-Rad50 nuclease (MR). This protection is essential for the stability of hairpin-ended linear plasmid and chromosomes in bacteria but had not been explored before. The authors demonstrate that TelN is necessary and sufficient to block MR-dependent DNA cleavage when bound to its specific telomere sequence. By combining elegant genetics and biochemical approaches, it convincingly shows that TelN-dependent inhibition likely involves a specific interaction between TelN and the MR complex. The manuscript is well written, easy to read and focused on the relevant information. The claims and the conclusions are supported by the data. There is no over-interpretation.

Comments: - Figure 1B, unnormalized transformation efficiency would be useful to show in SI

The unnormalized B. subtilis transformation efficiency has now been added as new figure panel S1B.

• Figures 2B, 2C, 3C, 3D, 4C, 5A and 5B: quantification of independent experiments should be added

While these DNA protection experiments show a clearly reproducible pattern of DNA degradation, the exact response to TelN titration varies somewhat between experimental replicates. We initially included the quantification of remaining full-length DNA because the corresponding band is hard to discern in the gel image due to pixel saturation. However, we realize now that this may mislead readers to think that the degradation occurs always with the exact same dosage response.

To avoid this, we have decided to remove the quantification and instead show the relevant part of the gel also at higher contrast to better visualize the loss of full-length DNA due to DNA degradation. In addition, we have included replicate experiments carried out at the same MR concentration (125 nM M₂R₂) or at higher concentration (500 nM M₂R₂) in the supplementary material. These examples demonstrate the general reproducibility of the assay.

**Referee cross-commenting**

Perfect for me. It seems that there is a consensus.

Reviewer #1 (Significance (Required)):

This pioneering study provides a very strong basis for a new understanding of telomeres in bacteria and offers fascinating evolutionary perspectives when compared to similar mechanisms active at telomeres in eukaryotic cells.

Reviewer #2 (Evidence, reproducibility and clarity (Required)):

The paper is well-presented and well-written throughout. The paper shows convincingly that TelN protects hairpin DNA ends from the activity of SbcCD, presumably providing a protection mechanism for N15 phage DNA in vivo. Furthermore, this protection activity is shown not to require the catalytic (resolvase) activity of TelN, nor its poorly characterised C-terminal domain. The paper also suggests that this inhibition acts both at the level of competition for the DNA hairpin end and at the level of a direct protein:protein interaction between TelN and MR. An (acknowledged) weakness is that there is no real insight into the protein:protein interaction suggested by the experiments shown in Figure 5. Ideally, the protein:protein interaction interface would be identified and mutations in this interface would be shown to reduce hairpin protection.

Specific comments/questions

(1) What pathway (in vivo) leads to inactivation of linear hairpin DNA - one suspects that cleavage by SbcCD at the hairpins is probably not the full story. Presumably SbcCD cleavage facilitates further processing by other long range resection systems such as RecBCD, Exo1, RecQ/J etc. Would it be appropriate to view the hairpin as an adaption to protect against these nucleases, which then must be complemented with a mechanism to suppress SbcCD?

The reviewer's suggestion that hairpin ends represent a first layer of adaptation against nucleolytic processing is compelling. Hairpin structures inherently resist many exonucleases due to their covalently closed nature (absence of free 3’ or 5’ ends) but remain vulnerable to MR processing (Connelly et al, 1998, 1999; Saathoff et al, 2018). This creates a scenario where effective telomere protection requires both the structural barrier provided by the hairpin and an active mechanism to suppress MR activity. We have added this perspective to the relevant paragraph in the discussion.

(2) Section starting "Direct inhibition of MR by TelN in vitro". What is the word direct supposed to convey here? To me it suggests that the inhibition is via direct interaction of TelN with MR (rather than, for example, a result of competition for the hairpin DNA end) which is not shown here. Suggest either defining or removing the word direct. This point gains more importance considering that differentiating between inhibition mechanisms becomes a focus of later parts of the paper.

By "direct inhibition," we meant that TelN blocks MR nuclease activity without requiring additional cofactors, as demonstrated in this minimal reaction system containing only TelN, MR complex, DNA substrate, and ATP. To avoid ambiguity, we have reworded the corresponding headline and paragraph.

(3) Figure 2B - Why no control lane without MR? - this is a basic control to show that he degradation we are seeing in the absence of TelN is MR-dependent. Formally, as shown, the degradation could be caused by the ATP stock.

We have now included ATP-only control lanes (without MR complex), which show no substrate degradation, confirming that ATP stocks do not contain contaminating nucleases and that the observed degradation is indeed MR-dependent. These controls are included in the supplementary data (Figure S3A) along with additional replicate experiments. Notably, the dose-dependent protection observed at low TelN concentrations (where MR activity is not fully inhibited) provides additional evidence for the specificity of the MR-TelN interaction system, as non-specific nuclease contamination would result in complete substrate degradation regardless of TelN concentration.

(4) Why not use B. subtilis SbcCD for the species specificity experiment? Also, is it not surprising that TelN yielded zero protection against MRX given that the DNA sequence specificity experiments above suggest competition for DNA substrate is part of the inhibition mechanism?

We agree that this would be a great addition. We attempted but were unable to purify active B. subtilis SbcCD protein despite multiple attempts. The yeast MRX experiment serves the same purpose of demonstrating species specificity and represents a more evolutionarily distant comparison, which strengthens our conclusions about bacterial-specific inhibition.

(5) If the authors felt it appropriate, I thought there was scope for further discussion/introductory material. There are strong parallels here with mechanisms used by phage to protect themselves from the activities of RecBCD, which include both proteins that protect DNA ends like T4 gene 2, we well as proteins that bind directly to RecBCD to inactivate it like lambda Gam. As such, the work here will appeal as much to those interested in bacterial defence systems / phage:host interactions as it does to those interested in telomere biology. Especially significant is the inhibition of DNA end processing factors by lambda Gam since this protein is reported to interact with both RecBCD and SbcCD (PMID: 2531105).

We agree that there are obvious parallels between lambda Gam and TelN as counter-defence factors. This was likely largely missed in previous work because the telomere resolution activity of TelN masked its function in counter-defence. We have added a statement on this matter at the end of the discussion.

(6) Just a gripe really: it seems to be 'de rigeur' at the moment to re-name bacterial proteins for their human orthologues, presumably to elevate the perceived importance of the work(?), but it is not a practice I think is terribly helpful as it causes issues when searching literature. Minimally it would be great if the authors could ensure they add SbcCD as a keyword for search purposes.

We appreciate the reviewer's concern about nomenclature inconsistencies in the literature. We have chosen MR over SbcCD as a more generic term that covers eukaryotes, archaea and lately also bacteria and will hopefully contribute to a more consistent terminology in the literature across the domains of life in the future. Our choice to use "Mre11-Rad50" (MR) for the E. coli SbcCD complex is also consistent with prominent recent publications (Käshammer et al., 2019; Gut et al., 2022), explicitly referring to the E. coli system as "Mre11-Rad50" while acknowledging the bacterial designation. To link to previous literature, we made sure that both "SbcCD" and "Mre11-Rad50" are mentioned in the abstract. And, as suggested, we have now also added “SbcCD” to our keyword list to facilitate comprehensive literature searches.

**Referee cross-commenting**

I have nothing to add. The reviewers' comments are all broadly positive and consistent.

Reviewer #2 (Significance (Required):

This is an excellent paper unveiling a phage encoded "counter-defence" mechanism designed to protect phage DNA from degradation. It will be of special interest to those studying telomere biology of phage:host interactions.

---

---

Reviewer #3

The authors investigate how the N15 phage protelomerase TelN protects linear chromosomes that terminate in hairpin structures (a sort of telomere). In E. coli and B. subtilis cells, removal or truncation of telN reduces transformation/survival of linear DNA, whereas complementation with full-length or a catalytically inactive TelN restores viability, consistent with TelN playing a nonenzymatic capping function.

In vitro, TelN binds hairpin substrates with moderate affinity and protects them from the nuclease activity of the Mre11/Rad50 complex. The authors propose that TelN originated as an early, sequence specific barrier against MR mediated DNA end processing, establishing fundamental principles of telomere protection that persist from bacteria to eukaryotes.

Major comments:

The manuscript convincingly shows that TelN can functionally block the Mre11Rad50 (MR) nuclease on a hairpin DNA end in a sequence specific manner (suggesting a physical interaction), but it doesn't directly demonstrate this. A simple pull-down or equilibrium binding method would be useful in proving a physical interaction.

We agree that this would be a valuable addition to the study. We have made several attempts to detect direct interaction by co-immunoprecipitation. However, without success so far. We do not have sufficient material for equilibrium binding methods (yet).__ ____ __

The MR complex requires ATP hydrolysis for resection of DNA ends. It would be a nice addition to the manuscript if the effect of TelN of Rad50 ATPase activity was tested.

We have tested the effect of TelN on Rad50 ATPase activity and found no significant impact under our experimental conditions, possible in line with the lack of stable interaction.

The bar plot on Fig 3B indicates that the experiments are performed in triplicate. The statistical significance of the differences between conditions should be determined. The same general comment could be made regarding the quantification of the polyacrylamide gels - how reproducible are these values?

We performed paired t-test analysis for the following figures and now indicate the p-values wherever significant (below 0.05): Figures 1D, 1E, 3B, 4B and S4B. We used paired t-tests to generally compare linear vs circular plasmid transformation efficiency for each condition. In Figure 4B, which included two different linear DNA constructs, we compared the two linear DNA constructs directly to each other. [Given that our experimental design included multiple control conditions with known expected outcomes to validate assay performance, rather than many independent exploratory comparisons, we report uncorrected p-values as the primary analysis. The inclusion of multiple controls with predictable outcomes reduces the likelihood of false positive interpretations.]

As stated in response to reviewer 1, while the exact values for the DNA degradation profile vary somewhat between experiments (likely due to variations in band quantification – see also response to comment below), the general trends are robust as for example indicated by similar experiments performed with higher MR concentration (500 nM instead of 125 nM M₂R₂ concentrations for all TelN variants) demonstrating reproducibility across different conditions. For Figure 5, however, we are unable to provide additional repeat experiments due to limitations in reagent availability. Considering the robust effect seen with Ec MR controls and the presence of multiple samples in the dilution series, we are nevertheless confident about the conclusion.

Minor comments:

A better explanation of how the gels were quantified should be provided. Were the products included in the analysis, or was it just the decrease in the substrate band that was measured?

As also stated above, we have removed the band quantification and instead show the bands also at different contrast settings.

In our original approach, gel band quantification was performed using ImageQuant TL software (version 8.2.0, GE Healthcare). For each gel, individual lanes were defined using either fixed-width boundaries (95-103 pixels) or automatic edge detection, depending on the gel quality and band definition. Band volumes were calculated using rolling ball background subtraction (radius 180 pixels) with automatic band detection. Substrate degradation was assessed by measuring the integrated density (volume) of the remaining full-length (or near full-length) substrate bands under different treatment conditions. The band volume values were plotted directly to compare substrate levels across treatment groups.

We now present the data as two gel panels: an exposure showing the full reaction profile, and another exposure focusing on the substrate bands to clearly demonstrate dose-dependent protection. Additional replicate experiments including ATP-only controls (confirming no contamination from ATP stocks) and experiments at 500 nM M₂R₂ concentrations, are provided in the supplementary data. This approach provides more direct visualization of the biological phenomenon with comprehensive control validation.

I felt like the Results jump rather abruptly from B. subtilis chromosome assays to E. coli plasmid experiments. Maybe the addition of a few linking sentences would improve this transition.

---

Upon re-reading the manuscript we agree with this assertion and have added further information to provide a smoother transition.

A comment on the stoichiometry of TelN and genome ends during phage replication would be useful.

Our in vitro data suggest that effective protection can be achieved at relatively low TelN:DNA ratios in vitro, consistent with the notion of formation of stable, protective nucleoprotein structures. We unfortunately do not currently have information on the copy number of TelN per cell or per hairpin end. It is not easy to obtain reliable values for these numbers. However, we can speculate that multiple TelN proteins are present due to the presence of three copies of a DNA sequence motif (binding to CTD1) in each telomeric DNA, consistent with the formation of stable, protective nucleoprotein structures.

Reviewer #3 (Significance (Required)):

General assessment:

Strengths: A nice combination of genetics and biochemistry convincingly demonstrates that TelN protects linear chromosomes/replicons from MR-dependent degradation independent of its cleavage-ligase activity. It does this by binding to the hairpin DNA ends in a sequence specific fashion and the species specificity suggests a direct physical interaction, which likely inhibits the nuclease activity of the MR complex

Limitations: The lack of characterization of the putative physical interaction between TelN and the MR complex is considered a weakness.

Advance: The manuscript fills in a mechanistic gap between protelomerase-mediated telomere formation and maintenance by demonstrating a protective/capping role. This is the first quantitative analysis of DNA-end protection from MR nuclease activity by TelN.

Audience: Readers interested in bacterial chromosome biology, DNA repair, the parallels to eukaryotic shelterin will be interesting to the broader telomere and genome stability communities.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #3

Evidence, reproducibility and clarity

Summary:

The authors investigate how the N15 phage protelomerase TelN protects linear chromosomes that terminate in hairpin structures (a sort of telomere). In E. coli and B. subtilis cells, removal or truncation of telN reduces transformation/survival of linear DNA, whereas complementation with full‑length or a catalytically inactive TelN restores viability, consistent with TelN playing a non‑enzymatic capping function.

In vitro, TelN binds hairpin substrates with  moderate affinity and protects them from the nuclease activity of the Mre11/Rad50 complex. The authors propose that TelN originated as an early, sequence‑specific barrier against MR‑mediated DNA end processing, establishing fundamental principles of telomere protection that persist from bacteria to eukaryotes.

Major comments:

The manuscript convincingly shows that TelN can functionally block the Mre11‑Rad50 (MR) nuclease on a hair‑pin DNA end in a sequence specific manner (suggesting a physical interaction), but it doesn't directly demonstrate this. A simple pull-down or equilibrium binding method would useful in proving a physical interaction.

The MR complex requires ATP hydrolysis for resection of DNA ends. It would be a nice addition to the manuscript if the effect of TelN of Rad50 ATPase activity was tested.

The bar plot on Fig 3B indicates that the experiments are performed in triplicate. The statistical significance of the differences between conditions should be determined. The same general comment could be made regarding the quantification of the polyacrylamide gels - how reproducible are these values?

Minor comments:

A better explanation of how the gels were quantified should be provided. Were the products included in the analysis, or was it just the decrease in the substrate band that was measured?

I felt like the Results jump rather abruptly from B. subtilis chromosome assays to E. coli plasmid experiments. Maybe the addition of a few linking sentences would improve this transition.

A comment on the stoichiometry of TelN and genome ends during phage replication would be useful.

Significance

General assessment:

Strengths: A nice combination of genetics and biochemistry convincingly demonstrates that TelN protects linear chromosomes/replicons from MR-dependent degradation independent of its cleavage-ligase activity. It does this by binding to the hairpin DNA ends in a sequence specific fashion and the species specificity suggests a direct physical interaction, which likely inhibits the nuclease activity of the MR complex

Limitations: The lack of characterization of the putative physical interaction between TelN and the MR complex is considered a weakness.

Advance: The manuscript fills in a mechanistic gap between protelomerase‑mediated telomere formation and maintenance by demonstrating a protective/capping role. This is the first quantitative analysis of DNA-end protection from MR nuclease activity by TelN.

Audience: Readers interested in bacterial chromosome biology, DNA repair, the parallels to eukaryotic shelterin will be interesting to the broader telomere and genome‑stability communities.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #2

Evidence, reproducibility and clarity

The paper is well-presented and well-written throughout. The paper shows convincingly that TelN protects hairpin DNA ends from the activity of SbcCD, presumably providing a protection mechanism for N15 phage DNA in vivo. Furthermore, this protection activity is shown not to require the catalytic (resolvase) activity of TelN, nor its poorly characterised C-terminal domain. The paper also suggests that this inhibition acts both at the level of competition for the DNA hairpin end and at the level of a direct protein:protein interaction between TelN and MR. An (acknowledged) weakness is that there is no real insight into the protein:protein interaction suggested by the experiments shown in Figure 5. Ideally, the protein:protein interaction interface would be identified and mutations in this interface would be shown to reduce hairpin protection.

Specific comments/questions

(1) What pathway (in vivo) leads to inactivation of linear hairpin DNA - one suspects that cleavage by SbcCD at the hairpins is probably not the full story. Presumably SbcCD cleavage facilitates further processing by other long range resection systems such as RecBCD, Exo1, RecQ/J etc. Would it be appropriate to view the hairpin as an adaption to protect against these nucleases, which then must be complemented with a mechanism to suppress SbcCD?

(2) Section starting "Direct inhibition of MR by TelN in vitro". What is the word direct supposed to convey here? To me it suggests that the inhibition is via direct interaction of TelN with MR (rather than, for example, a result of competition for the hairpin DNA end) which is not shown here. Suggest either defining or removing the word direct. This point gains more importance considering that differentiating between inhibition mechanisms becomes a focus of later parts of the paper.

(3) Figure 2B - Why no control lane without MR? - this is a basic control to show that he degradation we are seeing in the absence of TelN is MR-dependent. Formally, as shown, the degradation could be caused by the ATP stock.

(4) Why not use B. subtilis SbcCD for the species specificity experiment? Also, is it not surprising that TelN yielded zero protection against MRX given that the DNA sequence specificity experiments above suggest competition for DNA substrate is part of the inhibition mechanism?

(5) If the authors felt it appropriate, I thought there was scope for further discussion/introductory material. There are strong parallels here with mechanisms used by phage to protect themselves from the activities of RecBCD, which include both proteins that protect DNA ends like T4 gene 2, we well as proteins that bind directly to RecBCD to inactivate it like lambda Gam. As such, the work here will appeal as much to those interested in bacterial defence systems / phage:host interactions as it does to those interested in telomere biology. Especially significant is the inhibition of DNA end processing factors by lambda Gam since this protein is reported to interact with both RecBCD and SbcCD (PMID: 2531105).

(6) Just a gripe really: it seems to be 'de rigeur' at the moment to re-name bacterial proteins for their human orthologues, presumably to elevate the perceived importance of the work(?), but it is not a practice I think is terribly helpful as it causes issues when searching literature. Minimally it would be great if the authors could ensure they add SbcCD as a keyword for search purposes.

Referee cross-commenting

I have nothing to add. The reviewers comments are all broadly positive and and consistent.

Significance

This is an excellent paper unveiling a phage encoded "counter-defence" mechanism designed to protect phage DNA from degradation. It will be of special interest to those studying telomere biology of phage:host interactions.

Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

Learn more at Review Commons

Referee #1

Evidence, reproducibility and clarity

This study addresses how the bacterial telomere protein TelN protect telomere ends against the action of the Mre11-Rad50 nuclease (MR). This protection is essential for the stability of hairpin-ended linear plasmid and chromosomes in bacteria but had not been explored before. The authors demonstrates that TelN is necessary and sufficient to block MR-dependent DNA cleavage when bound to its specific telomere sequence. By combining elegant genetics and biochemical approaches, it convincingly shows that TelN-dependent inhibition likely involves a specific interaction between TelN and the MR complex. The manuscript is well written, easy to read and focused on the relevant information. The claims and the conclusions are supported by the data. There is no over-interpretation.

Comments:

• Figure 1B, unnormalized transformation efficiency would be useful to show in SI
• Figures 2B, 2C, 3C, 3D, 4C, 5A and 5B: quantification of independent experiments should be added

Referee cross-commenting

Perfect for me. It seems that there is a consensus.

Significance

This pioneering study provides a very strong basis for a new understanding of telomeres in bacteria and offers fascinating evolutionary perspectives when compared to similar mechanisms active at telomeres in eukaryotic cells.

eLife Assessment

This study provides valuable insights into the evolutionary conservation of sex determination mechanisms in ants by identifying a candidate sex-determining region in a parthenogenetic species. The strength of evidence is solid, using well-executed genomic analyses to identify differences in heterozygosity between females and diploid males, though not yet functional validation of the candidate locus.

Reviewer #1 (Public review):

The authors have implemented several clarifications in the text and improved the connection between their findings and previous work. As stated in my initial review, I had no major criticisms of the previous version of the manuscript, and I continue to consider this a solid and well-written study. However, the revised manuscript still largely reiterates existing findings and does not offer novel conceptual or experimental advances. It supports previous conclusions suggesting a likely conserved sex determination locus in aculeate hymenopterans, but does so without functional validation (i.e., via experimental manipulation) of the candidate locus in O. biroi. I also wish to clarify that I did not intend to imply that functional assessments in the Pan et al. study were conducted in more than one focal species; my previous review explicitly states that the locus's functional role was validated in the Argentine ant.

Reviewer #3 (Public review):

The authors have made considerable efforts to conduct functional analyses to the fullest extent possible in this study; however, it is understandable that meaningful results have not yet been obtained. In the revised version, they have appropriately framed their claims within the limits of the current data and have adjusted their statements as needed in response to the reviewers' comments.

Author response:

The following is the authors’ response to the original reviews

Reviewer #1 (Public review):

This study investigates the sex determination mechanism in the clonal ant Ooceraea biroi, focusing on a candidate complementary sex determination (CSD) locus-one of the key mechanisms supporting haplodiploid sex determination in hymenopteran insects. Using whole genome sequencing, the authors analyze diploid females and the rarely occurring diploid males of O. biroi, identifying a 46 kb candidate region that is consistently heterozygous in females and predominantly homozygous in diploid males. This region shows elevated genetic diversity, as expected under balancing selection. The study also reports the presence of an lncRNA near this heterozygous region, which, though only distantly related in sequence, resembles the ANTSR lncRNA involved in female development in the Argentine ant, Linepithema humile (Pan et al. 2024). Together, these findings suggest a potentially conserved sex determination mechanism across ant species. However, while the analyses are well conducted and the paper is clearly written, the insights are largely incremental. The central conclusion - that the sex determination locus is conserved in ants - was already proposed and experimentally supported by Pan et al. (2024), who included O. biroi among the studied species and validated the locus's functional role in the Argentine ant. The present study thus largely reiterates existing findings without providing novel conceptual or experimental advances.

Although it is true that Pan et al., 2024 demonstrated (in Figure 4 of their paper) that the synteny of the region flanking ANTSR is conserved across aculeate Hymenoptera (including O. biroi), Reviewer 1’s claim that that paper provides experimental support for the hypothesis that the sex determination locus is conserved in ants is inaccurate. Pan et al., 2024 only performed experimental work in a single ant species (Linepithema humile) and merely compared reference genomes of multiple species to show synteny of the region, rather than functionally mapping or characterizing these regions.

Other comments:

The mapping is based on a very small sample size: 19 females and 16 diploid males, and these all derive from a single clonal line. This implies a rather high probability for false-positive inference. In combination with the fact that only 11 out of the 16 genotyped males are actually homozygous at the candidate locus, I think a more careful interpretation regarding the role of the mapped region in sex determination would be appropriate. The main argument supporting the role of the candidate region in sex determination is based on the putative homology with the lncRNA involved in sex determination in the Argentine ant, but this argument was made in a previous study (as mentioned above).

Our main argument supporting the role of the candidate region in sex determination is not based on putative homology with the lncRNA in L. humile. Instead, our main argument comes from our genetic mapping (in Fig. 2), and the elevated nucleotide diversity within the identified region (Fig. 4). Additionally, we highlight that multiple genes within our mapped region are homologous to those in mapped sex determining regions in both L. humile and Vollenhovia emeryi, possibly including the lncRNA.

In response to the Reviewer’s assertion that the mapping is based on a small sample size from a single clonal line, we want to highlight that we used all diploid males available to us. Although the primary shortcoming of a small sample size is to increase the probability of a false negative, small sample sizes can also produce false positives. We used two approaches to explore the statistical robustness of our conclusions. First, we generated a null distribution by randomly shuffling sex labels within colonies and calculating the probability of observing our CSD index values by chance (shown in Fig. 2). Second, we directly tested the association between homozygosity and sex using Fisher’s Exact Test (shown in Supplementary Fig. S2). In both cases, the association of the candidate locus with sex was statistically significant after multiple-testing correction using the Benjamini-Hochberg False Discovery Rate. These approaches are clearly described in the “CSD Index Mapping” section of the Methods.

We also note that, because complementary sex determination loci are expected to evolve under balancing selection, our finding that the mapped region exhibits a peak of nucleotide diversity lends orthogonal support to the notion that the mapped locus is indeed a complementary sex determination locus.

The fourth paragraph of the results and the sixth paragraph of the discussion are devoted to explaining the possible reasons why only 11/16 genotyped males are homozygous in the mapped region. The revised manuscript will include an additional sentence (in what will be lines 384-388) in this paragraph that includes the possible explanation that this locus is, in fact, a false positive, while also emphasizing that we find this possibility to be unlikely given our multiple lines of evidence.

In response to Reviewer 1’s suggestion that we carefully interpret the role of the mapped region in sex determination, we highlight our careful wording choices, nearly always referring to the mapped locus as a “candidate sex determination locus” in the title and throughout the manuscript. For consistency, the revised manuscript version will change the second results subheading from “The O. biroi CSD locus is homologous to another ant sex determination locus but not to honeybee csd” to “O. biroi’s candidate CSD locus is homologous to another ant sex determination locus but not to honeybee csd,” and will add the word “candidate” in what will be line 320 at the beginning of the Discussion, and will change “putative” to “candidate” in what will be line 426 at the end of the Discussion.

In the abstract, it is stated that CSD loci have been mapped in honeybees and two ant species, but we know little about their evolutionary history. But CSD candidate loci were also mapped in a wasp with multi-locus CSD (study cited in the introduction). This wasp is also parthenogenetic via central fusion automixis and produces diploid males. This is a very similar situation to the present study and should be referenced and discussed accordingly, particularly since the authors make the interesting suggestion that their ant also has multi-locus CSD and neither the wasp nor the ant has tra homologs in the CSD candidate regions. Also, is there any homology to the CSD candidate regions in the wasp species and the studied ant?

In response to Reviewer 1’s suggestion that we reference the (Matthey-Doret et al. 2019) study in the context of diploid males being produced via losses of heterozygosity during asexual reproduction, the revised manuscript will include (in what will be lines 123-126) the highlighted portion of the following sentence: “Therefore, if O. biroi uses CSD, diploid males might result from losses of heterozygosity at sex determination loci (Fig. 1C), similar to what is thought to occur in other asexual Hymenoptera that produce diploid males (Rabeling and Kronauer 2012; Matthey-Doret et al. 2019).”

We note, however, that in their 2019 study, Matthey-Doret et al. did not directly test the hypothesis that diploid males result from losses of heterozygosity at CSD loci during asexual reproduction, because the diploid males they used for their mapping study came from inbred crosses in a sexual population of that species.

We address this further below, but we want to emphasize that we do not intend to argue that O. biroi has multiple CSD loci. Instead, we suggest that additional, undetected CSD loci is one possible explanation for the absence of diploid males from any clonal line other than clonal line A. In response to Reviewer 1’s suggestion that we reference the (Matthey-Doret et al. 2019) study in the context of multilocus CSD, the revised manuscript version will include the following additional sentence in the fifth paragraph of the discussion (in what will be lines 372-374): “Multi-locus CSD has been suggested to limit the extent of diploid male production in asexual species under some circumstances (Vorburger 2013; Matthey-Doret et al. 2019).”

Regarding Reviewer 2’s question about homology between the putative CSD loci from the (Matthey-Doret et al. 2019) study and O. biroi, we note that there is no homology. The revised manuscript version will have an additional Supplementary Table (which will be the new Supplementary Table S3) that will report the results of this homology search. The revised manuscript will also include the following additional sentence in the Results, in what will be lines 172-174: “We found no homology between the genes within the O. biroi CSD index peak and any of the genes within the putative L. fabarum CSD loci (Supplementary Table S3).”

The authors used different clonal lines of O. biroi to investigate whether heterozygosity at the mapped CSD locus is required for female development in all clonal lines of O. biroi (L187-196). However, given the described parthenogenesis mechanism in this species conserves heterozygosity, additional females that are heterozygous are not very informative here. Indeed, one would need diploid males in these other clonal lines as well (but such males have not yet been found) to make any inference regarding this locus in other lines.

We agree that a full mapping study including diploid males from all clonal lines would be preferable, but as stated earlier in that same paragraph, we have only found diploid males from clonal line A. We stand behind our modest claim that “Females from all six clonal lines were heterozygous at the CSD index peak, consistent with its putative role as a CSD locus in all O. biroi.” In the revised manuscript version, this sentence (in what will be lines 199-201) will be changed slightly in response to a reviewer comment below: “All females from all six clonal lines (including 26 diploid females from clonal line B) were heterozygous at the CSD index peak, consistent with its putative role as a CSD locus in all O. biroi.”

Reviewer #2 (Public review):

The manuscript by Lacy et al. is well written, with a clear and compelling introduction that effectively conveys the significance of the study. The methods are appropriate and well-executed, and the results, both in the main text and supplementary materials, are presented in a clear and detailed manner. The authors interpret their findings with appropriate caution.

This work makes a valuable contribution to our understanding of the evolution of complementary sex determination (CSD) in ants. In particular, it provides important evidence for the ancient origin of a non-coding locus implicated in sex determination, and shows that, remarkably, this sex locus is conserved even in an ant species with a non-canonical reproductive system that typically does not produce males. I found this to be an excellent and well-rounded study, carefully analyzed and well contextualized.

That said, I do have a few minor comments, primarily concerning the discussion of the potential 'ghost' CSD locus. While the authors acknowledge (line 367) that they currently have no data to distinguish among the alternative hypotheses, I found the evidence for an additional CSD locus presented in the results (lines 261-302) somewhat limited and at times a bit difficult to follow. I wonder whether further clarification or supporting evidence could already be extracted from the existing data. Specifically:

We agree with Reviewer 2 that the evidence for a second CSD locus is limited. In fact, we do not intend to advocate for there being a second locus, but we suggest that a second CSD locus is one possible explanation for the absence of diploid males outside of clonal line A. In our initial version, we intentionally conveyed this ambiguity by titling this section “O. biroi may have one or multiple sex determination loci.” However, we now see that this leads to undue emphasis on the possibility of a second locus. In the revised manuscript, we will split this into two separate sections: “Diploid male production differs across O. biroi clonal lines” and “O. biroi lacks a tra-containing CSD locus.”

(1) Line 268: I doubt the relevance of comparing the proportion of diploid males among all males between lines A and B to infer the presence of additional CSD loci. Since the mechanisms producing these two types of males differ, it might be more appropriate to compare the proportion of diploid males among all diploid offspring. This ratio has been used in previous studies on CSD in Hymenoptera to estimate the number of sex loci (see, for example, Cook 1993, de Boer et al. 2008, 2012, Ma et al. 2013, and Chen et al., 2021). The exact method might not be applicable to clonal raider ants, but I think comparing the percentage of diploid males among the total number of (diploid) offspring produced between the two lineages might be a better argument for a difference in CSD loci number.

We want to re-emphasize here that we do not wish to advocate for there being two CSD loci in O. biroi. Rather, we want to explain that this is one possible explanation for the apparent absence of diploid males outside of clonal line A. We hope that the modifications to the manuscript described in the previous response help to clarify this.

Reviewer 2 is correct that comparing the number of diploid males to diploid females does not apply to clonal raider ants. This is because males are vanishingly rare among the vast numbers of females produced. We do not count how many females are produced in laboratory stock colonies, and males are sampled opportunistically. Therefore, we cannot report exact numbers. However, we will add the highlighted portion of the following sentence (in what will be lines 268-270) to the revised manuscript: “Despite the fact that we maintain more colonies of clonal line B than of clonal line A in the lab, all the diploid males we detected came from clonal line A.”

(2) If line B indeed carries an additional CSD locus, one would expect that some females could be homozygous at the ANTSR locus but still viable, being heterozygous only at the other locus. Do the authors detect any females in line B that are homozygous at the ANTSR locus? If so, this would support the existence of an additional, functionally independent CSD locus.

We thank the reviewer for this suggestion, and again we emphasize that we do not want to argue in favor of multiple CSD loci. We just want to introduce it as one possible explanation for the absence of diploid males outside of clonal line A.

The 26 sequenced diploid females from clonal line B are all heterozygous at the mapped locus, and the revised manuscript will clarify this in what will be lines 199-201. Previously, only six of those diploid females were included in Supplementary Table S2, and that will be modified accordingly.

(3) Line 281: The description of the two tra-containing CSD loci as "conserved" between Vollenhovia and the honey bee may be misleading. It suggests shared ancestry, whereas the honey bee csd gene is known to have arisen via a relatively recent gene duplication from fem/tra (10.1038/nature07052). It would be more accurate to refer to this similarity as a case of convergent evolution rather than conservation.

In the sentence that Reviewer 2 refers to, we are representing the assertion made in the (Miyakawa and Mikheyev 2015) paper in which, regarding their mapping of a candidate CSD locus that contains two linked tra homologs, they write in the abstract: “these data support the prediction that the same CSD mechanism has indeed been conserved for over 100 million years.” In that same paper, Miyakawa and Mikheyev write in the discussion section: “As ants and bees diverged more than 100 million years ago, sex determination in honey bees and V. emeryi is probably homologous and has been conserved for at least this long.”

As noted by Reviewer 2, this appears to conflict with a previously advanced hypothesis: that because fem and csd were found in Apis mellifera, Apis cerana, and Apis dorsata, but only fem was found in Mellipona compressipes, Bombus terrestris, and Nasonia vitripennis, that the csd gene evolved after the honeybee (Apis) lineage diverged from other bees (Hasselmann et al. 2008). However, it remains possible that the csd gene evolved after ants and bees diverged from N. vitripennis, but before the divergence of ants and bees, and then was subsequently lost in B. terrestris and M. compressipes. This view was previously put forward based on bioinformatic identification of putative orthologs of csd and fem in bumblebees and in ants [(Schmieder et al. 2012), see also (Privman et al. 2013)]. However, subsequent work disagreed and argued that the duplications of tra found in ants and in bumblebees represented convergent evolution rather than homology (Koch et al. 2014). Distinguishing between these possibilities will be aided by additional sex determination locus mapping studies and functional dissection of the underlying molecular mechanisms in diverse Aculeata.

Distinguishing between these competing hypotheses is beyond the scope of our paper, but the revised manuscript will include additional text to incorporate some of this nuance. We will include these modified lines below (in what will be lines 287-295), with the additions highlighted:

“A second QTL region identified in V. emeryi (V.emeryiCsdQTL1) contains two closely linked tra homologs, similar to the closely linked honeybee tra homologs, csd and fem (Miyakawa and Mikheyev 2015). This, along with the discovery of duplicated tra homologs that undergo concerted evolution in bumblebees and ants (Schmieder et al. 2012; Privman et al. 2013) has led to the hypothesis that the function of tra homologs as CSD loci is conserved with the csd-containing region of honeybees (Schmieder et al. 2012; Miyakawa and Mikheyev 2015). However, other work has suggested that tra duplications occurred independently in honeybees, bumblebees, and ants (Hasselmann et al. 2008; Koch et al. 2014), and it remains to be demonstrated that either of these tra homologs acts as a primary CSD signal in V. emeryi.”

(4) Finally, since the authors successfully identified multiple alleles of the first CSD locus using previously sequenced haploid males, I wonder whether they also observed comparable allelic diversity at the candidate second CSD locus. This would provide useful supporting evidence for its functional relevance.

As is already addressed in the final paragraph of the results and in Supplementary Fig. S4, there is no peak of nucleotide diversity in any of the regions homologous to V.emeryiQTL1, which is the tra-containing candidate sex determination locus (Miyakawa and Mikheyev 2015). In the revised manuscript, the relevant lines will be 307-310. We want to restate that we do not propose that there is a second candidate CSD locus in O. biroi, but we simply raise the possibility that multi-locus CSD *might* explain the absence of diploid males from clonal lines other than clonal line A (as one of several alternative possibilities).