DOI: 10.1038/s41593-025-01999-y

Resource: (IMSR Cat# JAX_008199,RRID:IMSR_JAX:008199)

Curator: @scibot

SciCrunch record: RRID:IMSR_JAX:008199

What is this?

DOI: 10.1016/j.devcel.2025.07.010

Resource: (Abcam Cat# ab124964, RRID:AB_11129103)

Curator: @scibot

SciCrunch record: RRID:AB_11129103

What is this?

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

REVIEWER 1

This is an important and solid study that identified sequences that can improve circRNA translation and that as or more importantly are very short and hence are suitable for generating of efficient protein expressing circRNAs. This manuscript fills an important gap in the field, and it is highly significant. The study is well controlled, the rationale clear and the results conclusive with no major flaws.

• While this is a minor concern as the vector has been used before, it will greatly improve the quality of the paper if the authors could just verify that the vector only generates circRNA molecules and not linear concatenamers. To do so the authors can focus only in their control and the most optimal transcripts and perform northern blot or well controlled RNAseR experiments to show that all RNA molecules containing the back splicing junction are circular We thank the reviewer for raising this point. As suggested, we performed RNaseR resistance assays on our three most efficient candidates driving cGFP translation (VCIP, T3-glo, and T3-U3) to confirm that all derived RNA molecules containing the back-splicing junction are circular. As proof of this, cGFP proved strongly resistant to RNase R (new Fig. S1N), confirming its circular structure. We further ruled out the possibility that molecules other than the circRNA encoding GFP serve as templates for translation from our vectors. Specifically, ad hoc PCR amplifications performed for this purpose (new Fig. S1M) showed no bands that would indicate the presence of concatemers. Indeed, ad hoc PCR amplifications (new Fig. S1M) revealed no bands indicative of concatemer formation. The primers used and the expected sizes of the amplicons are schematically represented in new Fig. S1M. In brief, we used a divergent primers set spanning the BSJ (3-4) to specifically detect the mature circRNA and a set of convergent primers (1-2) pairing on the GFP ORF, thus detecting both the circRNA and its linear precursor as well as the putative concatemer expected. Although a ~1 kb band was expected if a trans-splicing by-product was present, no such band was observed (new Fig. S1M). Moreover, RT-PCR amplification of the cGFP back-splice junction was markedly more efficient when reverse transcription was primed with random hexamers than with oligo(dT), priming total RNA or preferentially polyA+ RNA, respectively. These results are expected for a circRNA, as also indicated by the fact that the circZNF609 positive control behaves in a similar manner. Collectively, these results confirm the circular nature of our transcript and exclude translation originating from possible concatemers.

• These results are shown in new Fig. S1M and S1N and described in the text as follows: “Importantly, we ruled out the possibility that templates other than the GFP-encoding circRNA drive translation from our best performing constructs (V-cGFP, T3-glo-cGFP and T3-U3-cGFP). Ad hoc PCRs amplifications (Fig. S1M) revealed no bands indicative of concatemer formation. The left panel of Fig. S1M schematically illustrates the primer sets and expected amplicons sizes. In particular, we used a divergent primers set spanning the BSJ (3-4) to specifically detect the mature circRNA and a set of convergent primers (1-2) pairing on the GFP ORF detecting both the circRNA and its linear precursor as well as the putative concatemer expected. Although a ~1 kb band was expected if a trans-splicing by-product was present, no such band was observed. Moreover, RT-PCR amplification of the cGFP back-splice junction was markedly more efficient when reverse transcription was primed with random hexamers than with oligo(dT), priming total RNA or preferentially polyA+ RNA, respectively (Fig. S1M). These results are consistent with the circularity of the transcripts tested and coherent with the results obtained for circZNF609, used as control (Fig. S1M). Finally, cGFP resulted resistant to RNAseR treatment (Fig. S1N), further supporting its circular nature.”*

• There is a repetition of the world "a" in the abstract. We thank the reviewer for the attention paid to our text, we removed the extra “a” from the abstract.

• All circRNA translation studies should be cited when describing translation of circRNAs. We thank the reviewer for the suggestions, we corrected the mistake present in the text and included extra referenced about circRNA translation.

*Specifically, we included: *

• Fan, X., Yang, Y., Chen, C. et al. Pervasive translation of circular RNAs driven by short IRES-like elements. Nat Commun 13, 3751 (2022). https://doi.org/10.1038/s41467-022-31327-y
• Chen CK et al. Structured elements drive extensive circular RNA translation. Mol Cell. 2021 Oct 21; 81(20):4300-4318.e13.doi: 10.1016/j.molcel.2021.07.042. Epub 2021 Aug 25. PMID: 34437836; PMCID: PMC8567535.
• Obi P, Chen YG. The design and synthesis of circular RNAs. Methods. 2021 Dec;196:85-103. doi: 10.1016/j.ymeth.2021.02.020. Epub 2021 Mar 2. PMID: 33662562; PMCID: PMC8670866.
• Fukuchi, K., Nakashima, Y., Abe, N. et al. Internal cap-initiated translation for efficient protein production from circular mRNA. Nat Biotechnol (2025). https://doi.org/10.1038/s41587-025-02561-8
• Du, Y., Zuber, P.K., Xiao, H. et al. Efficient circular RNA synthesis for potent rolling circle translation. Nat. Biomed. Eng 9, 1062–1074 (2025). https://doi.org/10.1038/s41551-024-01306-3
• Wang F, Cai G, Wang Y, Zhuang Q, Cai Z, Li Y, Gao S, Li F, Zhang C, Zhao B, Liu X. Circular RNA-based neoantigen vaccine for hepatocellular carcinoma immunotherapy. MedComm (2020). 2024 Jul 29;5(8):e667. doi: 10.1002/mco2.667. PMID: 39081513; PMCID: PMC11286538.
• Andries O, Mc Cafferty S, De Smedt SC, Weiss R, Sanders NN, Kitada T. N(1)-methylpseudouridine-incorporated mRNA outperforms pseudouridine-incorporated mRNA by providing enhanced protein expression and reduced immunogenicity in mammalian cell lines and mice. J Control Release. 2015 Nov 10;217:337-44. doi: 10.1016/j.jconrel.2015.08.051. Epub 2015 Sep 3. PMID: 26342664.
• Yang Y, Fan X, Mao M, Song X, Wu P, Zhang Y, Jin Y, Yang Y, Chen LL, Wang Y, Wong CC, Xiao X, Wang Z. Extensive translation of circular RNAs driven by N6-methyladenosine. Cell Res. 2017 May;27(5):626-641. doi: 10.1038/cr.2017.31. Epub 2017 Mar 10. PMID: 28281539; PMCID: PMC5520850. REVIEWER 2

Circular RNAs (circRNAs) have attracted significant interest due to their unique properties, which make them promising tools for expressing exogenous proteins of therapeutic value. However, several limitations must be addressed before circRNAscan become a biologically and economically viable platform for the biotech industry.One of the main challenges is the reliance on large, highly structured sequences withinternal ribosome entry site (IRES) activity to initiate translation of the downstream open reading frame. In this study, the authors propose an alternative strategy that combines the 5′ untranslated region (5′UTR) of a previously characterized natural circRNA(circZNF609) with a short 13-nt nucleotide sequence shown to act as a translational enhancer. By evaluating the activity of various constructs containing a reporter geneacross multiple cell lines, they identify the most efficient and compact sequence, 63-nt long, capable of boosting translation within a circular RNA context.

Major Comments:

• This study is well-executed and relies on standard in vitro molecular biology techniques, which are adequate to support the conclusions drawn. *We thank the reviewer for the very positive opinion on the execution of our study. *

• The experimental procedures are clearly described, and the statistical analyses have been performed according to accepted standards. *We thank the reviewer for the very positive comment about the analyses we performed. *

Minor Comments:

• The manuscript would greatly benefit from a comprehensive revision to improve clarity and language. Involving a native English speaker during the editing process could significantly enhance the manuscript's readability and overall quality. The Results section would benefi t from closer attention, as certain parts of the description are attimes confusing and could be clarifi ed for better reader comprehension. We thank the reviewer for the input. We performed a huge revision of the text to improve language quality and enhance readability. We extended the descriptions in the results sections in order to explicit and clarify our data.

• The references should be carefully reviewed for accuracy and consistency-forinstance, references 9 and 10 appear to require correction or clarifi cation. We thank the reviewer for the careful reading of our paper. We amended the reference section, and we expanded it.

Reviewer #2 (Significance (Required)):

This study addresses a critical bottleneck in RNA therapeutics. The use of the proposed short sequences could significantly enhance the in vivo activity of protein-encoding circular RNAs. A highly efficient, compact translational enhancer has thepotential to substantially improve the therapeutic applicability of circRNAs and broaden their range of applications. Given the potential utility of these findings, we would anticipate pursuing intellectual property (IP) protection. To further strengthen the study, future work should include additional data on polysome association and a detailed analysis of the secondary structure of the 66-nt enhancer sequence. This work should be of broad interest to molecular biologists working on RNA biology, translation, and RNA-based therapeutics. I expect the identified sequence will betested by multiple laboratories to evaluate its strength and versatility, further underscoring the potential impact of this study. For context, I am actively engaged in research on non-coding RNAs.

• *

REVIEWER 3

In this brief report, the authors take advantage of circular RNA expression plasmids to define elements that can be used to enable efficient translation. They test a handful of known IRES elements as well as short translation enhancing elements (TEEs) for their ability to promote translation of circular GFP and c-ZNF609 reporters. They focus on one particular element that is of a short length and seems to work as well as longer IRES elements. My major concern relates to possible alternative sources of the translated proteins, which the authors have not ruled out (see below). I find themanuscript to be too preliminary in its current state.

• Work from the Meister group (Ho-Xuan et al 2020 Nucleic Acids Res 48:10368) has shown that apparent translation from circRNA over-expression plasmids is not from circular RNAs, but instead from trans-splicing linear by-products. The authors have not ruled out such alternative explanations here, e.g. by using deletion constructs that prevent backsplicing. We thank the reviewer for raising this point. *We ruled out the possibility that molecules other than the circRNA encoding GFP serve as templates for translation from our vectors. Specifically, ad hoc PCR amplifications performed for this purpose (new Fig. S1M) showed no bands that would indicate the presence of concatemers. Indeed, ad hoc PCR amplifications (new Fig. S1M) revealed no bands indicative of concatemer formation. The primers used and the expected sizes of the amplicons are schematically represented in new Fig. S1M. In particular, we used a divergent primers set spanning the BSJ (3-4) to specifically detect the mature circRNA and a set of convergent primers (1-2) pairing on the GFP ORF detecting both the circRNA and its linear precursor as well as the putative concatemer expected. Although a ~1 kb band was expected if a trans-splicing by-product was present, no such band was observed. Moreover, RT-PCR amplification of the cGFP back-splice junction was markedly more efficient when reverse transcription was primed with random hexamers than with oligo(dT), priming total RNA or preferentially polyA+ RNA, respectively. These results are expected for a circRNA, as also indicated by the fact that the circZNF609 positive control behaves in a similar manner. Collectively, these results confirmed the circular nature of our transcript and excluded translation originating from possible concatemers. *

These results are shown in new Fig. S1M and S1N and described in the text as follows: “Importantly, we ruled out the possibility that templates other than the GFP-encoding circRNA drive translation from our top constructs (V-cGFP, T3-glo-cGFP and T3-U3-cGFP). Ad hoc PCRs amplifications (Fig. S1M) revealed no bands indicative of concatemer formation. The left panel of Fig. S1M schematically illustrates the primer sets and expected amplicons sizes. In brief, we used a divergent primers set spanning the BSJ (3-4) to specifically detect the mature circRNA and a set of convergent primers (1-2) pairing on the GFP ORF detecting both the circRNA and its linear precursor as well as the putative concatemer expected. Although a ~1 kb band was expected if a trans-splicing by-product was present, no such band was observed (new Fig. S1M). Moreover, RT-PCR amplification of the cGFP back-splice junction was markedly more efficient when reverse transcription was primed with random hexamers than with oligo(dT), priming total RNA or preferentially polyA+ RNA, respectively (Fig. S1M). These results are consistent with the circularity of the transcripts tested (Fig. S1M). Importantly, cGFP PCR amplifications showed similar results as a validated endogenous circRNA, namely circZNF609, used as control (Fig. S1M, right panel), confirming the circular nature of cGFP. Finally, cGFP resulted resistant to RNAseR treatment (Fig. S1N), further supporting its circular nature.”* *

• Echoing the point above, the overall results would be stronger if the authors couldconfirm IRES activity using highly pure, in vitro transcribed RNAs that are transfected into cells * We thank the reviewer for this suggestion. Unfortunately, we are currently unable to produce synthetic circular molecules in-house, and the cost and time for purchasing synthetic ones are prohibitive. Nevertheless, we have performed the experiments described above to ensure the circularity of the transcripts tested.*

• The authors should also confirm their IRES activity using standard dual luciferase reporter (linear) constructs which have long been a standard approach in the field. We thank the reviewer for raising this point. As recommended, we cloned our three best candidates (VCIP, T3-glo, and T3-U3) into the pRL-TK/pGL3 dual-luciferase vector to assess their IRES activity (producing the vectors VCIP-Luc, T3-glo-Luc, and T3-U3-Luc), transfected them into RD cells, and, after 24 h of incubation, measured luciferase activity to assess the IRES performance of each candidate. From our analyses, VCIP and T3-U3 confirmed their IRES activity, although showing different relative efficiency, whereas T3-glo was inactive in the linear luciferase context. This finding is consistent with previous observations (Legnini et al., 2017) showing that the performance of IRES sequences in a linear luciferase reporter may differ from their activity when driving translation from a circRNA template. Overall, these results highlight the need for further investigation into the sequences and contexts specifically governing circRNA translation, rather than relying solely on knowledge derived from linear RNAs. *The results are shown below. We did not include them in the text to not overcomplicate the readability. However, we are happy to add and discuss them if required. *

***

***

Bar plot representing the relative luciferase activity deriving from VCIP-Luc (“V”), T3-glo-Luc (“T3-glo”), and T3-U3-Luc (“T3-U3”)*. Dual luciferase assay was performed and Renilla luciferase activity from each candidate was normalized against the Firefly luciferase. An empty ptKRL-pgl3 vector was used as reference. The ratio of each sample versus its experimental control was tested by two-tailed Student’s t test. * indicates a Student’s t test-derived p-value * *

• Methods, Plasmids Construction Section: Rather than including long lists of oligos and forcing a reader to figure out the final product that was cloned, it would be more intuitive if the authors provided the full sequences of the ORF and IRES sequencesthat were tested. We thank the reviewer for the comment, we added the sequences to the methods (Supplementary Table 1).

• The manuscript needs extensive English editing. Parts of it are also formatted in anunusual style, especially the introduction where it seems like each paragraph is a single sentence. As requested by the reviewer, we edited the text to make the language and content more accessible to readers.

• References included by the authors are selective and surprisingly do not include Chen et al (2021) Mol Cell 20:4300-4318 which already defined IRES elements for circRNAs that are fairly small. *Thank you for pointing this out. We have now cited the elegant work of Chen et al. (2021, Mol Cell 20:4300–4318) in the revised manuscript. While Chen and colleagues screened IRES-like elements of roughly 200 nt, our study was designed to uncover an even more minimal motif. The elements we report are therefore markedly shorter, highlighting a complementary, rather than overlapping, aspect of IRES available for driving circRNA translation. However, we now refer to Chen et al. in our text. *

• Error bars in Fig 2, especially Fig 2B, are huge. It seems impossible to make any conclusion given the large variety across these experiments. Thank you for your input. Although the error bars appear relatively large, the overall conclusions remain robust, as also noted by the other reviewers: both T3-glo and T3-U3 are intrinsically compact elements, yet they drive translation as efficiently as larger canonical IRESs. The error bars largely reflect the inherent variability of transient transfection assays, which naturally increases with the number of constructs examined. To strengthen our dataset without discarding existing replicates, we chose not to repeat experiments in the previously tested lines. Instead, we assessed our vectors in an additional model, the D283 medulloblastoma cell line. In this setting, we unexpectedly observed that the EMCV IRES surpasses the VCIP IRES, opposite to what we saw in the other lines, yet even here the short elements we identified remain strong competitors (new Fig. 2C, S2G, S2H). The evaluation of multiple CDSs across several cell lines, make our findings to be solid and well supported.

Hay gram+ y gram- igual que en bacterias?

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

Response to the Reviews

We thank the reviewers for their input and detailed feedback, which has helped us improve both the manuscript and the Microscopy Nodes software. Based on the comments, we have implemented new features, currently available as version 2.2.1 of Microscopy Nodes. We have edited the text and figures of the manuscript to reflect these changes and add clarification where needed.

Reviewer #1

Evidence, reproducibility and clarity

*The work by Gros et al. presents a paper introducing Microscopy Nodes, a new plugin for Blender 3D visualization software designed to import and visualize multi-dimensional (up to 5D) light and electron microscopy datasets. Given that Blender is not directly suited for such tasks, this plugin significantly simplifies the process, making its visualization engine accessible to a wide range of researchers without prior knowledge of Blender. The plugin supports importing volumes and labels from generic TIF or modern OME-Zarr image formats and includes supplementary video tutorials on YouTube to facilitate basic understanding of the visualization workflows.

Major comments: - The manuscript suggests that Microscopy Nodes can easily handle large datasets, as evidenced by the showcases. However, in my personal tests, I was unable to import a moderate TIF stack of about 5GB, which is considerably smaller than the showcased datasets. Post-import, a data cube was displayed, but the Blender interface became unresponsive. The manuscript should include a section stating limitations and addressing issues and providing suggestions for visualization of large datasets.*

We want to thank the reviewer for this valuable comment, which led us to find a core issue in Blender’s large data handling. Specifically, Blender’s rasterized pipeline causes issues with > 4 GiB of data loaded. This issue does not occur in the raytraced (Cycles) renderer, which is why we had not previously encountered it.

To address this, we have extended the reloading workflow of Microscopy Nodes to provide a workaround for this. If the data is larger than 4 Gibibytes (GiB) (per timepoint, or per timepoint per channel), Microscopy Nodes now automatically downsamples these data during import. While using these downsampled options is recommended for adjusting the visualization settings, the user can then still make their animation and reload their data to the largest scale for the final render by using the raytraced (Cycles) renderer. Additionally, we have raised this bug with the core Blender developers, and hope to work this out in the long term (blender/blender#136263).

We reflect these changes in the manuscript in the segment:

“Blender currently has a notable limitation that its default ‘quick’ rasterized rendering engines (such as ‘EEVEE’, but also the viewport ‘Surface’ and ‘Wireframe’ modes) do not support more than 4 Gibibytes (GiB) of volumetric data. The raytracing render mode ‘Cycles’, however, can handle large volumetric data. To allow users with large data to flexibly use Microscopy Nodes, we implemented a reloading scheme, where one first loads a smaller version of the data (under 4 GiB per timeframe for all loaded channels combined) - and only upon final render in Cycles, exchange it for the full/larger scale copy (Fig 3A). This downscaling of data offers additional benefits as it allows for fast adjustment of the render settings on e.g. a personal computer which can eventually be transferred to a larger workstation or HPC cluster for the final render at full resolution. This feature is critical as working in Cycles with larger files requires sufficient RAM to fit the (temporary) VDB files comfortably. For example, multiple figures in this manuscript were made on a 32GB RAM M1 Macbook Pro (Fig 1A, Video SV1, Fig 1D, Figure 2A-D, Fig S2A-B), but for larger data or long movies the movies were made on workstations or prepared on a laptop and then transferred to an HPC cluster for final rendering.”

* - The feature of importing Zarr-datasets over HTTP is great, but the import process was very slow in my tests, even on a robust network. For reference, loading 1.8 GB of the PRPE1_4x dataset at s1 level took 52 minutes. This raises concerns about potential code issues and general usability of the suggested workflow.*

We believe that this loading time may have been caused by the same issue that plagued all of our datasets of >4GB outside of the raytraced mode, as we have not seen loading issues like that. Moreover, Microscopy Nodes now supports Zarr version to Zarr 3/OME-Zarr 0.5, which allows ‘sharded’ Zarr datasets, which should be even faster at loading large blocks of data at the same time, as Microscopy Nodes does.

- The onsite documentation is a bit outdated and fails to fully describe the plugin settings.

We have updated our documentation to offer new written tutorials, which include full start-up tutorials, but also for some key extra instructions.

- The YouTube tutorials feature an outdated version of the plugin, which could confuse the general microscopy audience. These should be updated to better align with the current plugin functionality. Additionally, using smaller, easily accessible datasets for these tutorials would improve user testing experiences. Hosting complete (downsampled) demo project folder on platforms like zenodo.org could also enhance usability of such tutorials.

We have made a new series of YouTube tutorials that align with the current interface of Microscopy Nodes. These tutorials include public datasets, allowing users to follow along easily. We have chosen to also retain the older tutorials for users running legacy versions of the plugin, as they cover different workflows.

- The manuscript describes a novel dataset used in Fig. 2, but no reference is provided. Additionally, practical implementation of the coloring description for Fig. 2D can be unclear for inexperienced users, necessitating either step-by-step instructions or the provision of downsampled Blender files to aid understanding.

We have now shared the OME-Zarr address in the text (https://uk1s3.embassy.ebi.ac.uk/idr/share/microscopynodes/FIBSEM_dino_masks.zarr), and included this both in the manuscript and the tutorials. Additionally, to guide the implementation and explain the logic behind the coloring we introduced additional panels in Fig S1 and Fig S2 to showcase the shader setups used for this image.

[OPTIONAL] When importing labels, they can be assigned to individual materials only if initially split into multiple color channels. It would be great if the same logic is implemented when those materials are provided as indices within a single color channel. There can be a switch to define the logic used during the import process: e.g. the current one, when the objects are just colored based on a color map, or when they are arranged as individual materials as done when labels are imported from multiple color channels.

We agree with the reviewer and to address this concern with the update to version 2.2, we have implemented a new colorpicking system (See Fig 3B, inset 3, Fig 3C), this allows users to choose between a single color, various continuous, or categorical color maps.

Minor comments: - The manuscript shows nice visualizations of time series, light, and electron microscopy datasets, but in its current state, it is targeted more for light microscopy, where the signal is white. On the other hand, many EM datasets are rendered in inverted contrast (TEM-like), where the signal is black. To render such volume properly, it is needed to go into the Shading tab and flip the color ramp. Would it be possible to perhaps define the data type during import to accommodate various data types or perhaps select the flipped color ramp when the emission mode is switched off? It could make it easier for inexperienced EM users to use the plugin.

To address this, we include new default settings, with ‘invert colormaps on load’ option in the preferences, and default colors per channel (See Fig S4). We have also implemented a new color picking system in version 2.2 (See Fig 3B, inset 3, Fig 3C) that hopefully makes it easier before and after load to change colors.

- It was not completely clear to me whether it is possible to render a single/multiple EM slices using the inverted (TEM-like) contrast. For example, XY, XZ, YZ ortho slices across the volume. The manuscript contains: "This visualization is also supported in Blender, allowing for arbitrary selections of viewing angles (Fig 2B).", but it is not clear how to achieve that.

We introduced an additional explanation in Fig S1A and added a separate density window in the default shader to make this opaque view easier. To get a single slicing plane, users can reduce the scale of the slicing cube in one axis, at it is now also explained in Fig S2B.

- In 3D microscopy, it is quite common to have data with anisotropic voxels. As a result, the surfaces may require smoothing. I was not able to quickly find a way to smooth the surfaces (at least smooth modifiers for surfaces did not work for me). Is it possible to apply smoothing during the import of labels, or alternatively, smoothing of the generated surfaces can be a topic for an additional YouTube video.

The smoothness of the loaded masks can be indirectly affected in the preferences by changing the mesh resolution (changing the relative amount of vertices per pixel), but can be further affected by operations such as the Blender “Smooth” or e.g. the “Smooth by Laplacian” modifiers. To guide the users in doing so, we have included instructions for smoothing in the written tutorials on the website https://aafkegros.github.io/MicroscopyNodes/tutorials/surface_smoothing/ .

- It is also typical to have somewhat custom color maps for materials. It would be great if the plugin remembers the previously used color map for labels.

We have implemented new Preference settings, which include default colors and colormaps per channel, improving customization and reproducibility. This new option is described in Figure S4.

* - The pixel size edit box rounds up the values to 2 digits after the dot. Could it be changed to accommodate 3 or 4 digits as the units are um.*

Blender’s interface truncates the display, but stores higher-precision values internally, and become visible when users click or edit the values. We have added support for alternative pixel units to reduce the impact of the truncation.

- Import is not working when: - Start Blender - Select Data storage: with project - Overwrite files: on, set env: on, chunked: on - Select a file to import - Save Blender file - Pressing the Load button gives an error: "Empty data directory - please save the project first before using With Project saving."

We thank the reviewer for finding this bug which is now fixed in version 2.2.

- I was not able to play the downloaded supplementary video 3 using my VLC media player, while it was working fine in a browser. The video can be opened but looks distorted and heavily zoomed in. It may need to be re-saved from a video editor.

We have recompiled this video.

- References 12 and 16 are URL links instead of proper references to articles.

Thanks for catching this mistake in our bibliography. We have corrected this.

Significance

*This work effectively bridges a gap in the availability of tools for 3D microscopy dataset visualization. While many visualization programs exist, the high-quality ones are often expensive and thus not accessible to all researchers. The integration of Blender with Microscopy Nodes democratizes access to high-quality 3D visualization, enabling researchers to explore datasets and models from multiple perspectives, potentially leading to new discoveries and enhancing the understanding of key study findings. Despite its limitations, my experience with the plugin was engaging and useful. I would like to thank the authors for such useful work!

Limitations: - There remains a steep learning curve associated with using Microscopy Nodes, primarily due to Blender's complexity. More comprehensive tutorials could help mitigate this. - The conversion of imported images to Blender's internal 32-bit format results in a 4x increase in data size for 8-bit datasets. - Managing moderate-sized volumes (5-10 GB) can be challenging without clear strategies for effective handling. - The import of Zarr-datasets over the net is notably slow.

Audience: The plugin is suitable for a broad audience with a basic understanding of 3D visualization concepts, providing a solid foundation for exploring Blender's extensive features and options for optimal visualizations.

Reviewer expertise: Light microscopy, electron microscopy, image segmentation and analysis, software development, no experience with Blender*

Reviewer #2

*Evidence, reproducibility and clarity *

*Summary:

The article introduces Microscopy Nodes, a Blender add-on designed to simplify the loading and visualization of 3D microscopy data. It supports TIF and OME-Zarr images, handling datasets with up to five dimensions. The authors present different visualization modes, including volumetric rendering, isosurfaces, and label masks, demonstrating the application in light and electron microscopy. They provide examples using expansion microscopy, electron microscopy, and real-time imaging, highlighting how the tool enhances scientific communication and interactive visualization.

Comments:

However, some key aspects could be improved to enhance usability and reproducibility:

Example datasets: The images used in the YouTube tutorials were not accessible, making it difficult to reproduce the workflows shown in the figures and tutorials. It would be helpful if the authors provided direct links to the datasets or ensured that the same examples used in the tutorials were readily available for replication.*

We created new and updated tutorials and for all new tutorials, the data is now easily available from an S3 server.

Input file specifications: The article does not clearly detail how input files should be formatted. Many users will pre-visualize images in Fiji to convert their original images to a compatible format. It would be beneficial to specify which formats are supported for hyperstack creation, including details on bit depth, dimension ordering, label formats, and metadata compatibility, if applicable.

We have added new documentation on this on the website and in the manuscript. The addon can take 8, 16, and 32 bit data, and any dimension order (with the letters tzcyx) and pixel size. Dimension order and pixel size can be edited in the GUI. This is reflected in the manuscript in the rewritten section in Design and Implementation:

“It can handle 8bit to 32bit integer and floating point data, although all data types will be resaved into 32bit floating point VDB files, which can cause temporary files to take up more space than the original. Microscopy Nodes loads 2D to 5D files of containing data across time, z, y, x and channels, in arbitrary order (can be remapped in the user interface as well, Fig 3B, inset 2). To focus on relevant data, users can clip the time axis, which can be useful for long videos.”

* Hardware requirements: The article does not discuss RAM or hardware constraints in detail. In testing, attempting to load two images into the same project caused the program to freeze (tested on Mac M1). Specifying hardware requirements and limitations would help users manage expectations when working with large datasets.*

We have since found a limitation in the Blender engine that indeed limits the amount of data loaded (see also comment by Reviewer 1). Currently, rasterized engines are capped at 4 GiB, and only the raytraced engine can handle larger data. As such, the Microscopy Nodes pipeline, where one works with small images until it is time to render a final version, and the data is only exchanged for the final render, is still viable. To make this easier, we now also included optional downscaling for Tif images. This is described in the rewritten section on Design and Implementation:

“Blender currently has a notable limitation that its default ‘quick’ rasterized rendering engines (such as ‘EEVEE’, but also the viewport ‘Surface’ and ‘Wireframe’ modes) do not support more than 4 Gibibytes (GiB) of volumetric data. The raytracing render mode ‘Cycles’, however, can handle large volumetric data. To allow users with large data to flexibly use Microscopy Nodes, we implemented a reloading scheme, where one first loads a smaller version of the data (under 4 GiB per timeframe for all loaded channels combined) - and only upon final render in Cycles, exchange it for the full/larger scale copy (Fig 3A). This downscaling of data offers additional benefits as it allows for fast adjustment of the render settings on e.g. a personal computer which can eventually be transferred to a larger workstation or HPC cluster for the final render at full resolution. This feature is critical as working in Cycles with larger files requires sufficient RAM to fit the (temporary) VDB files comfortably. For example, multiple figures in this manuscript were made on a 32GB RAM M1 Macbook Pro (Fig 1A, Video SV1, Fig 1D, Figure 2A-D, Fig S2A-B), but for larger data or long movies the movies were made on workstations or prepared on a laptop and then transferred to an HPC cluster for final rendering.”

Significance

*General Assessment:

One of the major strengths of this work is its seamless compatibility with Blender, a powerful and widely used animation and 3D rendering tool. Integrating advanced visualization techniques from the animation and graphics industry into scientific imaging opens new possibilities for presenting complex microscopy data in an intuitive and accessible way. Additionally, the support for OME-Zarr is particularly valuable, as this format represents a major shift in bioimaging towards scalable, cloud-compatible, and standardized data storage solutions. The adoption of OME-Zarr facilitates large-scale data handling and improves interoperability across imaging platforms, making this integration a significant step forward for the field. Overall, the greatest strength of the tool lies in its flexibility for rendering microscopy data, but its accessibility for users without Blender experience might be a challenge.

Advance in the Field This work introduces a novel solution to the visualization challenges in microscopy by leveraging Blender's advanced rendering capabilities.

Audience This paper will be of interest to: Bioimage researchers seeking to enhance their microscopy data visualization. Image analysis tool developers interested in integrating advanced visualization into their workflows.

Field of Expertise This review is based on expertise in image analysis, segmentation, and 3D biological data visualization.*

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

The paper "Microscopy Nodes: Versatile 3D Microscopy Visualization with Blender" presents an easy and accessible approach for microscopists and microscopy users to visualize their data in a different and more controlled way. The authors have developed a plug-in script that enables the integration of complex 3D datasets into Blender, a widely used software for 3D visualization and illustration. By leveraging Blender's advanced rendering engine, the plug-in provides greater control over the scene, enviromint and presentation of the 3D data.

I believe that this development, especially when combined with additional analysis tools can be of a great value for microscopist and advanced users to presenting their 3D data sets.

However, at this stage, the paper does not seem to fully demonstrate the benefits of using Microscopy Nodes. To enhance the paper impact, it would be helpful for the authors to further emphasize and provide examples of how Blender's rendering specifically improves data presentation and, in turn, enhances the understanding of the data compared to existing solutions. Specifically, the authors claim at the end of the introduction that their development provides powerful tools for high-quality, visually compelling presentations, enabling "more effective communication of 3D biological data." I believe this statement should be supported by a figure comparing currently available visualization methods and demonstrating how using Blender enhances data presentation and by which enhances the communication of the results. *

*Additionally, at the end of the first paragraph of the results, the authors say: "These options allow us to combine the data and its analyzed interpretation in the same representation with Microscopy Nodes." However, this capability already exists in currently available software. Aside from now being able to achieve this in Blender, what additional benefits does it offer? *

We now include a new Table 1, to showcases which requirements for visualizing complex biological data are available in different visualization software, and discuss this in the text:

“Although several tools for 3D visualization of bioimages already exist and offer essential features for microscopy data (Table 1), many are proprietary, and open-source alternatives often struggle to deliver a comprehensive user experience, such as advanced animation and annotation controls. Proprietary solutions may offer some of these capabilities, but they are frequently limited by licensing costs, platform restrictions, and a lack of customizability. In contrast, Blender is a mature, well-supported open-source platform with a large community of developers that excels in both animation and visualization. By integrating microscopy-specific functionality through Microscopy Nodes, Blender becomes a uniquely powerful solution that bridges the gap between high-end graphics capabilities and the specialized needs of bioimage visualization.”

Additionally, we attempted to remake Figure 2C and 2D in the EM-field standard software Amira, but were not able to. This is because without an advanced light scattering algorithm, it is very hard to see the depth in the nucleus, and the semi-transparent masks do show each other behind them, but cannot interact with the volume.

We chose not to include this in the actual manuscript, as we are not experts at the Amira software, and will, by the nature of this manuscript, present a challenge that Blender is especially good at, such as here the combination of scattering light and semitransparent masks.

* In the last sentence of the second paragraph of the results, it is stated: "Blender powered by Microscopy Nodes: the ability to combine microscopy data with any 3D illustration in the same 3D environment." Could you please elaborate on the accuracy of the models that can be built and provide guidelines for achieving this using the data coordinates imported by Microscopy Nodes? If the illustrations are purely freehand and do not require specific accuracy, it would be helpful to clarify the advantages of creating them within the same environment rather than separately, as many scientists currently do. Additionally, if the inclusion of 3D model illustrations is one of the key advantages of using Blender, I believe it would be beneficial to present this in a figure rather than only in the supplementary video. *

We thank the reviewer for this comment and agree that in the previously submitted version of Microscopy Nodes, it was very difficult to align objects accurately, as the coordinate space was not transparent. A hurdle in this was the fact that Blender only works well with the unit ‘meters’. To address this issue, we now provide a choice of mapping the physical size to meters, as shown in the new interface (See Fig 3B, inset 5). Here the user can choose from the default ‘px -> cm’ (this will always look fine for a quick look) to options such as ‘nm -> m’ or ‘µm -> m’, which, combined with the new choice for adjusting the object origin upon load, allow users to treat the Blender coordinate space as based on the actual physical scales. Additionally, other Blender addons, such as Molecular Nodes (Reference 25 of the manuscript), also allow for accurate localization for cryo-EM datasets.

We appreciate the note that we should more clearly display the ability to show our illustrations and the data together in the figure and have added a visualization to show this in Figure 1C.

* Reviewer #3 (Significance (Required)):

The significance of the paper at this stage is primarily technical and mainly relevant to the field of microscopy

My field of expertise is microscopy and 3D visualization of models using mainly Maya3D and AMIRA.*

Test again

Author response:

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

Reviewer #1 (Public review):

Weaknesses:

(1) The manuscript's logical flow is challenging and hard to follow, and key arguments could be more clearly structured, particularly in transitions between mechanistic components.

We have revised our manuscript so as to make it easy for readers to follow the logical flow in transitions between mechanistic components by adding the descriptions of Figure S1E-J, Figure S2F-K, Figure S3A-H, Figure S4A-F, Figure S5, and Figure S6 in the revised manuscript.

(2) The causality between stress-induced α2A-AR internalization and the enhanced MAO-A remains unclear. Direct experimental evidence is needed to determine whether α2A-AR internalization itself or Ca2+ drives MAO-A activation, and how they activate MAO-A should be considered.

We believe that the causality between stress-induced α2A-AR internalization and the enhancement of MAO-A is clearly demonstrated by our current experiments, while our explanations may be improved by making them easier to understand especially for those who are not expert on electrophysiology.

Firstly, it is well established that autoinhibition in LC neurons is mediated by α2A-AR coupled-GIRK (Arima et al., 1998, J Physiol; Williams et al., 1985, Neuroscience). We found that spike frequency adaptation in LC neurons was also mediated by α2A-AR coupled GIRK-I (Figure 1A-I), and that α2A-AR coupled GIRK-I underwent [Ca²⁺]_i dependent rundown (Figures 2, S1, S2), leading to an abolishment of spike-frequency adaptation (Figures S4). [Ca²⁺]_i dependent rundown of α2A-AR coupled GIRK-I was prevented by barbadin (Figure 2G-J), which prevents the internalization of G-protein coupled receptor (GPCR) channels.

Abolishment of spike frequency adaptation itself, i.e., “increased spike activity” can increase [Ca²⁺]_i because [Ca²⁺]_i is entirely dependent on the spike activity as shown by [Ca²⁺]_i imaging method in Figure S3.

Thus, α2A-AR internalization can increase [Ca²⁺]_i through the abolishment of autoinhibition or spike frequency adaptation, and a [Ca²⁺]_i increase drives MAO-A activation as reported previously (Cao et al., 2007, BMC Neurosci). The mechanism how Ca²⁺ activates MAO-A is beyond the scope of the current study.

Our study just focused on the mechanism how chronic or sever stress can cause persistent overexcitation and how it results in LC degeneration.

(3) The connection between α2A-AR internalization and increased cytosolic NA levels lacks direct quantification, which is necessary to validate the proposed mechanism.

Direct quantification of the relationship between α2A-AR internalization and increased cytosolic NA levels may not be possible, and may not be necessarily needed to be demonstrated as explained below.

The internalization of α2A-AR can increase [Ca²⁺]_i through the abolishment of autoinhibition or spike frequency adaptation, and [Ca²⁺]_i increases can facilitate NA autocrine (Huang et al., 2007), similar to the transmitter release from nerve terminals (Kaeser & Regehr, 2014, Annu Rev Physiol).

Autocrine released NA must be re-uptaken by NAT (NA transporter), which is firmly established (Torres et al., 2003, Nat Rev Neurosci). Re-uptake of NA by NAT is the only source of intracellular NA, and NA re-uptake by NAT should be increased as the internalization of NA biding site (α2A-AR) progresses in association with [Ca²⁺]_i increases (see page 11, lines 334-336).

Thus, the connection between α2A-AR internalization and increased cytosolic NA levels is logically compelling, and the quantification of such connection may not be possible at present (see the response to the comment made by the Reviewer #1 as Recommendations for the authors (2) and beyond the scope of our current study.

(4) The chronic stress model needs further validation, including measurements of stress-induced physiological changes (e.g., corticosterone levels) to rule out systemic effects that may influence LC activity. Additional behavioral assays for spatial memory impairment should also be included, as a single behavioral test is insufficient to confirm memory dysfunction.

It is well established that restraint stress (RS) increases corticosterone levels depending on the period of RS (García-Iglesias et al., 2014, Neuropharmacology), although we are not reluctant to measure the corticosterone levels. In addition, there are numerous reports that showed the increased activity of LC neurons in response to various stresses (Valentino et al., 1983; Valentino and Foote, 1988; Valentino et al., 2001; McCall et al., 2015), as described in the text (page 4, lines 96-98). Measurement of cortisol levels may not be able to rule out systemic effects of CRS on the whole brain.

We had already done another behavioral test using elevated plus maze (EPM) test.By combining the two tests, it may be possible to more accurately evaluate the results of Y-maze test by differentiating the memory impairment from anxiety. However, the results obtained by these behavioral tests are just supplementary to our current aim to elucidate the cellular mechanisms for the accumulation of cytosolic free NA. Therefore, we have softened the implication of anxiety and memory impairment (page 13, lines 397-400 in the revised manuscript).

(5) Beyond b-arrestin binding, the role of alternative internalization pathways (e.g., phosphorylation, ubiquitination) in α2A-AR desensitization should be considered, as current evidence is insufficient to establish a purely Ca²⁺ -dependent mechanism.

We can hardly agree with this comment. 

It was clearly demonstrated that repeated application of NA itself did not cause desensitization of α2A-AR (Figure S1A-D), and that the blockade of b-arrestin binding by barbadin completely suppressed the Ca^2a-dependent downregulation of GIRK (Figure 2G-K). These observations can clearly rule out the possible involvement of phosphorylation or ubiquitination for the desensitization.

Not only the barbadin experiment, but also the immunohistochemistry and western blot method clearly demonstrated the decrease of α2A-AR expression on the cell membrane (Figure 3).

Ca²⁺-dependent mechanism of the rundown of GIRK was convincingly demonstrated by a set of different protocols of voltage-clamp study, in which Ca²⁺ influx was differentially increased. The rundown of GIRK-I was orderly potentiated or accelerated by increasing the number of positive command pulses each of which induces Ca²⁺ influx (compare Figure S1E-J, Figure S2A-E and Figure S2F-K along with Figure 2A-F). The presence or absence of Ca²⁺ currents and the amount of Ca²⁺ currents determined the trend of the rundown of GIRK-I (Figures 2, S1 and S2). Because the same voltage protocol hardly caused the rundown when it did not induce Ca²⁺ currents in the absence of TEA (Figure S1F; compare with Figure 2B), blockade of Ca²⁺ currents by nifedipine would not be so beneficial.

We believe the series of voltage-clamp protocols convincingly demonstrated the orderly involvement of [Ca²⁺]_i in accelerating the rundown of GIRK-I.

(6) NA leakage for free NA accumulation is also influenced by NAT or VMAT2. Please discuss the potential role of VMAT2 in NA accumulation within the LC in AD. 

It has been demonstrated that reduced VMAT2 levels increased susceptibility to neuronal damage: VMAT2 heterozygote mice displayed increased vulnerability to MPTP as evidenced by reductions in nigral dopamine cell counts (Takahashi et al, 1997, PNAS). Thus, when the activity of VMAT2 in LC neurons were impaired by chronic restraint stress, cytosolic NA levels in LC neurons would increase. We have added such discussion in the revised manuscript (page 12, lines 381-384).

(7) Since the LC is a small brain region, proper staining is required to differentiate it from surrounding areas. Please provide a detailed explanation of the methodology used to define LC regions and how LC neurons were selected among different cell types in brain slices for whole-cell recordings.

LC neurons were identified immunohistochemically and electrophysiologically as we previously reported (see Fig. 2 in Front. Cell. Neurosci. 16:841239. doi: 10.3389/fncel.2022.841239). We have added this explanation in the method section of the revised manuscript (page 15, lines 474-475). A delayed spiking pattern in response to depolarizing pulses (Figure S10 in the revised manuscript) applied at a hyperpolarized membrane potential was commonly observed in LC neurons in many studies (Masuko et al., 1986; van den Pol et al., 2002; Wagner-Altendorf et al., 2019).

Reviewer #2 (Public review):

Weaknesses:

(1) The manuscript reports that chronic stress for 5 days increases MAO-A levels in LC neurons, leading to the production of DOPEGAL, activation of AEP, and subsequent tau cleavage into the tau N368 fragment, ultimately contributing to neuronal damage. However, the authors used wild-type C57BL/6 mice, and previous literature has indicated that AEP-mediated tau cleavage in wild-type mice is minimal and generally insufficient to cause significant behavioral alterations. Please clarify and discuss this apparent discrepancy.

In our study, normalized relative value of AEP-mediated tau cleavage (Tau N368) was much higher in CRS mice than non-stress wild-type mice. It is not possible to compare AEP-mediated tau cleavage between our non-stress wild type mice and those observed in previous study (Zhang et al., 2014, Nat Med), because band intensity is largely dependent on the exposure time and its numerical value is the normalized relative value. In view of such differences, our apparent band expression might have been intensified to detect small changes.

(2) It is recommended that the authors include additional experiments to examine the effects of different durations and intensities of stress on MAO-A expression and AEP activity. This would strengthen the understanding of stress-induced biochemical changes and their thresholds.

GIRK rundown was almost saturated after 3-day RS and remained the same in 5-day RS mice (Fig. 4A-G), which is consistent with the downregulation of α2A-AR and GIRK1 expression by 3-day RS (Fig. 3C, F and G; Fig. 4J and K). However, we examined the protein levels of MAO-A, pro/active-AEP and Tau N368 only in 5-day RS mice without examining in 3-day RS mice. This is because we considered the possibility that a high [Ca²⁺]_i condition may have to be sustained for some period of time to induce changes in MAO-A, AEP and Tau N368, and therefore 3-day RS may be insufficient to induce such changes. We have added this in the revised manuscript (page 17, lines 521-525).

(3) Please clarify the rationale for the inconsistent stress durations used across Figures 3, 4, and 5. In some cases, a 3-day stress protocol is used, while in others, a 5-day protocol is applied. This discrepancy should be addressed to ensure clarity and experimental consistency.

Please see our response to the comment (2).

(4) The abbreviation "vMAT2" is incorrectly formatted. It should be "VMAT2," and the full name (vesicular monoamine transporter 2) should be provided at first mention.

Thank you for your suggestion. We have revised accordingly.

Reviewer #3 (Public review):

Weaknesses:

Nevertheless, the manuscript currently reads as a sequence of discrete experiments rather than a single causal chain. Below, I outline the key points that should be addressed to make the model convincing.

Please see the responses to the recommendation for the authors made by reviewer #3.

Reviewer #1 (Recommendations for the authors):

(1) Improve the clarity and organization of the manuscript, ensuring smoother transitions between concepts and mechanisms.

Please see the response to the comment raised by Reviewer #1 as Weakness

(2) Adjust any quantifying method for cytosolic NA levels under different conditions to support the link between receptor internalization and NA accumulation.

If fluorescent indicator of cytosolic free NA is available, it would be possible to measure changes in cytosolic NA levels. However, at present, there appeared to be no fluorescence probe to label cytosolic NA. For example, NS521 labels both dopamine and norepinephrine inside neurosecretory vesicles (Hettie & Glass et al., 2014, Chemistry), and BPS3 fluorescence sensor labels NA around cell membrane by anchoring on the cell membrane (Mao et al., 2023, Nat Comm). Furthermore, the method reported in “A Genetically Encoded Fluorescent Sensor for Rapid and Specific In Vivo Detection of Norepinephrine” is limited to detect NA only when α2AR is expressed. In the present study, increases in cytosolic NA levels are caused by internalization of α2AR. Cytosolic NA measurements with GRAB NE photometry may not be applicable in the present study. However, we have discussed the availability of such fluorescent methods to directly prove the increase in cytosolic NA as a limitation of our study (page 14, lines 429-436 in the revised manuscript).

(3) Include validation of the chronic stress model with physiological and behavioral measures (e.g., corticosterone levels and another behavioral test).

Please see the response to the comment raised by Reviewer #1 as Weakness (4).

(4) All supplemental figures should be explicitly explained in the Results section. Specifically, clarify and describe the details of Figure S1G-K, Figure S2F-K, Figure S3A-H, Figure S4A-F, Figure S5, and Figure S6 to ensure all supplementary data are fully integrated into the main text.

We have more explicitly and clearly described the details of Figure S1E-J, Figure S2F-K, Figure S3A-H, Figure S4A-F, Figure S5, and Figure S6 and fully integrated those explanations into the main text in the revised manuscript.

(5) In Figure 3, the morphology of TH-positive cells differs between panels D and E. Additionally, TH is typically expressed in the cytosol, but in the provided images, it appears to be localized only to the membrane. Please clarify this discrepancy and provide a lower-magnification image to display a larger area, not one cell.

In a confocal image, TH is not necessarily expressed homogenously in the cytosol, but is expressed in a ring-shaped pattern inside the plasma membrane, avoiding the cell nucleus and its surrounding Golgi apparatus and endoplasmic reticulum (ER) (Henrich et al., 2018, Acta Neuropathol Commun; see Fig. 4a and 6e), especially when the number of z-stack of confocal images is small. This is presumably because LC neurons are especially enriched with numerous Golgi apparatus and ER (Groves & Wilson, 1980, J Comp Neurol).

In Figure S7, we showed a lower-magnification image of LC and its adjacent area (mesencephalic trigeminal nucleus). In the LC area, there are a variety of LC neurons, which include oval shaped neurons (open arrowhead; similar to Figure 3D) and also rhombus-like shaped neurons (open double arrowheads, similar to Figure 3E). A much lower-magnification image of LC neurons constituting LC nucleus was shown in Figure 5A.

(6) In Figure 5, the difference in MAO-A expression is not clearly visible in the fluorescence images. Enzymatic assays for AEP and MAO-A should be included to demonstrate the increased activity better.

In the current study, we did not elaborate to detect the changes in TH, MAO-A and AEP in terms of immunohistochemical method. Instead, we elaborated to detect such changes in terms of western blot method. The main conclusions in the current study were drawn primarily by electrophysiological techniques as we have expended much effort on electrophysiological experiments. Because the relative quantification of active AEP and Tau N368 proteins by western blotting analysis may accurately reflect changes in those enzyme activities, enzymatic assay may not be necessarily required but is helpful to better demonstrate AEP and MAO-A activity. We have described the necessity of enzymatic assay to better demonstrate the AEP and MAO-A activities (page 10, lines 314-315).

Reviewer #3 (Recommendations for the authors):

(1) Causality across the pathway

Each step (α2A internalisation, GIRK rundown, Ca²⁺ rise, MAO-A/AEP upregulation) is demonstrated separately, but no experiment links them in a single preparation. Consider in vivo Ca²⁺ or GRAB NE photometry during restraint stress while probing α2A levels with i.p. clonidine injection or optogenetic over excitation coupled to biochemical readouts. Such integrated evidence would help to overcome the correlational nature of the manuscript to a more mechanistic study.

It is not possible to measure free cytosolic NA levels with GRAB NE photometry when α2A AR is internalized as described above (see the response to the comment made by reviewer #1 as the recommendation for the authors).

(2) Pharmacology and NE concentration

The use of 100 µM noradrenaline saturates α and β adrenergic receptors alike. Please provide ramp measurements of GIRK current in dose-response at 1-10 µM NE (blocked by atipamezole) to confirm that the rundown really reflects α2A activity rather than mixed receptor effects.

It is true that 100 µM noradrenaline activates both α and β adrenergic receptors alike. However, it was clearly showed that enhancement of GIRK-I by 100 µM noradrenaline was completely antagonized by 10 µM atipamezole and the Ca²⁺ dependent rundown of NA-induced GIRK-I was prevented by 10 µM atipamezole. Considering the Ki values of atipamezole for α2A AR (=1~3 nM) (Vacher et al., 2010, J Med Chem) and β AR (>10 µM) (Virtanen et al., 1989, Arch Int Pharmacodyn Ther), these results really reflect α2A AR activity but not β AR activity (Figure S5). Furthermore, because it is already well established that NA-induced GIRK-I was mediated by α2A AR activity in LC neurons (Arima et al., 1998, J Physiol; Williams et al., 1985, Neuroscience), it is not necessarily need to re-examine 1-10 µM NA on GIRK-I.

(3) Calcium dependence is not yet definitive

The rundown is induced with a TEA-enhanced pulse protocol. Blocking L-type channels with nifedipine (or using Cd²⁺) during this protocol should show whether Ca²⁺ entry is necessary. Without such a control, the Ca²⁺ link remains inferential.

The Ca²⁺ link was precisely demonstrated by a series of voltage clamp experiment, in which Ca²⁺ influx was orderly potentiated by increasing the number of positive voltage pulses (Figures S1 and S2). As the number of positive voltage pulses was increased, the rundown of GIRK-I was accelerated or enhanced more. The relationship between the number of spikes and the Ca²⁺ influx detected as Ca²⁺ transients was well documented in Ca2+ imaging experiments using fura-2 (Figure S3).

The presence or absence of Ca²⁺ currents and the amount of Ca²⁺ currents determined the trend of the rundown of GIRK-I (Figs. 2, S1 and S2). The same voltage protocol hardly caused the rundown when it did not induce Ca²⁺ currents in the absence of TEA (Fig. S1F; compare with Fig. 2B), and the series of voltage-clamp protocols convincingly demonstrated the orderly involvement of [Ca²⁺]_i in accelerating the rundown of GIRK-I. Therefore, blockade of Ca²⁺ currents by nifedipine may not be so beneficial.

(4) Age mismatch and disease claims

All electrophysiology and biochemical data come from juvenile (< P30) mice, yet the conclusions stress Alzheimer-related degeneration. Key endpoints need to be replicated in adult or aged mice, or the manuscript should soften its neurodegenerative scope.

As described in the section of Conclusion, we never stress Alzheimer-related degeneration, but might give such an impression. To avoid such a misunderstanding, we have added a description “However, the present mechanism must be proven to be valid in adult or old mice, to validate its involvement in the pathogenesis of AD.” (page 14, lines 448-450).

(5) Direct evidence for extracellular/cytosolic NE

The proposed rise in reuptake NA is inferred from electrophysiology. Modern fluorescent sensors (GRAB NE, nLight) or fast scan voltammetry could quantify NE overflow and clearance during stress, directly testing the model.

Please see the response to the comment made by Reviewer #1 as the Recommendations for the authors (2) as described above.

(6) Quantitative histology

Figure 5 presents attractive images but no numerical analysis. Please provide ROI-based fluorescence quantification (with n values) or move the images to the supplement and rely on the Western blots.

We have moved the immunohistochemical results in Fig. 5 to the supplement as we believe the quantification of immunohistochemical staining is not necessarily correct.

DOI: 10.1111/bph.70139

Resource: (Proteintech Cat# 22734-1-AP, RRID:AB_2879158)

Curator: @scibot

SciCrunch record: RRID:AB_2879158

What is this?

The essense of the proof is to note that the open neighbourhood family of \(y\) form a cauchy and open filter base, the inverse of it is also cauchy and open, thus y is in the image of f(B(0,1))

Suman 100%?

Sin grupo control!

Author response:

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

Reviewer #1 (Public review):

Summary:

Migration of the primordial germ cells (PGCs) in mice is asynchronous, such that leading and lagging populations of migrating PGCs emerge. Prior studies found that interactions between the cells the PGCs encounter along their migration routes regulates their proliferation. In this study, the authors used single cell RNAseq to investigate PGC heterogeneity and to characterize their niches during their migration along the AP axis. Unlike prior scRNAseq studies of mammalian PGCs, the authors conducted a time course covering 3 distinct stages of PGC migration (pre, mid, and post migration) and isolated PGCs from defined somite positions along the AP axis. In doing so, this allowed the authors to uncover differences in gene expression between leading and lagging PGCs and their niches and to investigate how their transcript profiles change over time. Among the pathways with the biggest differences were regulators of actin polymerization and epigenetic programming factors and Nodal response genes. In addition, the authors report changes in somatic niches, specifically greater non-canonical WNT in posterior PGCs compared to anterior PGCs. This relationship between the hindgut epithelium and migrating PGCs was also detected in reanalysis of a previously published dataset of human PGCs. Using whole mount immunofluorescence, the authors confirmed elevated Nodal signaling based on detection of the LEFTY antagonists and targets of Nodal during late stage PGC migration. Taken together, the authors have assembled a temporal and spatial atlas of mouse PGCs and their niches. This resource and the data herein provide support for the model that interactions of migrating mouse PGCs with their niches influences their proliferation, cytoskeletal regulation, epigenetic state and pluripotent state.

Overall, the findings provide new insights into heterogeneity among leading and lagging PGC populations and their niches along the AP axis, as well as comparisons between mouse and human migrating PGCs. The data are clearly presented, and the text is clear and well-written. This atlas resource will be valuable to reproductive and developmental biologists as a tool for generating hypotheses and for comparisons of PGCs across species.

Strengths:

(1) High quality atlas of individual PGCs prior to, during and post migration and their niches at defined positions along the AP axis.

(2) Comparisons to available datasets, including human embryos, provide insight into potentially conserved relationships among PGCs and the identified pathways and gene expression changes.

(3) Detailed picture of PGC heterogeneity.

(4) Valuable resource for the field.

(5) Some validation of Nodal results and further support for models in the literature based on less comprehensive expression analysis.

Weaknesses:

(1) No indication of which sex(es) were used for the mouse data and whether or not sex-related differences exist or can excluded at the stages examined. This should be clarified.

We have added: “Embryos of both sexes were pooled without genotyping, as the timepoints analyzed were prior to sex specification” to both the Animals section of the Materials and Methods and the Figure 1 legend. In addition, bioinformatic evaluation of potential sex biases in Nodal-Lefty signaling using Y-chromosome gene expression is reported in supplementary figure 4 and discussed in Discussion paragraph 2.

Reviewer #2 (Public review):

Summary:

This work addresses the question of how 'leading' and 'lagging' PGCs differ, molecularly, during their migration to the mouse genital ridges/gonads during fetal life (E9.5, E10.5, E11.5), and how this is regulated by different somatic environments encountered during the process of migration. E9.5 and E10.5 cells differed in expression of genes involved in canonical WNT signaling and focal adhesions. Differences in cell adhesion, actin cytoskeletal dynamics were identified between leading and lagging cells, at E9.5, before migration into the gonads. At E10.5, when some PGCs have reached the genital ridges, differences in Nodal signaling response genes and reprogramming factors were identified. This last point was verified by whole mount IF for proteins downstream of Nodal signaling, Lefty1/2. At E11.5, there was upregulation of genes associated with chromatin remodeling and oxidative phosphorylation. Some aspects of the findings were also found to be likely true in human development, established via analysis of a dataset previously published by others.

Strengths:

The work is strong in that a large number of PGCs were isolated and sequenced, along with associated somatic cells. The authors dealt with problem of very small number of migrating mouse PGCs by pooling cells from embryos (after ascertaining age matching using somite counting). 'Leading' and 'lagging' populations were separated by anterior and posterior embryo halves and the well-established Oct4-deltaPE-eGFP reporter mouse line was used.

Weaknesses:

The work seems to have been carefully done, but I do not feel the manuscript is very accessible, and I do not consider it well written. The novel findings are not easy to find. The addition of at least one figure to show the locations of putative signaling etc. would be welcome.

Thank you for the excellent suggestion. Fig. 6 has been added to highlight the main novel findings of this work and integrate them among contributions of earlier studies to provide a more complete view of signaling pathways and cell behaviors governing PGC migration.

(1) The initial discussion of CellRank analysis (under 'Transcriptomic shifts over developmental time...' heading) is somewhat confusing - e.g. If CellRank's 'pseudotime analysis' produces a result that seems surprising (some E9.5 cells remain in a terminal state with other E9.5 cells) and 'realtime analysis' produces something that makes more sense, is there any point including the pseudotime analysis (since you have cells from known timepoints)? Perhaps the 'batch effects' possible explanation (in Discussion) should be introduced here. Do we learn anything novel from this CellRank analysis? The 'genetic drivers' identified seem to be genes already known to be key to cell transitions during this period of development.

Thank you for this important observation. We have clarified the text in this section and added “This discrepancy may reflect differences in differentiation potential of some E9.5 PGCs that end in a terminal state among anterior E9.5 PGCs, but could also result from technical batch effects generated during library preparation. These possible interpretations are further discussed in the Discussion section.” to the pertinent results section and added additional relevant thoughts on the implications of this finding in Discussion paragraphs 4 and 7. We feel that it is important to include both results to the reader, as it is challenging to differentiate between heterogeneous developmental and migratory potential among E9.5 anterior PGCs and differential influence of batch effects across sequencing libraries with the data available.

(2) In Discussion - with respect to Y-chromosome correlation, it is not clear why this analysis would be done at E10.5, when E11.5 data is available (because some testis-specific effect might be more apparent at the later stage).

Since we had identified autocrine Nodal signaling primarily in anterior late migratory PGCs at E10.5 and knew that Nodal signaling was involved in sex specification of testicular germ cells into prospermatogonia by E12.5, we wanted to determine whether the Nodal signaling in late migratory PGCs at E10.5 was likely to be a sex-specific effect or was common to PGCs in both sexes. This was assessed in supplementary figure 4 and determined unlikely to be related to sex specification of PGCs as Nodal signaling was not strongly correlated with Y-chromosome transcripts in migratory PGCs. Assessing the relationship between Nodal signaling and Y-chromsome transcription at E11.5, when migration is complete, would be unlikely to help us further understand the dynamics of Nodal signaling during late PGC migration.

(3) Figure 2A - it seems surprising that there are two clusters of E9.5 anterior cells

Thank you for the interesting observation! One possibility is that the two states represent differential developmental competence as is suggested by the presence of one E9.5 anterior cluster along the differentiation trajectory in Fig 2A and one not within this differentiation trajectory. Another is that technical aspects of generating these sequencing libraries affected some cells more than others, resulting in clustering of highly affected and less affected cells, which would also be consistent with some E9.5 anterior cells lying within the differentiation trajectory and some not. Since it is challenging to differentiate between these possibilities with the data available, we have intentionally avoided overstating interpretations of this result in the manuscript text. We have included discussion of the potential implications of the transcriptional divergence you identify in Discussion paragraphs 4 and 7.

(4) Figure 5F - there does seem to be more LEFTY1/2 staining in the anterior region, but also more germ cells as highlighted by GFP

This is true; based on our selected anatomic landmarks for “anterior” and “posterior” as indicated in Methods, the “anterior” compartment typically contains more PGCs. Thus, we have included violin plots with all data points shown of signal intensities of both LEFTY1/2 and pSMAD2/3 in Fig. 5G and 5I so that the reader can evaluate the entire distribution of PGC signal intensities for each embryo.

Reviewer #3 (Public review):

Summary:

The migration of primordial germ cells (PGCs) to the developing gonad is a poorly understood, yet essential step in reproductive development. Here, the authors examine whether there are differences in leading and lagging migratory PGCs using single-cell RNA sequencing of mouse embryos. Cleverly, the authors dissected embryonic trunks along the anterior-to-posterior axis prior to scRNAseq in order to distinguish leading and lagging migratory PGCs. After batch corrections, their analyses revealed several known and novel differences in gene expression within and around leading and lagging PGCs, intercellular signaling networks, as well as number of genes upregulated upon gonad colonization. The authors then compared their datasets with publicly available human datasets to identify common biological themes. Altogether, this rigorous study reveals several differences between leading and lagging migratory PGCs, hints at signatures for different fates among the population of migratory PGCs, and provides new potential markers for post-migratory PGCs in both humans and mice. While many of the interesting hypotheses that arise from this work are not extensively tested, these data provide a rich platform for future investigations.

Strengths:

The authors have successfully navigated significant technical challenges to obtain a substantial number of mouse migratory primordial germ cells for robust transcriptomic analysis. Here the authors were able to collect quality data on ~13,000 PGCs and ~7,800 surrounding somatic cells, which is ten times more PGCs than previous studies.

The decision to physically separate leading and lagging primordial germ cells was clever and well-validated based on expected anterior-to-posterior transcriptional signatures.

Within the PGCs and surrounding tissues, the authors found many gene expression dynamics they would expect to see both along the PGC migratory path as well as across developmental time, increasing confidence in the new differentially expressed genes they found.

The comparison of their mouse-based migratory PGC datasets with existing human migratory PGC datasets is appreciated.

The quality control, ambient RNA contamination elimination, batch correction, cell identification and analysis of scRNAseq data were thorough and well-done such that the new hypotheses and markers found through this study are dependable.

The subsetting of cells in their trajectory analysis is appreciated, further strengthening their cell terminal state predictions.

Weaknesses:

Although it is useful to compare their mouse-based dataset with human datasets, the authors used two different analysis pipelines for each dataset. While this may have been due to the small number of cells in the human dataset as mentioned, it does make it difficult to compare them.

Direct comparisons between findings in human and mouse focused on CellChat cell-cell communication prediction results, which were conducted in an identical fashion using the same analysis methods for both datasets.

There were few validation experiments within this study. For one such experiment, whether there is a difference in pSMAD2/3 along the AP axis is unclear and not quantified as was nicely done for Lefty1/2.

Additional validation of the pSMAD2/3 signal intensity along the AP axis was performed and is now included in Fig. 5.

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

Learn more at Review Commons

Reply to the reviewers

Reviewer #1 (Evidence, reproducibility and clarity):

In Arabidopsis, DNA demethylation is catalyzed by a family of DNA glycosylases including DME, ROS1, DML2, and DML3. DME activity in the central cell leads to the hypomethylation of maternal alleles in endosperm. While ROS1, DML2, and DML3 function in vegetative tissues to prevent spreading DNA methylation from TE boundaries, their function in the endosperm was unclear.<br /> Using whole genome methylome analysis, the authors showed that ROS1 prevents hypermethylation of paternal alleles in the endosperm thus promotes epigenetic symmetry between maternal and paternal genomes.<br /> The approach and experimental desighs are appropriate, and the key conclusions are adequately supported by the results.<br /> However, there is not sufficient evidence to support the claim that DME demethylates the maternal allele at ROS1-dependent biallelically-demethylated regions. To clarify the issue, the authors could analyze if there is an overlap between DMRs identified in ros1 endosperm and those identified in dme endosperm using published data. If there is any, the authors could show a genome browser example of DMR including dme data.

Response: Thank you for your insight on our work. To address your concern and further test our model that DME prevents methylation of the maternal allele at regions where ROS1 is prevents methylation of the paternal allele, we turned to the allele-specific bisulfite-sequencing data published in Ibarra et al 2012. These data were from endosperm isolated at 7-8 DAP from aborting seeds of dme-2 +/- (Col-gl) plants pollinated by L_er_. Our analysis of these data is now included in Figures 6 and 7 and Supplemental Figures 13-17. We show that when the loss-of-function allele dme-2 is inherited maternally, average methylation of the maternal allele increases at ROS1-dependent regions (in the revised version of the paper now referred to as ROS1 paternal, DME maternal regions) from less than 10% CG methylation to approximately 40% CG methylation (Fig. 6D), consistent with our previous analysis using the non-allelic Hsieh et al 2009 data (now moved to Supplemental Figure 15). These results thus provide additional evidence that DME removes maternal allele methylation at regions where ROS1 removes paternal allele methylation (compare Fig. 6B and 6D). We included relevant genome browser examples in Figure 7E and Supplemental Figure 14. In the revised version, the relationship between ROS1 and DME is further expanded upon in the text.

Reviewer #1 (Significance):

Endosperm is a tissue unique to flowering plants. Though it is an ephemeral tissue, the endosperm plays essential roles for seed development and germination. The endosperm is also the site genomic imprinting occurs, and it has a distinct epigenomic landscape. This work provides a new insight that ROS1 may antagonize imprinted gene expression in the endosperm. However, it was not shown whether imprinted gene expression is indeed affected in ros1, or whether the ros1 mutation has phenotypic consequences. These results would be useful to discuss the evolution and significance of genomic imprinting.

Response: We agree that the biological significance of ROS1-mediated paternal allele demethylation is presently unknown. We performed RNA-seq on wild-type and ros1 3C and 6C endosperm nuclei, but these data were unfortunately not of high enough quality to include in the manuscript. In the Discussion we suggest that disrupting ROS1-mediated paternal allele demethylation might lead to a gain of imprinting over evolutionary time. In future work we are planning to address potential relationships to gene imprinting using a molecular, RNA-sequencing approach as well as an evolutionary comparative approach. As expected, given the expectation that imprinted genes are associated with a parent-of-origin specific epigenetic mark, we did not find any relationship between known imprinted genes and ROS1-dependent regions that are biallelically-demethylated regions in wild-type endosperm (see lines 362-372).

Reviewer #2 (Evidence, reproducibility and clarity):

SUMMARY

Hemenway and Gehring present evidence that the paternal genome in Arabidopsis endosperm is demethylated at several hundred loci by the DNA glycosylase/lyase ROS1. The evidence is primarily based on analysis of DNA methylation of ros1 mutants and of hybrid crosses where each parental genome can be differentiated by SNPs. I have some comments/questions/concerns, two of them potentially serious, but I think Hemenway and Gehring can address them through additional analyses of data that they already have available and a bit of clarification in writing.

Response: Thank you for your thoughtful review of this study. Your insight and suggestions have helped add clarity to the paper.

MAJOR COMMENTS:

1. Could the excess methylation in ros1-3 relative to ros1-7 shown in Figures 1A and 1C be explained by a second mutation in the ros1-3 background that elevates methylation at some loci? Any mutation that increased RdDM at these loci, for example could have this effect. This could confound the identification and interpretation of biallelicly demethylated loci.

Response: We propose a simpler explanation for the additional hypermethylation observed in ros1-3: ros1-3 is a loss-of-function (null) allele whereas ros1-7 is likely a hypomorphic allele. For clarity, we have added a diagram of all of the alleles used in this study as Supplemental Figure 1B. The ros1-3 allele was first described in Penterman et al, PNAS, 2007. It is a T-DNA insertion allele that was isolated in the Ws accession and then backcrossed 6 times to Col-0, greatly minimizing the risk of unlinked secondary mutations being present. There is no genetic evidence that there is another T-DNA insertion in this line. The ros1-7 allele was described in Williams et al, Plos Genet, 2015. It was isolated from the Arabidopsis Col-0 TILLING population and is missense mutation (E956K) in a residue in the glycosylase domain that is conserved among the four DNA glycosylases. It is known that ROS1 transcripts are produced from the ros1-7 allele (Williams et al 2015). We observe less hypermethylation in the ros1-7 background compared to the ros1-3 background, and thus propose that the ros1-7 allele is a hypomorphic allele of ROS1. The use of two independent ros1 mutant alleles for initial endosperm methylation profiling strengthens the findings of our study. Importantly, regions that are hypermethylated in ros1-3 are also hypermethylated in ros1-7, but to a lesser extent, and vice versa (Fig 1D, Supplemental Figs. 3 and 4).

We also use a third allele in this study, ros1-1, which is a nonsense allele in the C24 accession. Notably, we find that the regions are demethylated on both maternal and paternal alleles in wild-type C24 gain DNA methylation primarily on the paternal allele in ros1-1 endosperm (Figure 4C,D and Supplemental Figure 10). This is discussed further in response to your second point.

Given these lines of evidence, a gain-of-function mutation in a methylation pathway, like RdDM, in the ros1-3 background is an unlikely explanation for increased hypermethylation compared to ros1-7. The use of three independent ros1 alleles for methylation profiling, all of which lead to the same conclusions, is a major strength of our study.

1. It appears that the main focus of the manuscript, the existence of loci that are paternally demethylated by ROS1, is supported by a set of 274 DMRs. This is a small number relative to the size of the genome and raises suspicions of rare false positives. Even the most stringent p-values that DMR-finding tools report do not guarantee that the DMRs are actually reproducible in an independent experiment. Demonstrating overlap between these 274 DMRs and an independently defined set using a different WT control and different ros1 allele would suffice to remove this concern. It appears that authors already have the needed raw data with ros1-1 and ros1-7 alleles.

Response: First, we should clarify that paternal demethylation by ROS1 is supported by more than the 274 DMRs. All ros1 CG hyperDMRs show an increase in paternal allele methylation in ros1 (Fig. 4B,D). The 274 DMRs are a distinct subset defined as having less methylation on the maternal allele than the paternal allele in ros1 endosperm and where there is no maternal allele hypomethylation in wild-type endosperm (refer to Fig. 5B).

We agree with your sentiments about DMR-finders and we are cautious of relying exclusively on DMR calls when making conclusions. We verify the nature of identified DMRs using metaplots and weighted average comparisons throughout the paper, which we think increases confidence in the conclusions and goes beyond a simple DMR-calling approach.

We argue that we have replicated the major conclusion of the paper, that ROS1 prevents paternal allele hypermethylation at target regions in the endosperm, in the following ways:

1. In the dataset without allelic-specific methylation information (Figures 1-3), we found that both ros1-3 and ros1-7 CG hyperDMRs have a limited capacity for hypermethylation in the endosperm relative to leaf or sperm (Table 1, Fig 3, Supplemental Fig. 4). In the allele-specific dataset, ros1-3 CG hyperDMRs were revealed to have particularly low maternal mCG relative to paternal mCG in ros1 mutant endosperm (Fig 4A-B, Supplemental Fig. 10).
2. We found that ros1-3 and ros1-1 hyperDMRs, which we identified using non-allelic data, are biased for paternal allele hypermethylation in the endosperm of F1 hybrids (Fig 4B,D). The replicability of the paternal bias in hypermethylation in both ros1-3 in the Col-0 ecotype and ros1-1 in the C24 ecotype is a critical result, and we have moved the ros1-1 hyperDMR plots from the supplement to main figure 4C-D in the revised version of the manuscript as a result of your comment.
3. The 274 DMRs identified as “biallelically-demethylated, ROS1-dependent” are by definition replicated between reciprocal cross directions. (Note that we now refer to these regions as ROS1 paternal, DME maternal regions in the revision.) Regions in this category had to be called as maternally-hypomethylated in both ros1-1 x ros1-3 and ros1-3 x ros1-1 endosperm. These regions also had to not be identified as maternally-hypomethylated in both C24 x Col-0 and Col-0 x C24. We hope this is clarified for readers by Table 1, which we have included based on your suggestion in comment #3, as well as other clarifying edits we made in this section of the paper.comparisons between maternal and paternal methylation in endosperm, DMRs defined by comparison between mutants and wildtype, and more. These need clearer descriptions of which sets are being referred to throughout the main text and in figure legends. A table summarizing them might help (not in the supplement). Use of consistent and precisely defined terms would help. Stating the number of DMRs along with the name for each set would help a lot, even though this would make for some redundancy. (The number of DMRs in each set not only helps with interpretation but also act as a sort of ID). The reason I put this as a major concern is because the text and figures are difficult to understand, and it is currently hard to evaluate both the results and the authors' conclusions from those results.

Response: Thank you for your feedback and suggestions. We have edited the main text so that only one descriptive name is used for each DMR type throughout the paper. We have also renamed regions for greater clarity. The previous “ROS1-independent, maternally demethylated regions” are now referred to as “DME maternal regions”. The previous “ROS1-_independent, biallelically-demethylated regions” are now referred to as “_ROS1 paternal, DME maternal regions”. These changes provide greater clarity and also emphasize the role of DME at regions that are paternally hypermethylated in ros1. We have added Table 1 to summarize the DMR classes of interest.

MINOR COMMENTS

1. The sRNA results in Figure 2B are difficult to interpret because they do not reveal anything about the number of TEs that have siRNAs overlapping them or their flanks. While the magnitude of some of the highest endosperm sRNA peaks is higher than the embryo peaks, that could be explained by a small number of TEs with large numbers of sRNAs. To make this result more interpretable, we also need some information about how many TEs have a significant number of sRNAs associated with them in endosperm and embryo in each region (e.g., middle, 5', 3', and flanks of TEs). What a "significant number of sRNAs" is would be up to the authors to decide based on the distribution of sRNA counts they observe for TEs. Perhaps the top quartile of TEs? Combined with the same analysis done in parallel with non-ROS1 target TEs, this would reveal whether there is any evidence for ROS1 counteracting sRNA-driven methylation spread from TEs.

Response: Thank you for the suggestion. We now present these data and the data for individual TEs underlying the metaplots in Supplemental Figure 7. As suggested by the reviewer, ROS1 TEs do not have uniformly higher levels of sRNA in their flanks in the endosperm compared to the embryo. We have modified our interpretations accordingly.

1. The statement "we are likely underestimating the true degree of differential methylation among genotypes" should be validated and partially quantified using a methylation metaplot like Figure 2A, but substitute DMRs for TEs. Related to that, Figure 1B needs an indicator of scale in bp.

Response: We have now included a methylation metaplot over ros1-3 hyperDMRs and ros1-7 hyperDMRs as Supplemental Figure 3 These plots show that indeed there is additional hypermethylation in DMR-proximal regions. We have added a scale bar to Figure 1B and other browser examples in the paper.

1. The statement "Over half of ROS1 target regions identified in the ros1-3 mutant endosperm were within 1 kb or intersecting a TE (Fig. 1D)" is hard to interpret without some kind of ROS1 non-target regions or whole-genome control comparison. How different are the numbers in Fig. 1D from a random expectation?

Response: We have now included a control for random regions in Figure 1E. We define these as regions where there was sufficient methylation data coverage and a low enough methylation level in wild-type to detect hypermethylation if it existed.

1. The sentence at line 262 is confusing. Is the comparison between dme mutant and ros1 mutant or between different types of regions? And it appears that the comparison value is missing in the "3-5% CG methylation gain..." e.g., "3-5% CG methylation vs 10-20%" or something like that.

Response: This section has been re-written as we now focus on allele-specific dme endosperm methylation data for our comparisons.

1. The dme mutant data in Figure 5C appear to be key to the model in Figure 7. The relative impact of the dme mutant in the two types of regions should be quantified.

Response: Thank you for this comment. To further probe our model that DME prevents hypermethylation of the maternal allele at regions where ROS1 is preventing hypermethylation of the paternal allele, we turned to the allele-specific bisulfite-sequencing data published in Ibarra et al 2012 (see also response to reviewer #1). Using these data, we show that when the loss-of-function allele dme-2 is inherited maternally, ROS1 paternal, DME maternal regions (previous referred to as ROS1-_dependent, biallelically-demethylated regions) are CG hypermethylated on the maternal allele (Figure 6D). Thus, these results both replicate the observations made with the Hsieh et al 2009 data, and provide additional evidence that _DME prevents maternal allele hypermethylation at regions were ROS1 is preventing paternal allele hypermethylation. These results have replaced the Hsieh et al 2009 results in Figure 6, and we have moved the analysis of Hsieh et al 2009 data to Supplemental Figure 15.

1. Looks like sRNA methods are missing.

Response: Thank you for identifying this. We previously included the reference for the analyzed dataset we used and the method for plotting under an unclear section header. These methods are now in the section “Analysis of average methylation and 24-nt sRNA patterns for features of interest”, and we have added additional reference to the specific dataset we used.

1. Supplemental Figure 1 is hard to interpret since it only list gene IDs, not gene names.

Response: As suggested, we have added gene names to this figure.

The last comments are suggestions for increasing the impact of this study:

1. Figure 2A and 3B suggest that ROS1 target TEs show demethylation in their flanks but not in the TE themselves. This is an interesting result. If it is true, more DMRs would be expected in the ROS1 target flanks than in the ROS1 target TEs. Reporting how many ROS1 target TEs have DMRs in them and what proportion have DMRs in their flanking 1-Kb regions would answer this question. Given the significance of this result, it also deserves a bit more context: Is the magnitude of increased methylation flanking TEs in ros1 mutant endosperm different than in ros1 mutant leaves or other tissue? Does methylation in TE flanks behave the way in dme mutant endosperm?

Response: We define “ROS1 target TEs” (now referred to more simply as ROS1 TEs) as TEs within 1kb or intersecting a ros1-3 hyperDMR. Consistent with your interpretation, 80% of the TEs in this category do not have a DMR overlapping them, instead they have a TE within 1kb. We now mention this in the text on line 150.

The total level of DNA methylation at ROS1 TEs is lower in the endosperm than in leaf, as DNA methylation levels are overall lower in endosperm than in leaf. The magnitude of increased methylation flanking TEs in ros1 mutant endosperm is not different between the two tissues. This is observable in Supplemental Fig. 5 in the revised version of the paper, and we report this result in the revised text. In the revision we also present methylation profiles of DME TEs in WT and ros1 endosperm (Fig. 7B-D). DME TEs are hypomethylated in both the body and flanks in WT and ros1.

1. The idea of biallelic demethylation has been theoretically suggested in maize to explain weak overlap between endosperm DMRs and imprinting (Gent et al 2022). If that were true in Arabidopsis, then ROS1 target, biallelicly demethylated loci would be less likely to have imprinted expression than maternally demethylated loci. This prediction could be tested using available data in Arabidopsis.

Response: Indeed, as you hypothesize, there are no known imprinted genes (Pignatta et al 2014) associated with biallelically-demethylated, ROS1-dependent regions (now referred to as ROS1 paternal, DME maternal regions). Expectedly, there are imprinted genes associated with maternally-demethylated regions (now referred to as DME regions). 23 imprinted genes identified in the Pignatta et al 2014 study are within 1 kb or intersecting a DME region. This is discussed on lines 364-374.

1. There is currently no evidence for biological significance of biallelicly demethylated loci. Knowing where they are in the genome might give some hints. A figure like Fig. 1D but specifically showing the biallelicly demethylated DMRs would be valuable.

Response: This is now included in Figure 7A.

1. It is hard to make the comparisons between genotypes and parental genomes in Figure 6 and know what they mean. Maybe a different way of displaying the data would help. Or maybe even a different labeling system could make it a little more accessible.

Response: We have revised this figure (now Fig. 8) in the following ways, which we believe address your comments and clarify the main conclusions:

Figure 8C is now a boxplot comparing methylation of the paternal allele of ROS1 paternal, DME maternal regions (previously referred to as biallelically-demethylated, ROS1-dependent regions) across endosperm ROS1 genotypes. This plot shows increased methylation of paternal alleles when the paternal parent is a ros1 mutant, regardless of whether the resultant F1 endosperm is homozygous or heterozygous for ros1 (columns 3, 4, 6).

Figure 8B remains as a scatterplot, where we can observe significant correlation between individual ROS1 paternal, DME maternal regions in homozygous ros1 endosperm and heterozygous ros1/+ endosperm. Note that paternal allele methylation is higher in homozygous ros1 endosperm for most regions.

Reviewer #2 (Significance):

Demethylation of the maternal genome in endosperm has been the subject of much research because it can result in genomic imprinting of gene expression. The enzymes responsible, DNA glycosylases/lyases, also demethylate DNA in other cell types as well, where DNA methylation is not confined to one parental genome (biallelic or biparental as opposed to uniparental demethylation). To the best of my knowledge, the extent or even existence of biallelelic demethylation in endosperm has not been studied until now (except for a superficial look in a bioRxiv preprint, https://www.biorxiv.org/content/10.1101/2024.07.31.606038v1). Hemenway and Gehring have carried out a thoughtful and detailed analysis of the topic in Arabidopsis at least as far as it depends on the DNA glycosylase ROS1.

A limitation is that the study design would miss biallelic demethylation by any of the other three DNA glycosylases in Arabidopsis. A second limitation is that there is no clear biological significance, just some conjecture about evolution. Nonetheless, given the novelty of the topic, biological significance may follow.

The audience for biallelic DNA demethylation in Arabidopsis endosperm is certainly in the "specialized" category, but its relevance to the larger topic of gene regulation in endosperm will attract a larger audience.

Response: With regard to the other demethylases, note that we also profiled methylation in ros1 dml2 dml3 triple mutant endosperm. We did not find evidence for many DMRs that were present in the triple mutant that were not present in the ros1 single mutant. We do not rule out a function for DML2 or DML3 in the endosperm, but this is not observed at the level of bulk endosperm.

The reviewer is correct that we have shown a molecular phenotype (paternal allele hypermethylation) and not a developmental or morphological phenotype. A function that occurs in one parent but not the other is, to us, exciting. Our thoughts about how this finding might relate to imprinting are indeed speculative, but not wildly so.

Reviewer #3 (Evidence, reproducibility and clarity):

DNA demethylases play a key role in DNA methylation patterning during flowering plant reproduction. The demethylase DME, in particular, is critical for proper endosperm development. While the function of DME in endosperm development has been explored, the contributions of the other demethylases in the same family, ROS1, DML2 and DML3 in Arabidopsis, have not yet been investigated. In vegetative tissues, ROS1 prevents hypermethylation of some loci. In this work, Hemenway and Gehring explore whether ROS1, DML2 and DML3 also affect DNA methylation patterns in endosperm. Using EM-seq of sorted endosperm nuclei, they show that loss of ROS1 indeed causes hypermethylation of a number of loci, particularly the flanks of methylated transposons, while loss of DML2 and DML3 has minimal additional effect. By obtaining allele-specific EM-seq data through crosses of Col and C24, the authors show that ros1 endosperm hypermethylation is mostly restricted to the paternal allele. The authors propose that at some sites, ROS1 helps bring down paternal methylation levels to match maternal methylation levels, which are typically reduced in endosperm due to DME activity in the female gametophyte prior to fertilization. In a ros1 mutant with paternal hypermethylation, these sites become differentially methylated on the maternal and paternal alleles, resembling imprinted loci. This work convincingly establishes a function for ROS1 in DNA methylation patterning in endosperm. However, I struggled with the clarity of the writing and reasoning in a few places, and would suggest clarification of a few points and additional analyses below.

Response: Thank you for your thoughtful review of our paper. Your questions and suggestions have been invaluable in revising the work.

I think making a few simple changes to streamline nomenclature would improve readability. For example, in the section starting on line 129, the same set of genomic features are called ROS1 target-proximal TEs, TEs that are near a ROS1 target region, and ROS1 target-associated TE regions. Also for example in line 254 "regions that are maternally-demethylated in wild-type endosperm, and are not dependent on ROS1 for proper demethylation" - are these the same as the "ROS1-independent, maternally-demethylated" regions in Fig. 5a? Given how complex these terms are, being consistent throughout the manuscript really helps the reader.

Response: We edited the text and figures so that only one descriptive name is used for each DMR class or region throughout the paper. Thank you for this feedback; these edits have made the paper much clearer.

Is there any notable effect of ros1 on gene expression in endosperm? Endosperm is a terminal tissue, so maintaining DNA methylation boundaries as ROS1 does in vegetative tissues seems less important. It begs the question of why ROS1 is doing this in endosperm, is it just because it's there, or is there an endosperm-specific function? Exploring effects on imprinting would be particularly interesting (does loss of ROS1 'create' imprinted loci at these newly asymmetrically methylated sites?) but probably beyond the scope of the present work.

Response: We agree, the question of the functional consequence of ROS1 activity in the endosperm is something we are keen to address in future work. We performed RNA-seq on wild-type and ros1 3C and 6C endosperm nuclei, but these data were unfortunately not of high enough quality to include in the manuscript. We are in particular interested in this question you have proposed – if loss of ROS1 can ‘create’ imprinted loci. We are planning to address this both using a molecular, RNA-sequencing approach as well as an evolutionary comparative approach. This is an important and exciting future direction.

Is DME expressed in sperm, or is expression of DME affected in ros1 sperm or endosperm? One other explanation for ros1 hypermethylation occurring primarily on the paternal allele is that, potentially, DME can substitute for ROS1 in the central cell where DME is already very active, but not in sperm cells. Related, how well expressed is ROS1 vs. DME in sperm cells?

Response: This is an important series of questions, and something we are very interested in as well. Studies of Arabidopsis pollen have shown that both ROS1 and DME, while they prevent some hypermethylation in sperm, are more active in the vegetative nucleus of pollen than in sperm. ROS1 is expressed at a low level in the microspore and bicellular pollen and DME is expressed at a low level throughout pollen development. We have included Supplemental Fig. 17 with available expression data to make this point in the paper. Likely, any effects of loss of ROS1 or DME on sperm DNA methylation are inherited from precursor cells (Ibarra et al 2012, Calarco et al 2012, Khouider et al 2021). Your proposal that perhaps DME can sub in for ROS1 in the central cell but not in sperm is intriguing. Unfortunately there’s not enough data in the central cell to convincingly address this at this time.

To investigate the relationship between DME and ROS1 in the male germline, we used the bisulfite-sequencing data generated in sperm cells in Khouider et al 2021. We calculated average DNA methylation levels in dme/+, ros1, dme/+;ros1, and wild-type Col-0 sperm cells at ROS1 paternal, DME maternal regions, shown in Supplemental Fig. 18A. We observed little increase in mCG methylation in dme/+ sperm relative to wild-type Col-0 sperm. This is consistent with your proposed model that DME is unable to demethylate these regions outside of the female germline. As expected, there is increased mCG in ROS1 paternal, DME maternal regions in ros1-3 mutant sperm relative to wild-type Col-0 sperm. DME maternal regions are highly methylated in wild-type Col-0 sperm.

Fig 2b shows that ROS1 target-associated TEs are enriched for sRNAs in endosperm relative to embryo, whereas the reverse is true for non-ROS1-assoc TEs. Since TEs are not always well annotated and some may be missing from this analysis, what about trying the reverse analysis - are regions enriched for 24nt sRNAs in endosperm significantly hypermethylated in ros1 endosperm? All regions or only some?

Response: We performed an analysis to address your inquiry and observed a low magnitude increase in DNA methylation in ros1 mutant endosperm at regions defined by Erdmann et al as more sRNA producing in the endosperm relative to the embryo (endosperm DSRs). Endosperm DSRs are generally lowly methylated in wild-type endosperm, as was observed originally in Erdmann et al 2017. Small increases in DNA methylation are observed at endosperm DSRs in all sequence contexts in ros1 endosperm. Overall, this is consistent with ROS1 targets being a subset of sRNA-producing regions in the endosperm. This analysis is now included in Supplemental Fig. 7C.

What is the relationship between previously-defined DME targets and ROS1 targets identified in this paper? DME tends to target small euchromatic TE bodies, whereas Fig. 3 suggests that ROS1 helps prevent methylation spreading on the outer edges of the TEs, rather than in the TE body. Do all DME targets tend to be adjacent to or flanked by ROS1 target sites? Or are the TEs affected by DME (in body) and by ROS1 (at edges) largely nonoverlapping? Fig. 5a suggests that the ROS1-dependent, biallelically-demethylated sites are both DME and ROS1 targets, but how often do these really appear to overlap? More than by chance?

Response: We have sought to address your comments through a series of analyses that we have included in Fig. 7 and Supplemental Fig. 16. We found that ROS1 paternal, DME maternal regions (formerly referred to as ROS1-dependent, biallelically-demethylated regions) and DME maternal regions (formerly referred to as ROS1-independent, maternally-demethylated regions) do not occupy the same genomic regions. However, we do observe some evidence for ROS1 activity in flanking regions of DME targets (Fig. 6A, Fig. 7B-D). To look at TEs specifically, as you suggest, we first identified TEs that were within 1kb or intersecting a DME maternal region. Based on our characterization of these regions, we assume these to be DME-targeted TEs. We then performed ends analysis to see if there was evidence of ROS1 activity at the ends of these TEs. Indeed, at a global level there is a slight hypermethylation of the paternal allele in a ros1 mutant at the end of these DME TEs (Fig. 7B). To better visualize how many DME TEs are showing ROS1 activity at their ends, we then plotted the difference between the median ros1-3 methylation and median Col-0 values in the non-allelic endosperm for each TE in a clustered heatmap (Fig. 7C). The parent-of-origin data does not have enough coverage for clustering in this way, so we used the non-allelic data. A small fraction of “DME TEs” gain methylation in the ros1 mutant endosperm relative to wild-type (Fig. 7C-D).

Are the TEs whose boundaries are demethylated by ROS1 more likely to be expressed in vegetative or endosperm tissues than TEs not affected by loss of ROS1? Expressed TEs likely produce more sRNAs, which would increase RdDM in a way that might need to be more actively countered by ROS1 than transcriptionally silent or evolutionarily older TEs.

Response: This is an interesting line of inquiry, although perhaps out of the scope of our present study. It has been shown that TEs demethylated by ROS1 are targeted by the RdDM pathway in Arabidopsis vegetative tissue (Tang et al 2016). Using data from Erdmann et al 2017, we looked at 24 nt sRNAs at ROS1-TEs in the endosperm and embryo (Supplemental Fig. 7). sRNA production at ROS1 TE-flanking regions is observed in both embryo and endosperm, but clearly not all ROS1 TEs produce 24 nt sRNA production in the seed. Future work comparing sRNA profiles in a ros1 mutant to those of wild-type could inform our understanding of TE spreading in a ros1 mutant, as would a comprehensive analysis of TE expression, again in both a ros1 mutant and in wild-type. It’s unclear to us if the endosperm would be the most informative or useful tissue to perform such analyses in.

Fig6 - as noted in the text, one way to test whether demethylation by ROS1 occurs before or after fertilization is to provide functional ROS1 through only one parent via reciprocal WT x ros-1 crosses, so that the endosperm always has ROS1 but either sperm or central cell does not, and see if this can rescue the paternal hypermethylation. If ROS1 acts prior to fertilization, then paternal ROS1 will rescue ros1 hypermethylation, but maternal ROS1 won't. If after fertilization, then either maternally or paternally supplied ROS1 will rescue the hypermethylation phenotype (assuming both are well expressed). Thus, to distinguish the two, it is sufficient to test whether maternally supplied ROS1 in an otherwise mutant background can rescue the hypermethylation phenotype, which is what is shown in Fig. 6. However, I think it's also important to show that paternally supplied ROS1 can also rescue the hypermethylation phenotype, which is not currently shown. The plots showing no effect on maternal mCG aren't as informative, since maternal methylation levels are mostly unaffected by ros1 anyway. Instead of comparing pairs of samples in a scatterplot, it might be clearer to show paternal mCG across all four comparisons (WT x WT, WT x ros1, ros1 x WT, and ros1 x ros1) side by side in a heatmap, using clustering to group similar behavior.

Response: We have revised this figure, now Fig. 8, in the following ways, which we believe addresses your comments and clarify the main conclusions (see same response to reviewer 2 for point 14):

Figure 8B remains as a scatterplot, where we observe significant correlation between individual ROS1 paternal, DME maternal regions in homozygous ros1 endosperm and heterozygous ros1/+ endosperm. Note that paternal allele methylation is higher in homozygous ros1 endosperm for most regions.

Figure 8C is now a boxplot comparing methylation of the paternal allele of ROS1 paternal, DME maternal regions (previously referred to as biallelically-demethylated, ROS1-dependent regions) across endosperm ROS1 genotypes. This plot shows increased methylation of paternal alleles when the paternal parent is a ros1 mutant, regardless of whether the resultant F1 endosperm is homozygous or heterozygous for ros1 (columns 3, 4, 6).

I would also suggest including a little more information in the main plots rather than only in the figure legends. For example, in Fig 2 including a label of 'ROS1-associated TE' for the two plots on the left, and 'TEs not associated with ROS1' on the right. Or for example in Fig. 3a indicating 'ros1-3 CG hyperDMRs' somewhere on the plot. This would just help make the figures easier to read at a glance. Please add common gene names to figures, instead just the ATG gene ID (Fig. S1a).

Response: Thank you for this feedback, we have made the suggested edits and additional edits of a similar nature.

Minor:<br /> - Fig. 1E is referenced in the text before Fig. 1D<br /> - Fig. S4 and S5 - there are more lines in the plot than the 6 genotypes listed in the legend, do these represent different replicates? If so that should be noted in the legend<br /> - Fig. 1B has no color legend for the different methylation sequence contexts (looks like same as 1A,C but should indicate either in plot or legend)<br /> - Line 42 should be "correspond to TE ends"<br /> - Line 93 "Based on previous studies..." should have references to those studies<br /> - When referring to the protein (rather than the genetic locus or mutant), ROS1 should not be italicized - for example line 130<br /> - Line 150 "we conclude that the loss"<br /> - Should add a y=x line to scatterplots, like those in Fig. 6<br /> - In fig. 1d, it's hard to evaluate the significance of the overlap of ROS1 targets with genes and TEs. Comparing these numbers to a control where the ROS1 targets have been randomly shuffled would help.

Response: We have made edits and additions where requested.

Reviewer #3 (Significance):

In this work, Hemenway and Gehring explore whether ROS1, DML2 and DML3 also affect DNA methylation patterns in endosperm. Using EM-seq of sorted endosperm nuclei, they show that loss of ROS1 indeed causes hypermethylation of a number of loci, particularly the flanks of methylated transposons, while loss of DML2 and DML3 has minimal additional effect. By obtaining allele-specific EM-seq data through crosses of Col and C24, the authors show that ros1 endosperm hypermethylation is mostly restricted to the paternal allele. The authors propose that at some sites, ROS1 helps bring down paternal methylation levels to match maternal methylation levels, which are typically reduced in endosperm due to DME activity in the female gametophyte prior to fertilization. In a ros1 mutant with paternal hypermethylation, these sites become differentially methylated on the maternal and paternal alleles, resembling imprinted loci. This work convincingly establishes a function for ROS1 in DNA methylation patterning in endosperm. However, I struggled with the clarity of the writing and reasoning in a few places, and would suggest clarification of a few points and additional analyses.

Response: Thank you for your comments. We have worked on streamlining the text and analysis.

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

DNA demethylases play a key role in DNA methylation patterning during flowering plant reproduction. The demethylase DME, in particular, is critical for proper endosperm development. While the function of DME in endosperm development has been explored, the contributions of the other demethylases in the same family, ROS1, DML2 and DML3 in Arabidopsis, have not yet been investigated. In vegetative tissues, ROS1 prevents hypermethylation of some loci. In this work, Hemenway and Gehring explore whether ROS1, DML2 and DML3 also affect DNA methylation patterns in endosperm. Using EM-seq of sorted endosperm nuclei, they show that loss of ROS1 indeed causes hypermethylation of a number of loci, particularly the flanks of methylated transposons, while loss of DML2 and DML3 has minimal additional effect. By obtaining allele-specific EM-seq data through crosses of Col and C24, the authors show that ros1 endosperm hypermethylation is mostly restricted to the paternal allele. The authors propose that at some sites, ROS1 helps bring down paternal methylation levels to match maternal methylation levels, which are typically reduced in endosperm due to DME activity in the female gametophyte prior to fertilization. In a ros1 mutant with paternal hypermethylation, these sites become differentially methylated on the maternal and paternal alleles, resembling imprinted loci. This work convincingly establishes a function for ROS1 in DNA methylation patterning in endosperm. However, I struggled with the clarity of the writing and reasoning in a few places, and would suggest clarification of a few points and additional analyses below.

I think making a few simple changes to streamline nomenclature would improve readability. For example, in the section starting on line 129, the same set of genomic features are called ROS1 target-proximal TEs, TEs that are near a ROS1 target region, and ROS1 target-associated TE regions. Also for example in line 254 "regions that are maternally-demethylated in wild-type endosperm, and are not dependent on ROS1 for proper demethylation" - are these the same as the "ROS1-independent, maternally-demethylated" regions in Fig. 5a? Given how complex these terms are, being consistent throughout the manuscript really helps the reader.

Is there any notable effect of ros1 on gene expression in endosperm? Endosperm is a terminal tissue, so maintaining DNA methylation boundaries as ROS1 does in vegetative tissues seems less important. It begs the question of why ROS1 is doing this in endosperm, is it just because it's there, or is there an endosperm-specific function? Exploring effects on imprinting would be particularly interesting (does loss of ROS1 'create' imprinted loci at these newly asymmetrically methylated sites?) but probably beyond the scope of the present work.

Is DME expressed in sperm, or is expression of DME affected in ros1 sperm or endosperm? One other explanation for ros1 hypermethylation occurring primarily on the paternal allele is that, potentially, DME can substitute for ROS1 in the central cell where DME is already very active, but not in sperm cells. Related, how well expressed is ROS1 vs. DME in sperm cells?

Fig 2b shows that ROS1 target-associated TEs are enriched for sRNAs in endosperm relative to embryo, whereas the reverse is true for non-ROS1-assoc TEs. Since TEs are not always well annotated and some may be missing from this analysis, what about trying the reverse analysis - are regions enriched for 24nt sRNAs in endosperm significantly hypermethylated in ros1 endosperm? All regions or only some?

What is the relationship between previously-defined DME targets and ROS1 targets identified in this paper? DME tends to target small euchromatic TE bodies, whereas Fig. 3 suggests that ROS1 helps prevent methylation spreading on the outer edges of the TEs, rather than in the TE body. Do all DME targets tend to be adjacent to or flanked by ROS1 target sites? Or are the TEs affected by DME (in body) and by ROS1 (at edges) largely nonoverlapping? Fig. 5a suggests that the ROS1-dependent, biallelically-demethylated sites are both DME and ROS1 targets, but how often do these really appear to overlap? More than by chance?

Are the TEs whose boundaries are demethylated by ROS1 more likely to be expressed in vegetative or endosperm tissues than TEs not affected by loss of ROS1? Expressed TEs likely produce more sRNAs, which would increase RdDM in a way that might need to be more actively countered by ROS1 than transcriptionally silent or evolutionarily older TEs.

Fig6 - as noted in the text, one way to test whether demethylation by ROS1 occurs before or after fertilization is to provide functional ROS1 through only one parent via reciprocal WT x ros-1 crosses, so that the endosperm always has ROS1 but either sperm or central cell does not, and see if this can rescue the paternal hypermethylation. If ROS1 acts prior to fertilization, then paternal ROS1 will rescue ros1 hypermethylation, but maternal ROS1 won't. If after fertilization, then either maternally or paternally supplied ROS1 will rescue the hypermethylation phenotype (assuming both are well expressed). Thus, to distinguish the two, it is sufficient to test whether maternally supplied ROS1 in an otherwise mutant background can rescue the hypermethylation phenotype, which is what is shown in Fig. 6. However, I think it's also important to show that paternally supplied ROS1 can also rescue the hypermethylation phenotype, which is not currently shown. The plots showing no effect on maternal mCG aren't as informative, since maternal methylation levels are mostly unaffected by ros1 anyway. Instead of comparing pairs of samples in a scatterplot, it might be clearer to show paternal mCG across all four comparisons (WT x WT, WT x ros1, ros1 x WT, and ros1 x ros1) side by side in a heatmap, using clustering to group similar behavior.

I would also suggest including a little more information in the main plots rather than only in the figure legends. For example, in Fig 2 including a label of 'ROS1-associated TE' for the two plots on the left, and 'TEs not associated with ROS1' on the right. Or for example in Fig. 3a indicating 'ros1-3 CG hyperDMRs' somewhere on the plot. This would just help make the figures easier to read at a glance. Please add common gene names to figures, instead just the ATG gene ID (Fig. S1a).

Minor:

• Fig. 1E is referenced in the text before Fig. 1D
• Fig. S4 and S5 - there are more lines in the plot than the 6 genotypes listed in the legend, do these represent different replicates? If so that should be noted in the legend
• Fig. 1B has no color legend for the different methylation sequence contexts (looks like same as 1A,C but should indicate either in plot or legend)
• Line 42 should be "correspond to TE ends"
• Line 93 "Based on previous studies..." should have references to those studies
• When referring to the protein (rather than the genetic locus or mutant), ROS1 should not be italicized - for example line 130
• Line 150 "we conclude that the loss"
• Should add a y=x line to scatterplots, like those in Fig. 6
• In fig. 1d, it's hard to evaluate the significance of the overlap of ROS1 targets with genes and TEs. Comparing these numbers to a control where the ROS1 targets have been randomly shuffled would help.

Significance

In this work, Hemenway and Gehring explore whether ROS1, DML2 and DML3 also affect DNA methylation patterns in endosperm. Using EM-seq of sorted endosperm nuclei, they show that loss of ROS1 indeed causes hypermethylation of a number of loci, particularly the flanks of methylated transposons, while loss of DML2 and DML3 has minimal additional effect. By obtaining allele-specific EM-seq data through crosses of Col and C24, the authors show that ros1 endosperm hypermethylation is mostly restricted to the paternal allele. The authors propose that at some sites, ROS1 helps bring down paternal methylation levels to match maternal methylation levels, which are typically reduced in endosperm due to DME activity in the female gametophyte prior to fertilization. In a ros1 mutant with paternal hypermethylation, these sites become differentially methylated on the maternal and paternal alleles, resembling imprinted loci. This work convincingly establishes a function for ROS1 in DNA methylation patterning in endosperm. However, I struggled with the clarity of the writing and reasoning in a few places, and would suggest clarification of a few points and additional analyses

Add IPFS! Desktop

Can do it manually

DOI: 10.1007/s13402-025-01089-y

Resource: (RRID:CVCL_0297)

Curator: @dhovakimyan1

SciCrunch record: RRID:CVCL_0297

What is this?

DOI: 10.1007/s13402-025-01089-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:CVCL_C8XN

What is this?

DOI: 10.1007/s13402-025-01089-y

Resource: (NCI-DTP Cat# OVCAR-8, RRID:CVCL_1629)

Curator: @scibot

SciCrunch record: RRID:CVCL_1629

What is this?

DOI: 10.1007/s13402-025-01089-y

Resource: (RRID:CVCL_0063)

Curator: @scibot

SciCrunch record: RRID:CVCL_0063

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: RRID:Addgene_104587

Curator: @scibot

SciCrunch record: RRID:Addgene_104587

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: (ATCC Cat# PTA-5077, RRID:CVCL_6911)

Curator: @scibot

SciCrunch record: RRID:CVCL_6911

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: RRID:Addgene_12259

Curator: @scibot

SciCrunch record: RRID:Addgene_12259

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: RRID:Addgene_64222

Curator: @scibot

SciCrunch record: RRID:Addgene_64222

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: RRID:Addgene_12260

Curator: @scibot

SciCrunch record: RRID:Addgene_12260

What is this?

DOI: 10.1038/s41523-025-00795-y

Resource: RRID:Addgene_51007

Curator: @scibot

SciCrunch record: RRID:Addgene_51007

What is this?

DOI: 10.1016/j.celrep.2025.116093

Resource: (BioLegend Cat# 107619, RRID:AB_493528)

Curator: @scibot

SciCrunch record: RRID:AB_493528

What is this?

DOI: 10.1007/s00401-025-02918-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_2802618

What is this?

DOI: 10.1007/s00401-025-02918-y

Resource: (Millipore Cat# MABN389, RRID:AB_2716647)

Curator: @scibot

SciCrunch record: RRID:AB_2716647

What is this?

DOI: 10.1007/s00401-025-02918-y

Resource: (Abcam Cat# ab108986, RRID:AB_10891773)

Curator: @scibot

SciCrunch record: RRID:AB_10891773

What is this?

If a wave is moving in one direction it is called a plane wave (dependent on only x & t OR y & t)

Reviewer #2 (Public review):

Summary:

The goal of this study was to find evidence for local adaptation in survival and fecundity of the model plant Arabidopsis thaliana. The authors grew a large set of Swedish Arabidopsis accessions at four common garden sites in northern and southern Sweden. Accessions were grown from seed in trays, which were laid on the ground at each site in late summer, screened for survival in fall and the following spring, and fecundity was determined from rosette size and seed production in spring. Experiments were complemented by 'selection experiments', in which seeds of the same accessions were sown in plots, and after two years of growth, plants were sampled to determine fitness from genotype frequencies, providing a more comprehensive evaluation of lifetime fitness than can be gleaned from fecundity alone.

As the main result, southern accessions had higher mortality in northern sites in one of two years, but also suffered more slug damage in southern sites in one year, indicating a potential link between frost tolerance and herbivore resistance. Fecundity of accession was highest when growing close to the 'home' environment, but while accessions from one sand dune population in southern Sweden had among the lowest fecundities overall, they consistently had the highest fitness in the selection experiment. Accessions from this population had large seed size and rapid root growth, which might be related to establishment success when arriving in a new, partially occupied habitat. However, neither trait could fully explain the very high fitness of this population, suggesting the presence of other, unmeasured traits.

Overall, the authors could provide clear evidence of local adaptation in different traits for some of their experiments, but they also highlight high temporal and spatial variability that makes prediction of microevolutionary change so challenging.

Strengths:

A major strength of this study is the highly comprehensive evaluation of different fitness-related traits of Arabidopsis under natural conditions. The evaluation of survival and fecundity in common garden experiments across four sites and two years provides an estimate of variability and consistency of results. The addition of the 'selection experiment' provides an extended view on plant fitness that is both original and interesting, in particular highlighting potential limitations of 'fitness-proxies' such as seed production that don't take into account seedling establishment and competitive exclusion.

Throughout the study, the authors have gone to impressive depths in exploring their data, and particularly the discovery of 'native volunteers' in selection experiment plots and their statistical treatment is very elegant and has resulted in compelling conclusions. Also, while the authors are careful in the interpretation of their GWAS results, they nonetheless highlight a few interesting gene candidates that may be underlying the observed plant adaptations, and which likely will stimulate further research.

Overall, the authors provide a rich new resource that is relevant and interesting both in the context of general evolutionary theory as well as more specifically for molecular biology.

Weaknesses:

While the repetition of the common garden experiments over two years is certainly better than no repetition (hence its mention also under 'strengths'), the very high variability found between the two years highlights the need for more extensive temporal replication. In this context, two temporal replicates are the bare minimum, and more repeats in time would be necessary to draw any kind of conclusion about the role of 'high mortality' and 'low mortality' years for the microevolution of Arabidopsis. It also seems that the authors missed an opportunity to explore potentially causal variation among years, as they did not attempt to relate winter mortality to actual climatic variables, even though they discuss winter harshness as a potential predictor.

The low temporal variation also makes the accidental slug herbivory appear somewhat random. Potted plants are notoriously susceptible to slug herbivory, and while it is certainly nice that slug damage predominantly affected one group of accessions, it nonetheless raises the question whether this reflects a 'real' selection pressure that plants commonly face in their respective local environments.

The addition of the 'selection experiment' is certainly original and provides valuable additional insights, but again, it seems a bit questionable which natural process really has affected this outcome. While the genetic and statistical analysis of this experiment seems to be state-of-the-art, the experimental design is rather rudimentary compared to more standard selection experiments. Specifically, the authors added seeds from greenhouse-grown mothers to experimental plots and only sampled plants two years later. This means that, potentially,y the first very big bottleneck was germination under natural conditions, which may have already excluded many of the accessions before they had a chance to grow. While this certainly is one type of selection, it is not exactly the type of selection that a 2-year selection experiment is set up to measure. Either initially establishing the selection experiment from plants instead of seeds, or genotyping the population over several generations, would have substantially strengthened the conclusions that could be drawn from this experiment. Also, the complete lack of information on population density is a bit problematic. It is not clear if there were other (non-Arabidopsis) plants present in the plots, how many Arabidopsis plants were established, if numbers changed over the year, etc. Given all of these limitations, calling this a 'selection experiment' is in fact somewhat misleading.

Despite these weaknesses, the authors could achieve their main goals, and despite the somewhat minimal temporal replication, they were lucky to sample two fairly distinct years that provided them with interesting variation, which they could partially explain using the variation among their accessions. Overall, this study will likely make an important contribution to the field of evolutionary biology, and it is another very strong example of how the extensive molecular tools in Arabidopsis can be leveraged to address fundamental questions in evolution and ecology, to an extent that is not (yet) possible in other plant systems.

Author response:

Reviewer #1 (Public review):

Summary: 

As a general phenomenon, adaptation of populations to their respective local conditions is well-documented, though not universally. In particular, local adaptation has been amply demonstrated in Arabidopsis thaliana, the focal species of this research, which is naturally highly selfing. Here, the authors report assays designed to evaluate the spatial scale of fitness variation among source populations and sites, as well as temporal variability in fitness expression. Further, they endeavor to identify traits and genomic regions that contribute to the demonstrated variation in fitness.  

Strengths: 

With many (200) inbred accessions drawn from throughout Sweden, the study offers an unusually fine sampling of genetic variation within this much-studied species, and through assays in multiple sites and years, it amply demonstrates the context-dependence of fitness expression. It supports the general phenomenon of local adaptation, with multiple nuances. Other examples exist, but it is of value to have further cases illustrating not only the context-dependence of fitness expression but also the sometimes idiosyncratic nature of fitness variation. I commend the authors on their cautionary language in relation to inferences about the roles of particular genomic regions (e.g.l.140-144; l.227)  

Weaknesses: 

To my mind, the manuscript is written primarily for the Arabidopsis community. This community is certainly large, but there are many evolutionary biologists who could appreciate this work but are not invited to do so. The authors could address the broader evolution community by acknowledging more of the relevant work of others (I've noted a few references in my comments to the authors). At least as important, the authors could make clearer the fact that A. thaliana is (almost) strictly selfing and how this feature of its biology both enables such a study and also limits inferences from it. Further, it seems to me that though I could be wrong, readers would appreciate a more direct, less discursive style of writing, and one that makes the broader import of the focal questions clearer. 

we agree that connecting the paper better to the broader field is desirable, and will try to do this in the revision. As for how selfing matters, there certainly are some things we can discuss, but a general discussion is probably a suitable topic for a review/opinion article!

As a reader, I would value seeing estimates of the overall fitness of the accessions in the different conditions, i.e., by combining the survival and fecundity results of the common garden experiments.

Combining estimates would be possible in the common garden experiments, and would bring us somewhat closer to total fitness estimates, although as noted by another reviewer (and also emphasized by us), the time scale of our experiment is not sufficient to evaluate the trade-off between survival and fecundity. Furthermore, we would still be missing the establishment component of fitness, which we found to be extremely important. Therefore little would be gained by combining the estimates, while at the same time losing resolution to disentangle the fitness components. We thus decided to focus on the individual fitness components and leave consideration of their joint effect for the Discussion.

Reviewer #2 (Public review):

Summary: 

The goal of this study was to find evidence for local adaptation in survival and fecundity of the model plant Arabidopsis thaliana. The authors grew a large set of Swedish Arabidopsis accessions at four common garden sites in northern and southern Sweden. Accessions were grown from seed in trays, which were laid on the ground at each site in late summer, screened for survival in fall and the following spring, and fecundity was determined from rosette size and seed production in spring. Experiments were complemented by 'selection experiments', in which seeds of the same accessions were sown in plots, and after two years of growth, plants were sampled to determine fitness from genotype frequencies, providing a more comprehensive evaluation of lifetime fitness than can be gleaned from fecundity alone. 

To clarify, fecundity was determined from total plant area using photos of the mature stems, not the rosettes or direct counting of seeds. That said, it is true that our fecundity estimate was well correlated with rosette area. Furthermore, we validate our fecundity estimates by showing they were highly correlated with seed production estimated by measuring and counting siliques on a separate set of plants grown under common garden conditions in one of our sites (Brachi et al.2022). 

As the main result, southern accessions had higher mortality in northern sites in one of two years, but also suffered more slug damage in southern sites in one year, indicating a potential link between frost tolerance and herbivore resistance. Fecundity of accession was highest when growing close to the 'home' environment, but while accessions from one sand dune population in southern Sweden had among the lowest fecundities overall, they consistently had the highest fitness in the selection experiment. Accessions from this population had large seed size and rapid root growth, which might be related to establishment success when arriving in a new, partially occupied habitat. However, neither trait could fully explain the very high fitness of this population, suggesting the presence of other, unmeasured traits. 

Overall, the authors could provide clear evidence of local adaptation in different traits for some of their experiments, but they also highlight high temporal and spatial variability that makes prediction of microevolutionary change so challenging. 

Strengths: 

A major strength of this study is the highly comprehensive evaluation of different fitness-related traits of Arabidopsis under natural conditions. The evaluation of survival and fecundity in common garden experiments across four sites and two years provides an estimate of variability and consistency of results. The addition of the 'selection experiment' provides an extended view on plant fitness that is both original and interesting, in particular highlighting potential limitations of 'fitness-proxies' such as seed production that don't take into account seedling establishment and competitive exclusion. 

Throughout the study, the authors have gone to impressive depths in exploring their data, and particularly the discovery of 'native volunteers' in selection experiment plots and their statistical treatment is very elegant and has resulted in compelling conclusions. Also, while the authors are careful in the interpretation of their GWAS results, they nonetheless highlight a few interesting gene candidates that may be underlying the observed plant adaptations, and which likely will stimulate further research. 

Overall, the authors provide a rich new resource that is relevant and interesting both in the context of general evolutionary theory as well as more specifically for molecular biology. 

Weaknesses:

While the repetition of the common garden experiments over two years is certainly better than no repetition (hence its mention also under 'strengths'), the very high variability found between the two years highlights the need for more extensive temporal replication. In this context, two temporal replicates are the bare minimum, and more repeats in time would be necessary to draw any kind of conclusion about the role of 'high mortality' and 'low mortality' years for the microevolution of Arabidopsis. It also seems that the authors missed an opportunity to explore potentially causal variation among years, as they did not attempt to relate winter mortality to actual climatic variables, even though they discuss winter harshness as a potential predictor.

We agree that two years is insufficient to understand how variation in selective pressures compound over time to generate micro-evolutionary change. The eight-year data in Oakley et al. (2023), which we discuss in the paper, support this. Our results are nonetheless sufficient to demonstrate the idiosyncratic nature of selection. In the revision, we will further emphasize that far longer time series would be needed for definitive conclusions.

Our short time series is also why we do not try to correlate with climate data, as this would amount to doing statistics with four data points (mostly two groups of accession N vs S, with mostly homogenous climates within groups, and two years).

The low temporal variation also makes the accidental slug herbivory appear somewhat random. Potted plants are notoriously susceptible to slug herbivory, and while it is certainly nice that slug damage predominantly affected one group of accessions, it nonetheless raises the question whether this reflects a 'real' selection pressure that plants commonly face in their respective local environments. 

We agree with this point as well. The evidence for selection on glucosinolates by generalist herbivores such as slugs is fairly strong, but the precise agent is not known, and probably varies over time and space. Our results merely demonstrate one possibility (and we will clarify this in the revision).

The addition of the 'selection experiment' is certainly original and provides valuable additional insights, but again, it seems a bit questionable which natural process really has affected this outcome. While the genetic and statistical analysis of this experiment seems to be state-of-the-art, the experimental design is rather rudimentary compared to more standard selection experiments. Specifically, the authors added seeds from greenhouse-grown mothers to experimental plots and only sampled plants two years later. This means that, potentially,y the first very big bottleneck was germination under natural conditions, which may have already excluded many of the accessions before they had a chance to grow. While this certainly is one type of selection, it is not exactly the type of selection that a 2-year selection experiment is set up to measure. Either initially establishing the selection experiment from plants instead of seeds, or genotyping the population over several generations, would have substantially strengthened the conclusions that could be drawn from this experiment.

We agree that more data would have been beneficial, and we do not make strong claims about the nature of selection. Among other phenotypes, we mention dormancy, and note that existing dormancy estimates do not predict fitness in our selection experiments. In addition the same seed batches germinated uniformly in the common-garden experiments with minimal stratification (we will note this in the revision).

Also, the complete lack of information on population density is a bit problematic. It is not clear if there were other (non-Arabidopsis) plants present in the plots, how many Arabidopsis plants were established, if numbers changed over the year, etc. Given all of these limitations, calling this a 'selection experiment' is in fact somewhat misleading. 

Seeds were introduced into sites that appeared appropriate for A. thaliana, leaving the background community intact. We provided information on sowing density; the density of plants (A. thaliana and other species) that we obtained during the course of the experiments varied considerably between sites, much like in natural populations, although we lack systematic measurements. We will provide more information (including photos) in the revision.  

Despite these weaknesses, the authors could achieve their main goals, and despite the somewhat minimal temporal replication, they were lucky to sample two fairly distinct years that provided them with interesting variation, which they could partially explain using the variation among their accessions. Overall, this study will likely make an important contribution to the field of evolutionary biology, and it is another very strong example of how the extensive molecular tools in Arabidopsis can be leveraged to address fundamental questions in evolution and ecology, to an extent that is not (yet) possible in other plant systems. 

Reviewer #3 (Public review)

Summary: 

The manuscript presents a large common garden experiment across Sweden using solely local germplasm. Additionally, there is a collection of selection experiments that begin investigating the factors shaping fecundity in these populations. This provides an impressive amount of data and analysis investigating the underlying factors involved. Together, this helps support the data showing that fluctuations and interactions are key components determining Arabidopsis fitness and are more broadly applicable across plant and non-plant species. 

Strengths: 

The field trials are well conducted with extensive effort and sampling. Similarly while the genetic analysis is complex it is well conducted and reflects the complexity of dealing with population structure that may be intricately linked to adaptive structure. This has no real solution and the option of presenting results with and without correction is likely the only appropriate option. 

Weaknesses: 

A significant finding from this study was that fecundity is shaped more by yearly fluctuations and their interaction with genotype than it is by the main effect of location or genotype. Another significant finding is that the strength of selection can be quite strong, with nearly 5x ranges across accessions. It should be noted that there are a number of other studies using Arabidopsis in the wild with multiple years and locations that found similar observations beyond the Oakley citation. In general, the context of how these findings relate to existing knowledge in Arabidopsis is a bit underdeveloped. 

We shall remedy this in the revision (see also comments by Reviewer #1).

The effects of the populations across the locations seem to rely on individual tests and PC analysis. It would seem to be possible to incorporate these tests more directly in the linear modeling analysis, and it isn't quite clear why this wasn't conducted. 

The fecundity estimates were modelled for all experiments simultaneously and the results are presented in Figure 6 to explore the relative importance of genotype effects and interaction terms including genotypes. For survival and fecundity, the BLUPS are generated from linear mixed models fitted for all experiments simultaneously including a random intercept effect for the genotypes within experiments. A principal component analysis is used to explore the pattern of accession effects (BLUPS) on fecundity (Figure 7); this will be explained in the Methods.  

I'm a bit puzzled by the discussion on how to find causative loci. This seems to focus solely on GWAS as the solution, with a goal to sequence vast individuals. But the loci that the manuscript discussed were found by a combination of structured mapping populations followed by molecular validation that then informed the GWAS. As such, I'm unsure if the proposed future approach of more sequencing is the best when a more balanced approach integrating diverse methods and population types will be more useful. 

We are puzzled by this comment in return. Our statement about more sequencing (penultimate sentence of discussion) was referring to achieving a better understanding of the history of migration and selection rather than identifying causative loci. Happy for clarification!

References

Brachi, Benjamin, Daniele Filiault, Hannah Whitehurst, Paul Darme, Pierre Le Gars, Marine Le Mentec, Timothy C. Morton, et al. 2022. “Plant Genetic Effects on Microbial Hubs Impact Host Fitness in Repeated Field Trials.” Proceedings of the National Academy of Sciences of the United States of America 119 (30): e2201285119.

Oakley, Christopher G., Douglas W. Schemske, John K. McKay, and Jon Ågren. 2023. “Ecological Genetics of Local Adaptation in Arabidopsis: An 8-Year Field Experiment.” Molecular Ecology, June. https://doi.org/10.1111/mec.17045.

This work has been peer reviewed in GigaScience (see https://doi.org/10.1093/gigascience/giaf076), which carries out open, named peer-review. These reviews are published under a CC-BY 4.0 license and were as follows:

Reviewer name: Daniel A. Skelly

Overall, this is a very nice writeup of a useful package that extends the Seurat package to expand possibilities for single cell analysts in R. I liked the visualization options, the ability to try certain python-based tools easily in R which was not previously easy, and some of the authors' new innovations like their use of pathway enrichment scores in broad ways. Kudos to the authors for releasing a package with really excellent documentation and tutorials!

I think this paper could be made better if the authors stressed with a little more clarity how specifically their work is innovative. The text in the present manuscript is fine but reads like a bit of a grab bag of functionality. For example, from the abstract: "SeuratExtend offers a user-friendly and intuitive interface for performing a wide range of analyses, including functional enrichment, trajectory inference, gene regulatory network reconstruction, and denoising. The package integrates multiple databases, … and incorporates popular Python tools … [We] showcase its novel applications in pathway-level analysis and cluster annotation. SeuratExtend enhances data visualization …"

How could they be more clear or specific? One example could be by categorizing what SeuratExtend can do that other packages can't. For example, I see innovations in perhaps three general areas: 1. Making single cell analyses easier/faster/prettier (i.e. visualizations, pathway enrichment) 2. Making previously published single cell tools more broadly accessible (e.g. first option to bring certain python tools to R) 3. New innovations (e.g. dimensionality reduction and clustering based on pathway enrichment scores; may not be completely new but I don't recall seeing this elsewhere) If this was added I feel the paper would more clearly communicate to readers the information necessary for them to choose whether they want to try the package.

I have the following additional significant comments: * Integration of multiple databases for GSEA — these methods are good, but what about in a few years when those databases have been updated? Do the authors intend to continue updating? Could they provide a function for users to use their own database (e.g. .gaf and .obo files, for example for another model organism)? Similar comment about gene identifer conversion, which may need to be updated every few years. * "While the Python ecosystem has benefited greatly from the comprehensive scverse project [7], which utilizes the universal AnnData format to connect various tools and algorithms, a comparable integrated solution has been lacking in the R community. SeuratExtend addresses this gap by providing a unified framework centered around the Seurat object, effectively becoming the R counterpart to scverse." —> some might argue that SeuratWrappers is this solution. The authors should more clearly and explicitly comment on what SeuratExtend does differently/better than SeuratWrappers. * I'm not particularly convinced by the authors' example studies that used SeuratExtend. For example, they describe Hua-Vella et al. (2022) and Hua et al. (2023). These are very nice studies and I have no doubt they made use of SeuratExtend in their analyses. But I don't see anything these authors describe those authors doing as being uniquely possible with SeuratExtend. Perhaps SeuratExtend made their analyses easier, or faster. But it would be better if we had some further concrete details. For example, something communicating a message like one of the following: (1) the authors only tested method X on a whim because it was so easy to run in SeuratExtend, and found that it revealed unexpected biology Y; or (2) the authors were able to bring together method X which runs in R and method Y which runs in python and the joint inference — not possible in other packages — revealed key result Z. If the authors of this manuscript can't point to those sorts of examples, then I'm not sure it adds much to include this discussion in the present paper. * I really liked the section "Novel Applications of SeuratExtend in Pathway-Level Analysis and Cluster Annotation", especially "Exploring and Analyzing Single-Cell Data at the Pathway Level". I thought these applications could perhaps be stressed a bit more strongly or made more prominent earlier in the paper. * Figures 2 and 3 are showing example plots from which we don't actually need to infer any important biology. I thought these figures could be combined and each individual plot type only shown once. (This is for clarity and I don't see anything incorrect about the authors' current plots. * There may be some issues with dependencies for some users. For example, it prompted me to install viridis and loomR as I went through the Quickstart. I ended up encountering an error there is no package called 'loomR' while trying. I had to manually install with remotes::install_github(repo = "mojaveazure/loomR"). Maybe provide an explicit dependencies list/list of recommended packages to install? * I had an error the first time calling Palantir.RunDM(). I hadn't created a seuratextend environment. I found that I could do this manually using create_condaenv_seuratextend(), but that this wasn't supported for Apple Silicon chips. I would suggest that the authors do try to find a way to get this working on newer Apple chips, because Mac machines are very common among bioinformaticians in my experience. * While the writing is largely quite clear, I found it to be a bit voluminous. If the authors are able to cut down on text length that may help in emphasizing the key points that make their package valuable to users.

I had these minor comments: * "Moreover, mainstream scRNA-seq analysis tools are primarily developed for either the R or Python platforms, with additional options like Nextflow and Snakemake" — I suggest revising this sentence. The tools are developed in R or python languages, which I would not call platforms. I would reword that Nextflow and Snakemake are workflow management systems that provide additional options for pipeline automation * "the R ecosystem surrounding Seurat appears relatively limited" — I'm not sure I would agree with this. I counted wrappers for 17 methods currently. Yes it is true that there are more packages in scverse. However, I suggest moderating your claims about Seurat being limited. * Suggest removing snakemake from Table 1 — it is really different from the other tools listed there

hay que poner una referencia acá, el lector no tiene por que estar enterado de esto

justificar esta elección

explicar como y con que criterio fue realizado el suavizado

como que no? abajo tenes series están planchadas y luego suben y tenes una que oscila

los títulos en los mapas lideres político y fuerzas armadas se pegan, emprolijar

si le sacas el titulo a la figura y lo pones en el caption le das mas lugar a los mapas

esto no debería estar en la parte en donde hablas de bert topic? Para ser más claro, conviene usar un párrafo para un tema, usar un párrafo para otro tema, y si tenes que compara dos temas empezar el diciendo algo así como "comparando entre los dos métodos ... "

punto y aparte

sale el ":", entra una "," y un punto al final de la ecuación

creo que en ningún momento se aclaro como es el sistema electoral en USA, chequear

Esto esta algo desprolijo, mejor algo más parecido a esto : "Los distintos paneles en la figura muestran la relación entre la cantidad de discursos recopilados (eje horizontal) y el "swing score" en los distintos estados para diferentes ventanas temporales"

esto es importante, sacarlo de adentro del paréntesis y darle el espació necesario en el texto

punto y aparte

no se entiende, se más claro y directo

citar referencia

dicho así parece que es la única alternativa y sabemos que hay otros approachs que vienen de las ciencias sociales igualmente validos. De hecho no estaría mal citar algunos trabajos de otras áreas a modo de ejemplo

siempre que se escribe una afirmación de este estilo, donde no es trivial lo que se afirma, conviene introducir una cita para justificarla

In other word,\(y \not \in \text{FreeVariables}(c) \cup {x}\)

Specify the impact of the API, e.g., improved data retrieval speed by X% or reduced storage costs by Y%.

Author response:

The following is the authors’ response to the previous reviews

Reviewer #1 (Public review):

Summary:

The major result in the manuscript is the observation of the higher order structures in a cryoET reconstruction that could be used for understanding the assembly of toroid structures. The cross-linking ability of ZapD dimers result in bending of FtsZ filaments to a constant curvature. Many such short filaments are stitched together to form a toroid like structure. The geometry of assembly of filaments - whether they form straight bundles or toroid like structures - depends on the relative concentrations of FtsZ and ZapD.

Strengths:

In addition to a clear picture of the FtsZ assembly into ring-like structures, the authors have carried out basic biochemistry and biophysical techniques to assay the GTPase activity, the kinetics of assembly, and the ZapD to FtsZ ratio.

Weaknesses:

The discussion does not provide an overall perspective that correlates the cryoET structural organisation of filaments with the biophysical data. The current version has improved in terms of addressing this weakness and clearly states the lacuna in the model proposed based on the technical limitations.

Future scope of work includes the molecular basis of curvature generation and how molecular features of FtsZ and ZapD affect the membrane binding of the higher order assembly.

Reviewer #3 (Public review):

Summary:

Previous studies have analyzed the binding of ZapD to FtsZ and provided images of negatively stained toroids and straight bundles, where FtsZ filaments are presumably crosslinked by ZapD dimers. Toroids without ZapD have also been previously formed by treating FtsZ with crowding agents. The present study is the first to apply cryoEM tomography, which can resolve the structure of the toroids in 3D. This shows a complex mixture of filaments and sheets irregularly stacked in the Z direction and spaced radially. The most important interpretation would be to distinguish FtsZ filaments from ZapD crosslinks, This is less convincing. The authors seem aware of the ambiguity: "However, we were unable to obtain detailed structural information about the ZapD connectors due to the heterogeneity and density of the toroidal structures, which showed significant variability in the conformations of the connections between the filaments in all directions." Therefore, the reader may assume that the crosslinks identified and colored red are only suggestions, and look for their own structural interpretations. But readers should also note some inconsistencies in stoichiometry and crosslinking arrangements that are detailed under "weaknesses."

Strengths.

This is the first cryoEM tomography to image toroids and straight bundles of FtsZ filaments bound to ZapD. A strength is the resolution, which. at least for the straight bundles. is sufficient to resolve the ~4.5 nm spacing of ZapD dimers attached to and projecting subunits of an FtsZ filament. Another strength is the pelleting assay to determine the stoichiometry of ZapD:FtsZ (although this also leads to weaknesses of interpretation).

Weaknesses

The stoichiometry presents some problems. Fig. S5 uses pelleting to convincingly establish the stoichiometry of ZapD:FtsZ. Although ZapD is a dimer, the concentration of ZapD is always expressed as that of its subunit monomers. Fig. S5 shows the stoichiometry of ZapD:FtsZ to be 1:1 or 2:1 at equimolar or high concentrations of ZapD. Thus at equimolar ZapD, each ZapD dimer should bridge two FtsZ's, likely forming crosslinks between filaments. At high ZapD, each FtsZ should have it's own ZapD dimer. However, this seems contradicted by later statements in Discussion and Results. (1) "At lower concentrations of ZapD, .. toroids are the most prominent structures, containing one ZapD dimer for every four to six FtsZ molecules." Shouldn't it be one ZapD dimer for every two FtsZ? (2) "at the high ZapD concentration...a ZapD dimer binds two FtsZ molecules connecting two filaments." Doesn't Fig. S5 show that each FtsZ subunit has its own ZapD dimer? And wouldn't this saturate the CTD sites with dimers and thus minimize crosslinking?

We thank the reviewer for these insightful comments. The affinity of ZapD for FtsZ is relatively low and a higher concentration of ZapD is required in solution to effectively saturate the binding sites of all FtsZ molecules forming macrostructures. It is important to clarify that the concentrations mentioned in the text refer to the amounts and ratios of protein added to the total volume of the sample, rather than the proteins actively interacting and forming bundles or macrostructures.

To differentiate, two aspects can be considered: the ratio of added protein (as mentioned in the text) and the fraction of proteins that contribute to the formation of the macrostructures. Under polymerization conditions, FtsZ-GTP recruits additional monomers to form polymers. Therefore, more FtsZ than ZapD would be involved in forming filaments and bundles. Our results support this hypothesis and show that a higher amount of ZapD is required in the sample to pellet with FtsZ bundles.

We propose that starting with the same initial concentration of FtsZ and ZapD in solution, only a small fraction of ZapD will bind to the structures, favoring the formation of toroidal structures despite the initial 1:1 ratio of proteins added to the sample. When considering a higher FtsZ:ZapD ratio (1:6), the increased amount of ZapD in solution would facilitate the saturation of all FtsZ binding sites, consistent with the observation of straight bundles. Analytical sedimentation velocity data further supported this finding, indicating a binding ratio of approximately 0.3-0.4, suggesting that one ZapD dimer binds for every 4-6 FtsZ monomers. The binding ratio indicates that two FtsZ monomers will bind to a single dimer of ZapD, but this only occurs when there is a significant excess of ZapD over FtsZ in the solution mixture. 

These findings align qualitatively with the relative intensities of the electrophoretic bands observed for FtsZ and ZapD in the pelleting assay with different FtsZ-ZapD mixtures, as shown in Suppl. Fig. 5 as % of FtsZ in the fractions. Without prior staining calibration of the gels, there is no simple quantitative relationship between gel band intensities after Coomassie staining and the amount of protein in a band (Darawshe et al. 1993 Anal Biochem - DOI: 10.1006/abio.1993.1581). This last point precludes a quantitative comparison between pelleting/SDS-PAGE data and analytical sedimentation measurements. For this reason, we have decided to present pelleting results as % of FtsZ in supernatant and pellet to avoid overestimations. 

A major weakness is the interpretation of the cryoEM tomograms, specifically distinguishing ZapD from FtsZ. The distinction of crosslinks seems based primarily on structure: long continuous filaments (which often appear as sheets) are FtsZ, and small masses between filaments are ZapD. The density of crosslinks seems to vary substantially over different parts of the figures. More important, the density of ZapD's identified and colored red seem much lower than the stoichiometry detailed above. Since the mass of the ZapD monomer is half that of FtsZ, the 1:1 stoichiometry in toroids means that 1/3 of the mass should be ZapD and 2/3 FtsZ. However, the connections identified as ZapD seem much fewer than the expected 1/3 of the mass. The authors conclude that connections run horizontally, diagonally and vertically, which implies no regularity. This seems likely, but as I would suggest that readers need to consider for themselves what they would identify as a crosslink.

The amount of ZapD in the toroids will be significantly less than one third. Although the theoretical addition of protein to the samples is at a 1:1 ratio, the actual amount of protein in the macrostructures containing ZapD is much lower, as shown by sedimentation velocity pelleting assays.

In contrast to the toroids formed at equimolar FtsZ and ZapD, thin bundles of straight filaments are assembled in excess ZapD. Here the stoichiometry is 2:1, which would mean that every FtsZ should have a bound ZapD DIMER. The segmentation of a single filament in Fig. 5e seems to agree with this, showing an FtsZ filament with spikes emanating like a picket fence, with a 4.5 nm periodicity. This is consistent with each spike being a ZapD dimer, and every FtsZ subunit along the filament having a bound ZapD dimer. But if each FtsZ has its own dimer, this would seem to eliminate crosslinking. The interpretative diagram in Fig. 6, far right, which shows almost all ZapD dimers bridging two FtsZs on opposite filaments, would be inconsistent with this 2:1 stoichiometry.

Assessing the precise stoichiometry of FtsZ and ZapD within the macrostructures is challenging. We interpret the spikes as ZapD dimers bridging two FtsZ filaments, implying a theoretical 1:1 stoichiometry in the straight bundle. However, ZapD may be enriched in certain areas, indicating that a single FtsZ monomer is binding to one side of the dimer. In contrast, the other side remains available for additional connections, resulting in a potential 2:1 stoichiometry. A combination of both scenarios is likely, although our resolution does not allow further characterization. Considering these complexities, we assume these connections represent a dimer of ZapD binding to two FtsZ monomers.

Figure 6 shows a simplified scheme illustrating how the bundles could be assembled based on the Cryo-ET data. We acknowledge the limitations of this diagram; its purpose is to depict the mesh formed by the stabilization of ZapD. We have not included interactions that do not lead to filament crosslinking, such as dimers binding to only one FtsZ filament. This focus enhances the interpretation of the scheme and the FtsZ-ZapD interaction. A sentence has been added to the caption to highlight the possibility of other interactions not considered in the scheme.

In the original review I suggested a control that might help identify the structures of ZapD in the toroids. Popp et al (Biopolymers 2009) generated FtsZ toroids that were identical in size and shape to those here, but lacking ZapD. These toroids of pure FtsZ were generated by adding 8% polyvinyl chloride, a crowding agent. The filamentous substructure of these toroids in negative stain seemed very similar to that of the ZapD toroids here. CryoET of these toroids lacking ZapD might have been helpful in confirming the identification of ZapD crosslinks in the present toroids. However, the authors declined to explore this control.

The mechanisms by which methylcellulose (MC) promotes the assembly of FtsZ macrostructures reported by Popp et al. involve more than simple excluded volume effects, as the low concentration of MC (less than 1 mg/ml) falls below the typical crowding regime. The latter suggests the existence of poorly characterized additional interactions between MC and FtsZ. These complexities preclude the use of FtsZ polymers formed in the presence of MC as a true control for the FtsZ toroidal structures reported here.

Finally, it should be noted that the CTD binding sites for ZapD should be on the outside of curved filaments, the side facing the membrane in the cell. All bound ZapD should project radially outward, and if it contacted the back side of the next filament, it should not bind (because the CTD is on the front side). The diagram second to right in Fig. 6 seems to incorporate this abortive contact.

The role of the flexible linker and its biological implications are still under debate in the field. The flexible linker allows ZapD-driven connections to be made in different directions. While these implications are not the primary focus of our manuscript, the flexible linker could allow connections between filaments in different orientations.

Reviewer #1 (Recommendations for the authors):

Most of the concerns which I had raised in the earlier version have been taken care of, as detailed in the response.

A few minor points, mostly related to re-phrasing are listed below:

Page 2: line 21: The use of the term 'C-terminal domain' for the C-terminal unstructured region of FtsZ is confusing. The term C-terminal domain or CTD for FtsZ is commonly used to describe part of the globular domain, while C-terminal tail or CCTP will be a more apt usage for all the instances in this manuscript.

We refer to the C-terminal domain as the carboxy-terminal region of the protein. This domain includes the C-terminal linker (CTL), which varies in length between species, followed by a conserved 11-residue sequence (CTC) and shorter, variable C-terminal sequences (CTV). We used the term "C-terminal domain" primarily to improve the readability of the manuscript, but we appreciate the reviewer's feedback. We have now adopted the term "CCTP" instead of "C-terminal domain" to improve the clarity of our manuscript.

On a related note, the schematic in Fig 1 shows the interaction with CCTP rather than the C-terminal domain of the globular FtsZ. Please provide an explanation.

We refer to the unstructured C-terminal domain of FtsZ as the C-terminal tail. To avoid confusion, we have introduced the term CCTP in this manuscript.

Supple Fig 2: "The FCS analysis demonstrated an increasing diffusion time of ZapD along with the FtsZ concentration as result of higher proportion of ZapD bound to FtsZ.

The increased diffusion time need not be interpreted as increased ZapD bound, it could also mean that FtsZ could polymerise in the presence of increasing ZapD, was this possibility ruled out? Including a comment on this aspect will be useful.

In these experiments, we monitored fluorescently labeled ZapD. Due to their interaction, we found that its diffusion time increased at high FtsZ concentrations. The data presented in Supplementary Figure 2 shows ZapD in the presence of FtsZ-GDP (i.e. under non-polymerization conditions).

Was it possible to get a molecular weight estimate based on the diffusion time?

It is possible to estimate hydrodynamic volumes using the Stokes-Einstein equation if the diffusion coefficient of the diffusing particles is known, assuming that the particles are small and spherical. A molecular weight can then be estimated using a standard density of 1.35 g/cm3 (Fisher et all. Protein science 2009 DOI: 10.1110/ps.04688204). This estimate is heavily dependent on the shape of the diffusing particle, as we assume that our protein of interest here is far from a spherical shape due to the interaction through the flexible linker, the hydrodynamic volumes are overestimated. This overestimation then leads to a further overestimation of the molecular weight. In addition, for a more accurate estimation of the sizes and thus molecular weights for proteins, a modified model of the Stokes-Einstein equation is required (Tyn and Gusek Biotechnology and Bioengineering DOI: 10/1002/bit.260350402), where additional information about the shape of the diffusing particle is estimated by measuring the radius of gyration of the particle. These calculations are complex and beyond the scope of our manuscript.

Supple Fig 4:

Does FtsZ GTPase activity (without ZapD) also vary with KCl concentrations? It will be useful to comment on this in Supplementary Figure 4.

Yes, it has been previously reported that moderate concentration of KCl is optimal for FtsZ GTPase activity. We added a comment to the caption.

Page 6, line 42: short filament segments arranged nearly 'parallel' to each other Since FtsZ filaments are polar, it is better to rephrase as 'parallel or antiparallel'.

Corrected.

Page 7, line 41: cross linking of short 'FtsZ' filaments and not ZapD?

It was a typo. Corrected

Page 8: delete 'from above' in the title?

Corrected

The use of the phrases such as 'cross linking from the top'; 'binds to FtsZ from above' is vague. (Figure 5b legend; discussion page 10, line 18; page 8, line 26; page 12, line 27). Similarly labelling on a schematic figure on the use of vertical, diagonal/lateral will be useful for the readers.

We thank the reviewer for the suggestions to improve the understanding of our data. We have simplified them by renaming these interactions as vertical.

Page 13, lines 6 -10

Rather than an orientation of top or from the side, just the presence of multiple crosslinks along coaxial filaments suffices for a straight bundle. The average spacing will be more uniform in such a straight bundle compared to a toroid where there might be regions without ZapD. I do not find the data on an upward orientation convincing. ZapD binding need not be above to have the C-terminal ends of FtsZ pointing towards the membrane. On the other hand, having ZapD bind above is likely to occlude membrane binding of FtsZ?

The flexibility of the FtsZ linker suggests that ZapD can bind filaments oriented in different directions. In a cellular environment, FtsZ molecules interact with other division proteins that compete with ZapD for binding sites. This competition could prevent the membrane from occluding and instead create binding sites between the filaments, stabilizing them.

Page 11, lines 32 - 34: Please rephrase the sentence, with focus on the main point to be conveyed. Do the authors want to say that the 'Same molecule contributes to variability in spacing based on the number of connections formed.'

Thank you for your comment. We have rephrased the sentence for clarity.

Page 11: paragraphs 1,2, and 3 appears to convey similar, related ideas and are redundant. Could these be shortened further into one paragraph highlighting how the ratio leads to differences in higher order FtsZ organisation?

These paragraphs discuss different ideas, and it is better to keep them separate.

In the response to reviewers, page 19, point 5 (iii), it is given that 5000 FtsZ molecules correspond to 2/3rd of the total, while in the manuscript text, it is given as one-third. Please correct the response text/manuscript text accordingly. The numbers in the cited reference appears to suggest 1/3rd.

Yes, it was 1/3rd. Thanks for pointing that out. 

Fig 1b. Y-axis: Absorbance spelling has a typo.

Page 14, line 11: Healthcare ('h' missing)

Page 14, line 15: HCl, KCl (L should be in small letter)

Page15, line 18: 43 - 48K rpm (not Krpm)

Supple Fig 1 legend: line 5: 's' missing for species

Corrected.

DOI: 10.1038/s41598-019-43327-y

Resource: RRID:BDSC_51710

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_51710

What is this?

DOI: 10.1038/s41598-019-43327-y

Resource: RRID:BDSC_7198

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_7198

What is this?

DOI: 10.1007/s11055-025-01814-y

Resource: ImageJ (RRID:SCR_003070)

Curator: @scibot

SciCrunch record: RRID:SCR_003070

What is this?

DOI: 10.1002/mc.70018

Resource: RRID:Addgene_212936

Curator: @scibot

SciCrunch record: RRID:Addgene_212936

What is this?

Variables

Author response:

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

Reviewer #1 (Public review):

Summary:

The authors aim to explore the effects of the electrogenic sodium-potassium pump (Na⁺/K⁺-ATPase) on the computational properties of highly active spiking neurons, using the weakly-electric fish electrocyte as a model system. Their work highlights how the pump's electrogenicity, while essential for maintaining ionic gradients, introduces challenges in neuronal firing stability and signal processing, especially in cells that fire at high rates. The study identifies compensatory mechanisms that cells might use to counteract these effects, and speculates on the role of voltage dependence in the pump's behavior, suggesting that Na⁺/K⁺-ATPase could be a factor in neuronal dysfunctions and diseases

Strengths:

(1) The study explores a less-examined aspect of neural dynamics-the effects of (Na⁺/K⁺-ATPase) electrogenicity. It offers a new perspective by highlighting the pump's role not only in ion homeostasis but also in its potential influence on neural computation.

(2) The mathematical modeling used is a significant strength, providing a clear and controlled framework to explore the effects of the Na+/K+-ATPase on spiking cells. This approach allows for the systematic testing of different conditions and behaviors that might be difficult to observe directly in biological experiments.

(3) The study proposes several interesting compensatory mechanisms, such as sodium leak channels and extracellular potassium buffering, which provide useful theoretical frameworks for understanding how neurons maintain firing rate control despite the pump's effects.

Weaknesses:

(1) While the modeling approach provides valuable insights, the lack of experimental data to validate the model's predictions weakens the overall conclusions.

(2) The proposed compensatory mechanisms are discussed primarily in theoretical terms without providing quantitative estimates of their impact on the neuron's metabolic cost or other physiological parameters.

We thank the reviewer for their concise and accurate summary and appreciate the constructive feedback on the article’s strengths and weaknesses. Experimental work is beyond the scope of our modeling-based study. However, we would like our work to serve as a framework for future experimental studies into the role of the electrogenic pump current (and its possible compensatory currents) in disease, and its role in evolution of highly specialized excitable cells (such as electrocytes).

Quantitative estimates of metabolic costs in this study are limited to the ATP that is required to fuel the pump. By integrating the net pump current over time and dividing by one elemental charge, one can find the rate of ATP that is consumed by the Na⁺/K⁺pump for either compensatory mechanism. The difference in net pump current is thus proportional to ATP consumption, which allows for a direct comparison of the cost efficiency of the Na⁺/K⁺ pump for each proposed compensatory mechanism. The Na⁺/K⁺ pump is, however, not the only ATP-consuming element in the electrocyte, and some of the compensatory mechanisms induce other costs related to cell

‘housekeeping’ or presynaptic processes. We now added a section in the appendix titled

‘Considerations on metabolic costs of compensatory mechanisms’ (section 11.4), where we provide ballpark estimates for the influence of the compensatory mechanisms on the total metabolic costs of the cell and membrane space occupation. Although we argue that according these estimates, the impact of discussed compensatory mechanisms could be significant, due to the absence of more detailed experimental quantification, a plausible quantitative cost approximation on the whole cell level remains beyond the scope of this article.

Reviewer #1 (Recommendations for the authors):

(1)  For the f-I curves in Figures 1 and 6, the firing rate increases as the input current increases. I am curious to know: (a) whether the amplitudes of the action potentials (APs) vary with increased input current; (b) whether the waveform of APs (such as in Fig. 1I) transitions into smaller amplitude oscillations at higher input currents; and (c) if the waveform does change at higher input currents, how do the "current contributions," "current," and "ion exchanges per action potential" in Figures 1HJ and 6AB respond?

To fully answer these questions, we added a supplemental figure with accompanied text in section 11.1 (Fig. A1). We also added a reference to this figure in the main text (section 4.1). Here, it is shown that, as previously illustrated in [1], AP amplitude decreases when the input current increases (Fig. A1 A, left). This effect remains upon addition of either a pump with constant pump rate and co-expressed sodium leak channels (Fig. A1 A, center), or a voltage-dependent pump (Fig. A1 A, right). Interestingly, even though the shape of the current contributions (Fig. A1 B) and the APs (Fig. A1 C) look very different for low (Fig. A1 C, top) and high inputs (Fig. A1 C, bottom), the total sodium and potassium displacement per AP, and thus the pump rate, is roughly the same (Fig. A1 D). Under the assumption that voltage-gated sodium channel (NaV) expression is adjusted to facilitate fixed-AP amplitudes, however, (as in [1]) more NaV channels would be expressed in fish with higher synaptic drives. This would then result in an additional sodium influx per AP and result in higher energetic requirements per AP for electrocytes with higher firing rates (also shown in [1]).

(2) Could the authors clarify what the vertical dashed line represents in Figures 1B and 1F? Does it correspond to an input current of 0.63uA?

(Reviewer comment refers to Fig. 1C and 1F in new version): Yes, it corresponds to the input current that is also used in figures 1D and 1G. We clarified this by adding an additional tick label on the x-axis in 1F. The current input of 0.63uA was chosen as a representative input for this cell as follows: we first modeled an electrocyte with a periodic synaptic drive as in [1]. The frequency of this drive was set to 400 Hz, which is an intermediate value in the range of reported EODfs (and thus presumably pacemaker firing rates) of 200-600Hz [2]. Then, acetylcholine receptor currents I_AChRNa and I_AChRNa were summed and averaged to obtain the average input current of 0.63uA. This is now also explained in new Methods section 6.2.1.

(3) What input current was used for Figures 1H, 1I, and 1J?

Response: In a physiological setting, where the electrocyte is electrochemically coupled to the pacemaker nucleus, stimulation of the electrocyte occurs through neurotransmitter release in the synaptic cleft, which then leads to the opening of acetylcholine receptor channels. As figures 1H-J concern different ion fluxes, we aimed to also include currents stemming from acetylcholine receptor channels. We therefore did not stimulate the electrocyte with a constant input current as in Fig. 1C and F, but simulated elevated constant neurotransmitter levels in the synaptic cleft, which then leads to elevated acetylcholine receptor currents. In the model, this neurotransmitter level, or ‘synaptic drive’ is represented by parameter syn_clamp. A physiologically relevant value for syn_clamp was deduced by averaging the synaptic drive during a 400 Hz pacemaker stimulus. This is now also explained in new Methods section 6.2.1.

(4) In Figure 4A, there is a slight delay between the PN spikes (driver) and the EO (receiver), and no EO spikes occur without PN spikes. However, the firing rate of EO (receiver) appears to decrease before the chirp initiations in Fig 4B; and this delay seems to disappear in Fig 4C. Could the authors explain these observations?

As shown in the bottom right of figure 4A, when plotting the instantaneous firing rate as one over the inter-spike-interval (1/ISI), the firing rate of a cell is only plotted at the end of every ISI. Therefore, even though the PN drives the electrocyte and thus spikes earlier in time than the electrocyte, when it initiates chirps, these will only be plotted as an instantaneous firing rate at the end of the chirp. If the electrocyte fires spontaneously within this chirp, its instantaneous firing rate will appear earlier in time than the initiation of the chirp of the PN. The PN did, however, initiate the chirp before that and causality between the PN and electrocyte is not disturbed.

(5) Regarding Figure 6, could the authors specify the input current used in Figures 6A and 6B?

Figure 6A and 6B have the same synaptic drive as Fig. 1 H, I and J (syn_clamp=0.13).

(6) In Section 6, I would recommend that the authors provide a table of parameters and their corresponding values for clarity.

Thank you for your suggestion. We now reorganized the method section and added two tables with parameters for clarity. Table 1 (see Methods 6.1) includes all parameters that differ from the parameters reported in [1], and parameters that arise from the additionally modeled equations to simulate ion concentration dynamics and pump. We also added the parameters used to simulate the different stimulus protocols (and corresponding tuned parameters) that are presented in the article in Table 2 (see Methods 6.2).

Reviewer #2 (Public review):

Summary:

The paper 'The electrogenicity of the Na⁺/K⁺-ATPase poses challenges for computation in highly active spiking cells' by Weerdmeester, Schleimer, and Schreiber uses computational models to present the biological constraints under which electrocytes-specialized highly active cells that facilitate electro-sensing in weakly electric fish-may operate. The authors suggest potential solutions these cells could employ to circumvent these constraints.

Electrocytes are highly active or spiking (greater than 300Hz) for sustained periods (for minutes to hours), and such activity is possible due to an influx of sodium and efflux of potassium ions into these cells for each spike. This ion imbalance must be restored after each spike, which in electrocytes, as with many other biological cells, is facilitated by the Na-K pumps at the expense of biological energy, i.e., ATP molecules. For each ATP molecule the pump uses, three positively charged sodium ions from the intracellular space are exchanged for two positively charged potassium ions from the extracellular volume. This creates a net efflux of positive ions into the extracellular space, resulting in hyperpolarized potentials for the cell over time. This does not pose an issue in most cells since the firing rate is much slower, and other compensatory mechanisms and other pumps can effectively restore the ion imbalances. In electrocytes of weakly electric fish, however, that operate under very different circumstances, the firing rate is exceptionally high. On top of this, these cells are also involved in critical communication and survival behaviors, emphasizing their reliable functioning.

In a computation model, the authors test four increasingly complex solutions to the problem of counteracting the hyperpolarized states that occur due to continuous NaK pump action to sustain baseline activity. First, they propose a solution for a well-matched Na leak channel that operates in conjunction with the NaK pump, counteracting the hyperpolarizing states naturally. Additionally, their model shows that when such an orchestrated Na leak current is not included, quick changes in the firing rates could have unexpected side effects. Secondly, they study the implication of this cell in the context of chirps - a means of communication between individual fishes. Here, an upstream pacemaking neuron entrains the electrocyte to spike, which ceases to produce a so-called chirp - a brief pause in the sustained activity of the electrocytes. In their model, the authors show that it is necessary to include the extracellular potassium buffer to have a reliable chirp signal. Thirdly, they tested another means of communication in which there was a sudden increase in the firing rate of the electrocyte followed by a decay to the baseline. For reliable occurrence of this, they emphasize that a strong synaptic connection between the pacemaker neuron and the electrocyte is warranted. Finally, since these cells are energy-intensive, they hypothesize that electrocytes may have energyefficient action potentials, for which their NaK pumps may be sensitive to the membrane voltages and perform course correction rapidly.

Strengths:

The authors extend an existing electrocyte model (Joos et al., 2018) based on the classical Hodgkin and Huxley conductance-based models of Na and K currents to include the dynamics of the NaK pump. The authors estimate the pump's properties based on reasonable assumptions related to the leak potential. Their proposed solutions are valid and may be employed by weakly electric fish. The authors explore theoretical solutions that compound and suggest that all these solutions must be simultaneously active for the survival and behavior of the fish. This work provides a good starting point for exploring and testing in in vivo experiments which of these proposed solutions the fish use and their relative importance.

Weaknesses:

The modeling work makes assumptions and simplifications that should be listed explicitly. For example, it assumes only potassium ions constitute the leak current, which may not be true as other ions (chloride and calcium) may also cross the cell membrane. This implies that the leak channels' reversal potential may differ from that of potassium. Additionally, the spikes are composed of sodium and potassium currents only and no other ion type (no calcium). Further, these ion channels are static and do not undergo any post-translational modifications. For instance, a sodium-dependent potassium pump could fine-tune the potassium leak currents and modulate the spike amplitude (Markham et al., 2013).

This model considers only NaK pumps. In many cell types, several other ion pumps/exchangers/symporters are simultaneously present and actively participate in restoring the ion gradients. It may be true that only NaK pumps are expressed in the weakly electric fish Eigenmannia virescens. This limits the generalizability of the results to other cell types. While this does not invalidate the results of the present study, biological processes may find many other solutions to address the non-electroneutral nature of the NaK pump. For example, each spike could include a small calcium ion influx that could be buffered or extracted via a sodium-calcium exchanger.

Finally, including testable hypotheses for these computational models would strengthen this work.

We thank the reviewer for the detailed summary and the identified weaknesses according to which we improved our article. Our model assumptions and simplifications are now mentioned in more detail in the introduction of the article (section 3), and justified in the Methods (section 6.1).

Furthermore, we added a discussion section (section 5.1) where we outline the conditions under which the present study can be extended to other cell types. We now also state more clearly that the pump current will be present for any excitable cell with significant sodium flux (assuming that the NaK pump carries out the majority of its active transport), but that compensatory mechanisms (if employed at all in a particular cell) could also be implemented via other ionic currents and transporters. We furthermore now highlight the testable hypotheses that we put forward with our computational study on the weakly electric fish electrocyte more explicitly in the first paragraph of the discussion.

Reviewer #2 (Recommendations for the authors):

Main text

Please explicitly state this model's assumptions in the introduction and elaborate on them in the discussion if necessary. For example, some assumptions that I find relevant to mention are: - The Na and K channels are classic HH conductance-based channels, with no post-translational modifications or beta subunit modifications as seen in other high-frequency firing cells (10.1523/JNEUROSCI.23-12-04899.2003).

Neither calcium nor chloride ions are considered in the spike generation. Nor are Na-dependent K channels (10.1152/jn.00875.2012).

Only the Na-K pump (and not the Na-Ca exchanger, Ca-pump, or Cl pumps) is modeled,

Calmodulin, which can buffer calcium, is highly expressed in electric eels, but it is not considered. If some of these assumptions have valid justifications in weakly electric fish electrocytes, please state so with the citations. I recognize that including these in your models is beyond the scope of the current paper.

We thank the reviewer for pointing out this issue. We now specified in the introduction that the model only contains sodium and potassium ions and only classic HH conductance-based channels. We there also explicitly specify the details on the Na⁺/K⁺-ATPase: it is the only active transporter in this model, thus solely responsible for maintaining ionic homeostasis; its activity is only modulated by intracellular sodium and extracellular potassium concentrations. In the discussion (6.1), we now elaborate on how ion-channel-related aspects (i.e., the addition of resurgent Na⁺ or Na⁺ -dependent K⁺ channels), additional ion fluxes (including some not relevant for the electrocyte but for other excitable cells), and additional active transporters and pumps would influence the results presented in the article.

In addition, there might be other factors that the authors and the reviewers have yet to consider. The model is a specific case study about the weakly electric fish electrocyte with high-frequency firing. It is almost guaranteed that biology will find other compensatory ways in different cell types, systems, and species (auditory nerve, for example). Given this, it would be prudent to use phrases such as 'this model suggests,' 'perhaps,' 'could,' 'may,' and 'eludes to,' etc., to accommodate other possible solutions to ion homeostasis in rapidly spiking neurons. The solutions the authors are proposing are some of many.

We rephrased some of the statements to highlight more the hypothetical nature of the compensatory mechanisms in specific cells and to draw attention to the fact that there can be many more such factors. This fact is now also explicitly mentioned in discussion section 5.2.

Figures

Some of my comments on the figures are stylistic, others are to improve clarity, and some are critical for accuracy.

The research problem concerns weakly electric fish E. virescens. I suggest introducing a picture of an electric fish in the beginning (such as that in Figure 3, but not exactly; see specific comments on this fish figure) along with a schema of the research question. 

We agree, and added an overview schema in Fig. 1A.

Font sizes change between the panels in all the figures. Please maintain consistency. The figure panel titles and axis labels should start with a capital letter.

Thank you for pointing this out, both issues have been resolved in the new version of the article.

Figure 1:

Please rearrange the figure - BCFG belong together and should appear in the same order. The x-axis labels could be better placed.

Consider using fewer pump current f-I curves (B, D, E, F). Five is sufficient to make the point. Having 10 curves adds to the clutter. The placement of the color bar could be better. Similarly, the placement of the panel titles 'without co-expression' and 'with co-expression' and the panel labeling (BCFG) makes it confusing. The panel labels should be above the panel title.

Response (C, D, F, G in new version): We improved the layout of figure 1. Panels B, C, F, G are now C, D, F, G. We opted to include panel E before panels F and G, because it shows the coexpression mechanism before its effect on the tuning curve. We did move the colorbar, added x-axis labels to B and C, and adjusted the location of the panel labels for clarity. We also plotted fewer pump currents.

B, F: What does the dashed line indicate?

Response (C, F in new version): The dashed line indicates the input current that was used in figures 1D and 1G. We now clarified this by adding this value on the x-axis.

C: Any reason not to show the lower firing rates?

Response (B in new version): In the previous version of the article, pump currents were estimated for electrocytes that were stimulated with the mean synaptic drive that stems from periodic stimulation in the 200-600 Hz regime. We now extended the range of synaptic inputs to obtain lower (and higher) firing rates. The linear relationship between firing rate and pump current also holds for these additional firing rates.

D: There is no difference between the curves at the top and the bottom. One fills the area between the curve and the zero line; the other shows the curve itself. Please use only one of the two representations.

Response (panel I in new version): In the previous version, the difference between the plots was that one showed the absolute values of the currents (the curves), and the other plot showed the contributions of the currents to the total (area between the curves). We now only depict the current contributions.

The I and H orders can be swapped.

Thank you, they are now swapped.

The colors used for Na and K are very dull (light blue and pink).

We now use darker colors in the new version of the article.

Figure 2:

Please verify that without the synaptic input perturbations (i.e., baseline in A, D), the firing rate (B, E) and pump current (C, F) converge to the baseline. There is a noticeable drift (downward for firing rate and upward for pump currents) at the 10-second time point.

Thanks to you noticing, we identified a version mismatch in the code that estimates the pump current required for ionic homeostasis (see Methods 6.1.2). We have now corrected the code and made sure to start the simulation in the steady state so that there is no drift at baseline firing. We also used this corrected code to present tuned parameters for different stimulus protocols in Table 2 (Methods 6.2).

Figure 3:

A. The dipole orientation with respect to the fish in panel B needs to be corrected. Consider removing this as this work is not about the dipole.

This panel has been removed.

B. This figure has already been overused in multiple papers; please redraw it. Localized expressions of different pumps and ion channels are present within each electrocyte, which generates the dipole. Either show this correctly or don't at all (the subfigure pointed out by the red arrow).

This panel has been moved to Fig. 1A. We opted to remove the localized expressions.

C and D belong together; please place them next to each other. Consider introducing panel D first since it follows a similar protocol to the last figure.

Response (A in new version): Panel placement has been adjusted. We opted to maintain the order to maintain the flow of the text, but we do now combine them in one panel.

E and F are very similar in that they are swapped on the x and y axes. Either that or I have severely misunderstood something, in which case it needs to be shown better.

Response (B and C in new version): We adjusted the placement of these panels. They are not the same, panel B shows the mean of physiological periodic inputs, and figure C shows that when this mean is fed to the electrocyte, it also induces tonic firing. The range of mean currents that result from periodic synaptic stimulation in the physiological regime (panel B, y-axis) is now indicated in panel C by a grey box along the x-axis.

G. Why show the lines with double arrow ends? The curves are diverging - that's enough.

Good point, we updated this panel accordingly (now panel D).

Figure 4

Please verify the time units in these plots. Something seems amiss. B and D lower plots-perhaps this is seconds? B could use an inset box/ background gray color (t1, t2) indicating the plots of the C panel (left, right). Likewise, for D (t1, t2), connect to E (left, right).

You are right, the x-axes were supposed to be in seconds, we updated this. We indicated the relations between D-C and D-E by gray backgrounds and by adding the corresponding panel label on the x-axis.

A: Indicate the perturbation in the schematic, i.e., extracellular K buffer.

The perturbation is now indicated.

D: Even with the extracellular K buffer, there is a decay (slower than in B) of the pump current over time. Please verify (you do not have to show in your paper) that this decay saturates.

After the ten chirps are initiated, pacemaker firing goes back to baseline. In both cases (panel B and panel D), the pump current goes back to baseline after some time. With extracellular potassium buffering, this happens more slowly due to a decreased reaction speed of the pump to changes in firing rate (in comparison to the case without extracellular potassium buffer).

The decrease in reaction speed however merely delays the effects of changes in firing rates on the pump current in time. Therefore, even with an extracellular potassium buffer, when more chirps are initiated in a short period of time, the pump current can still decrease to an extent that impairs entrainment. Using the same protocol as in panel B and D, we increased the number of chirps and found that with an extracellular potassium buffer, a maximum of 13 chirps could be encoded without entrainment failure (as opposed to 2 chirps without the buffer as shown in panel B).

Figure 5

Please verify the time units in these plots, as for Figure 4. B and E lower plots-perhaps this is seconds? B could use an inset box/ background gray color (t1, t2) indicating the plots of the panels C and D. Likewise, for E (t1, t2), connect to F and G.

The time axis in this figure was indeed also in seconds, which we corrected here. The relations between plots B-C/D and E-F/G are now indicated through gray backgrounds and corresponding panel references on the x-axis.

A: Indicate the perturbation in the schematic, i.e., the synapse's strength. There is no need to include the arrow or to mention freq. rise. The placement of the time scale can be misinterpreted as a current clamp. Instead, plot it as a zoomed inset.

The arrow is removed and we now also show a zoomed inset. Also, the perturbation is now indicated.

E: Verify that the pump current in the strong synapse case already starts at 1.25

We verified this and noticed that the pump current in the strong synapse case is indeed lower than that in the weak synapse case. This is because to ensure a fair comparison for this stimulation protocol, voltage-gated sodium channel conductance was tuned to maintain a spike amplitude of 13 mV in both cases (see Methods 6.2). In this case, a weak synapse leads to a lower influx of sodium via AChR channels, but a higher influx via voltage-gated sodium channels. The total sodium influx in this case is larger than that for a stronger synapse with relatively less voltage-gated sodium currents, and thus a larger pump current. In the previous version of the article, this was wrongly commented on in the figure captions, and we removed the erroneous statement.

This is not critical, but because the R-value here can be obtained as a continuous value, it would be appropriate to show it for the whole duration of the weak and strong synapses in B and E. Maybe consider including a schema that shows how R is calculated in panel A.The caption has a typo, 'during frequency rises before (D) and after (E)'. It should be before C) and after (D) instead.

The caption typo has been corrected. The R-value for the whole duration of the weak and strong synapses in B and E is 1.000. This is because the R-value is the variance of all phase relations between the PN and the electrocyte, and for the entire duration of the stimulus protocol, there are only a few outliers in phase relations at the maxima of the frequency rises. We decided to include this R-value to show that in general, synchronization between the PN and the electrocyte is very stable. The schema that explains how R is calculated has not been included in favor of not overcrowding the figure. We did add a reference in the figure caption to the methods section in which the calculation of R is explained.

Figure 6:

A: The top and bottom plots are redundant. Use one of the two. They show the same thing. It may be better to plot Na, K, pump, and net currents on the top panels and the Na leak, which is of smaller magnitude, in a different panel.

We now only show current contributions.

B: Please change the color schema. It is barely visible on my prints.

D: Pump current, instantaneous case, is barely visible

Color schemes were adjusted.

Figure A1: It's all good.

Methods:

Please provide some internal citations for where specific equations were used in the results/figures. You do this for sections 6.2.3, referencing Figure 5 (c,d,e,g), and 6.2.4, referencing Fig 5 C-E.

There are now internal references in each methods section to where in the figures they were used. We also included a table with stimulus parameters for each figure with a stimulus protocol (Table 2).

Also, the methods could be ordered in the same order as the results are presented. Please consider if some details in the methods could be moved to the appendix.

The ordering of the methods has now been changed to separately explain the model expansions (6.1) and the stimulus protocols (6.2). Both sections are in corresponding order of the figures presented in the article. We opted to maintain all details in the methods.

6.1.1 Please cite 26 after the first line. Where was this used? In Figure 3C, 4, 5?

We added the citation. The effects of co-expressed leak channels are shown in Fig. 1 EG, and were used to compensate for pump currents at baseline firing in figures 1 D, H-J (left, with pump), 2, 4, 5, and 6 A-B (left), C (top). This is now also added to the text for clarity.

Traditionally (Hodgkin, A. L. and Huxley, A. F. (1952). J. Physiol. (Lond.), 117:500-544. Table 3; & Hodgkin, A. L. and Huxley, A. F. (1952). J. Physiol. (Lond.), 116:473-496 Table 5 and the paragraph around it), leak potential is set such that it accounts for all leak from all ions. While in your work, this potential is equal to the reversal of potassium - it need not be so in the animal. There may be leaks from other ions as well, particularly sodium and chloride. Please verify that assuming the leak reversal is the same as that of potassium (Ek, in Equation 3) does not lead to having to model Na leak currents separately.

In the original model [1], it was assumed that the reversal potential of the leak was the same as that of potassium, which contains the implicit assumption that only potassium ions contribute to the leak. In our article, we also assume that sodium ions contribute to the leak. This can be modeled by adjusting the leak reversal potential accordingly, or by adding an additional leak current that solely models the sodium leak. We opted for the latter in order to track all sodium and potassium ions separately so that ion concentration dynamics could also be modeled properly. Chloride ions were neglected in this study; in our model they do not contribute to the leak. If one were to also model chloride currents and chloride concentration dynamics, it would be beneficial to model these as an additional separate leak current.

The notation of I_pump_0 needs to be more convenient. Please consider another notation instead of the _0 (pump at baseline). Similarly for [Na⁺]_in_0 [Na⁺]_out_0 and [K⁺]_in_0 and [K+]_out_0

We changed the notation for baseline similarly to [3], with ‘0’ as a superscript instead of a subscript.

Equation 11: Please mention why AChRs do not let calcium ions through. Please cite a justification for this. If this is an assumption of the model, please state this explicitly.

The AChR channels that were found in the E. virescence electrocytes are muscle-type acetylcholine nicotinic receptors [4], which are non-selective cation channels that could indeed support calcium flux [5]. No calcium currents were, however, modeled in the original electrocyte model [1], presumably due to the lack of significant contributions of calcium currents or extracellular calcium concentrations to electrocyte action potentials of a similar weakly electric electrogenic wave-type fish Sternopygus macrurus [6].

Due to the lack of calcium currents in the original electrocyte model, and due to the limitation of this study to sodium and potassium ions, we chose not to include calcium currents stemming from AChR channels. This assumption is now explicitly stated in Methods 6.1.

Equation 12, V_in, where the intracellular volume. If possible, avoid the notation of 'V' - you already use a small v for membrane potential.

We changed the notation for volume to ‘ω’ similarly to [3]. As we previously used ω as a notation for the firing rate, we changed the notation for firing rate to ‘r’.

Equation 17: Does this have any assumptions? Would the I_AchRNa, and thus Sum(mean(I_Na))) not change depending on the synaptic drive?

The assumptions of this equations are the following (now also mentioned in Methods 6.1.2):

The sum of all sodium currents also includes sodium currents through acetylcholine channels (I_AChRNa).

All active sodium transport (from intra- to extracellular space) is carried out by the Na⁺/K⁺-ATPase, and active sodium transport through additional transporters and pumps is negligible.

The time-average of sodium currents is either taken in a tonic firing regime where the timeinterval that is averaged over is a multiple of the spiking period, nT, or if it is taken for a more variable firing regime, the size of the averaging window should be sufficiently large to properly sample all firing statistics.

Under these assumptions, Eq. 17 can be used to compute suitable pump currents for different synaptic drives (as Sum(mean(I_Na))) and thus I_pump0 indeed change with the synaptic drive, see Table 2 in Methods 6.2). 

6.2: Please rewrite the first sentence of this paragraph.

The first sentence of this paragraph, which has been moved to section 6.2.2 for improved structuring of the text, has been rewritten.

6.2.1: The text section could use a rewrite.

Please elaborate on what t_p is. If it is not time, please do not use 't.' What is p here? What are the units of the equation (22), t_p < 0.05 (?)

This section has now also been moved to 6.2.2. It has been rewritten to improve clarity and t_p has been renamed to t_pn (as it does reflect time, which is now better explained). The units have now also been added to the equation (which is now Eq. 26).

6.2.4: Please rewrite this.

This section has been rewritten (and has been moved to section 6.1.4).

Bibliography

Some references are omitted (left anonymous) or inconsistent on multiple occasions.

Thank you for pointing this out! It is now rectified.

References used for author response

(1) Joos B, Markham MR, Lewis JE, Morris CE. A model for studying the energetics of sustained high frequency firing. PLOS ONE. 2018 Apr;13:e0196508.

(2) Hopkins CD. Electric communication: Functions in the social behavior of eigenmannia virescens. Behaviour. 1974;50(3-4):270–304.

(3) Hübel N, Dahlem MA. Dynamics from seconds to hours in hodgkin-huxley model with time-dependent ion concentrations and buer reservoirs. PLoS computational biology.ff2014;10(12):e1003941.

(4) BanY, Smith BE, Markham MR. A highly polarized excitable cell separates sodium channels from sodium-activated potassium channels by more than a millimeter. Journal of neurophysiology. 2015; 114(1):520–30.

(5) Vernino S, Rogers M, Radcliffe KA, Dani JA. Quantitative measurement of calcium flux through muscle and neuronal nicotinic acetylcholine receptors. Journal of Neuroscience. 1994;14(9):5514-5524.

(6) Ferrari M, Zakon H. Conductances contributing to the action potential of sternopygus electro-cytes. Journal of Comparative Physiology A. 1993;173:281–92.

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_34340

Curator: @maulamb

SciCrunch record: RRID:BDSC_34340

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_36303

Curator: @maulamb

SciCrunch record: RRID:BDSC_36303

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_44034

Curator: @scibot

SciCrunch record: RRID:BDSC_44034

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_26198

Curator: @scibot

SciCrunch record: RRID:BDSC_26198

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_11749

Curator: @scibot

SciCrunch record: RRID:BDSC_11749

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_5824

Curator: @scibot

SciCrunch record: RRID:BDSC_5824

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_91367

Curator: @scibot

SciCrunch record: RRID:BDSC_91367

What is this?

DOI: 10.1038/s41467-025-62043-y

Resource: RRID:BDSC_9306

Curator: @scibot

SciCrunch record: RRID:BDSC_9306

What is this?

DOI: 10.1038/s41598-025-10166-z

Resource: (RRID:CVCL_D615)

Curator: @scibot

SciCrunch record: RRID:CVCL_D615

What is this?

DOI: 10.1038/s41467-025-61697-y

Resource: RRID:BDSC_5905

Curator: @scibot

SciCrunch record: RRID:BDSC_5905

What is this?

Prueba de anotación

Author response:

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

Reviewer #1 (Recommendations for the authors):

Line 122: There were a number of qualitative descriptors in the paper. For instance, if the authors want to say massive campaign, how massive? How rapid? These are relative terms in this context.

We have revised the text to minimize qualitative descriptors and to provide concrete numbers where possible. The revised sentence (line 121) now reads “We began our structural investigation of nitrogenase evolutionary history by conducting on a large-scale structure prediction analysis of 5378 protein structures, a more than threefold increase compared to available nitrogenase structures in the PDB. We then analyzed our phylogenetic dataset to identify notable structural changes.”

Line 179: "massively scale up" How massive?

We agree with the reviewer’s observation, in response, we have removed the phrase “massively scale up” and revised the text.

Line 182: "no compromise on alignment depth and negligible cost to prediction accuracy". How do you know this? Is this shown somewhere? Was there a comparison between known structures and the predicted structure for those nitrogenases that have structures?

In response to this comment, we have made several clarifications and revisions in the manuscript:

We modified Figure S1, which now shows the pLDDT (per-residue confidence metric from Alphafold) values of all our predictions. These scores are consistently high (over 90 for the D and K subunits, and approximetly 90 for the H subunits) regardless of whether the recycling protocol or the bona-fide protocol was used.

The reviewer’s comment demonstrated to us that the Figure S1 needed to more clearly representing these values, we therefore updated it accordingly.

To prevent any misinterpretation of our claims about the accuracy and cost of the method , we have revised the text at line 179, as follows:

“In total, 2,689 unique extant and ancestral nitrogenase variants were targeted. All structures were generated in approximately 805 hours, including GPU computations and MMseqs2 alignments performed using two different protocols: one for extant or most likely ancestral sequences, and another for ancestral variants.”

To support our analyses further, Figure S10A compares our model predictions with available PDB structures for nitrogenases.

Additionally, Figure S10B compare our predicted structures with the experimental structures reported in this article. In all cases, we observe low RMSD values.

Line 220: "fall within 2 angstroms" instead of "fall 2A"?

We have updated it in the text.

Line 315: It is not clear how the binding affinities and other measurements in Figure 4 and S6C were measured, and it is not discussed in the material and methods.

We thank the reviewer for pointing out this lack of clarity. The binding affinity estimations were performed using Prodigy. We have updated the main text (see line 322) to explicitly state that binding affinities were estimated using Prodigy. In addition, we have expanded the Materials and Methods section to include additional information about the structure characterization methods (lines 745-749). Previously, these details were only noted in Supplementary Table S6.

Line 510-511: "Subtle, modular structural adjustments away from the active site were key to the evolution and persistence of nitrogenases over geologic time". This seems like a bit of an overstatement. While the authors see structural differences in the ancestral nitrogenase and speculate these differences could be involved in oxygen protection, there is no evidence that the ancestral nitrogenase is more sensitive to oxygen than the extant nitrogenase.

We appreciate the reviewer’s comment. Our intention was to emphasize that subtle, modular structural adjustments might have contributed to oxygen protection rather than to assert that ancestral nitrogenases are more oxygen-sensitive than their extant counterparts. We have revised the text to clarify.

Reviewer #2 (Recommendations for the authors):

What is the reference for the measured RMSDs in Fig 2A? What is the value on the y-axis? The range of 'Count' is unclear, given that there are 5000 structures predicted in the study.

Figure 2A presents a histogram of RMSD values from all pairwise alignments among 769 structures (385 extant and 384 ancestral DDKK), totaling 591,361 comparisons. We excluded ancestral DDKK variants due to computational limitations.  

Similarly, what is the sequence identity in Figure 2B calculated relative to?

In Figure 2B, sequence identities are derived from pairwise comparisons across all structures in our dataset. Each value represents the identity between two specific structures, rather than being measured against a single reference.

The claim that 'structural analysis could reproduce sequence-based phylogenetic variation' should probably be tempered or qualified, given that the RMSD differences calculated are so low.

We hope to have addressed the concerns about the low RMSD values in the previous comments. We have revised the text (line 204), which now reads: “it still strongly correlates with sequence identity (Figure 2B), indicating that even minor structural variations can recapitulate sequence-based phylogenetic distinctions.”

How are binding affinities (Figure 4) calculated?

We have now clarified the binding affinity calculations in the main text. The model used is now detailed at line 322, with additional information provided in the Methods section.

Presumably, crystallized proteins (Anc1A, Anc1B, Anc2) were also among those whose structures were predicted with AF. A comparison should be provided of the predicted and crystallized structures, as this is an excellent opportunity to further comment on the reliability of AlphaFold.

In the revised manuscript, Figure S10 now present structural comparisons between the crystallized proteins and their AlphaFold-predicted counterparts.

The labels in Figure 5B are not clear. Are the 3rd and 4th panels also comparative RMSD values? But only one complex name is provided.

We appreciate this feedback and now revised the Figure 5B for clarity.

Page 9 line 220, missing word: 'varaints fall within/under 2angstroms'

We thank the reviewer for the correction, we have updated the text.

Este concepto subraya la importancia de la claridad en la definición del problema de investigación. Al formular el problema con precisión, el investigador establece una base sólida que facilita la identificación de objetivos, preguntas y métodos, reduciendo ambigüedades y enfocando el estudio hacia resultados concretos. Esto refuerza la necesidad de dedicar tiempo a la revisión de literatura y al análisis del contexto para garantizar que el problema sea comprensible y relevante.

El título actúa como una "tarjeta de presentación" del proyecto, condensando la esencia de la investigación. Incluir especificidad, espacialidad y temporalidad asegura que el título sea claro y delimite el alcance del estudio. Por ejemplo, un título como "Conocimientos sobre COVID-19 en estudiantes de la UVG, 2022" define claramente qué, dónde y cuándo, ayudando a los lectores a comprender inmediatamente el enfoque y contexto del trabajo.

Este punto destaca la importancia de contextualizar el problema para darle relevancia. Describir quiénes se ven afectados y las implicaciones (causas y consecuencias) permite al investigador justificar la pertinencia del estudio y conectar con las necesidades reales de una población o situación. Esto refuerza que un buen enunciado no solo describe el problema, sino que lo sitúa en un marco social, cultural o práctico significativo.

Este principio resalta la importancia de la reproducibilidad en la investigación científica. Detallar los recursos materiales (como software, equipos o documentos) asegura que el estudio sea transparente y verificable. Este aprendizaje refuerza que una investigación bien planificada considera no solo la ejecución, sino también la posibilidad de que otros puedan replicarla para validar los resultados.

Usar un formato estandarizado como APA asegura que las fuentes sean citadas de manera clara y profesional, facilitando la trazabilidad de la información. Este aprendizaje refuerza la importancia de la integridad académica, ya que citar correctamente no solo da crédito a los autores originales, sino que también permite a otros investigadores acceder a las fuentes para profundizar en el tema.

El modelo SMART es una herramienta clave para garantizar que los objetivos sean prácticos y efectivos. Por ejemplo, un objetivo como "Demostrar los conocimientos de los estudiantes sobre Check4Covid en 2022" es específico (conocimientos), medible (a través de encuestas), alcanzable (dentro del contexto de la UVG), relevante (para la prevención de COVID-19) y temporal (en 2022). Este enfoque refuerza la importancia de diseñar objetivos que guíen la investigación sin desviarse.

Las preguntas auxiliares son esenciales para desglosar el problema en aspectos manejables. Esta pregunta específica guía la investigación hacia la evaluación de la efectividad de una herramienta, promoviendo un análisis crítico de sus fortalezas y limitaciones. Aprender a formular preguntas claras y enfocadas, como esta, ayuda a estructurar la investigación y a mantener el rumbo hacia el objetivo general.

La justificación es el "corazón" persuasivo de la investigación, ya que conecta el problema con su relevancia práctica o teórica. Al explicar por qué el estudio es necesario, el investigador no solo motiva su realización, sino que también convence a otros (como financiadores o académicos) de su valor. Este aprendizaje enfatiza la necesidad de alinear el proyecto con necesidades reales o vacíos de conocimiento.

Este concepto subraya la importancia de la claridad en la definición del problema de investigación. Al formular el problema con precisión, el investigador establece una base sólida que facilita la identificación de objetivos, preguntas y métodos, reduciendo ambigüedades y enfocando el estudio hacia resultados concretos. Esto refuerza la necesidad de dedicar tiempo a la revisión de literatura y al análisis del contexto para garantizar que el problema sea comprensible y relevante.

El título actúa como una "tarjeta de presentación" del proyecto, condensando la esencia de la investigación. Incluir especificidad, espacialidad y temporalidad asegura que el título sea claro y delimite el alcance del estudio. Por ejemplo, un título como "Conocimientos sobre COVID-19 en estudiantes de la UVG, 2022" define claramente qué, dónde y cuándo, ayudando a los lectores a comprender inmediatamente el enfoque y contexto del trabajo.

El modelo SMART es una herramienta clave para garantizar que los objetivos sean prácticos y efectivos. Por ejemplo, un objetivo como "Demostrar los conocimientos de los estudiantes sobre Check4Covid en 2022" es específico (conocimientos), medible (a través de encuestas), alcanzable (dentro del contexto de la UVG), relevante (para la prevención de COVID-19) y temporal (en 2022). Este enfoque refuerza la importancia de diseñar objetivos que guíen la investigación sin desviarse.

Las preguntas auxiliares son esenciales para desglosar el problema en aspectos manejables. Esta pregunta específica guía la investigación hacia la evaluación de la efectividad de una herramienta, promoviendo un análisis crítico de sus fortalezas y limitaciones. Aprender a formular preguntas claras y enfocadas, como esta, ayuda a estructurar la investigación y a mantener el rumbo hacia el objetivo general.

Este punto destaca la importancia de contextualizar el problema para darle relevancia. Describir quiénes se ven afectados y las implicaciones (causas y consecuencias) permite al investigador justificar la pertinencia del estudio y conectar con las necesidades reales de una población o situación. Esto refuerza que un buen enunciado no solo describe el problema, sino que lo sitúa en un marco social, cultural o práctico significativo.

Este principio resalta la importancia de la reproducibilidad en la investigación científica. Detallar los recursos materiales (como software, equipos o documentos) asegura que el estudio sea transparente y verificable. Este aprendizaje refuerza que una investigación bien planificada considera no solo la ejecución, sino también la posibilidad de que otros puedan replicarla para validar los resultados.

La justificación es el "corazón" persuasivo de la investigación, ya que conecta el problema con su relevancia práctica o teórica. Al explicar por qué el estudio es necesario, el investigador no solo motiva su realización, sino que también convence a otros (como financiadores o académicos) de su valor. Este aprendizaje enfatiza la necesidad de alinear el proyecto con necesidades reales o vacíos de conocimiento.

Usar un formato estandarizado como APA asegura que las fuentes sean citadas de manera clara y profesional, facilitando la trazabilidad de la información. Este aprendizaje refuerza la importancia de la integridad académica, ya que citar correctamente no solo da crédito a los autores originales, sino que también permite a otros investigadores acceder a las fuentes para profundizar en el tema.

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

1. Response to reviewers

We would like to thank the reviewers for carefully reading our manuscript and for their valuable comments in support for the publication of our investigation of rapid promoter evolution of accessory gland genes between Drosophila species and hybrids. We are glad to read that the reviewers find our work interesting and that it provides valuable insights into the regulation and divergence of genes through their promoters. We are encouraged by their acknowledgement of the overall quality of the work and the importance of our analyses in advancing the understanding of cis-regulatory changes in species divergence.

2. Point-by-point description of the revisions

Reviewer #

Reviewer Comment

Author Response/Revision

Reviewer 1

The authors test the hypothesis that promoters of genes involved in insect accessory glands evolved more rapidly than other genes in the genome. They test this using a number of computational and experimental approaches, looking at different species within the Drosophila melanogaster complex. The authors find an increased amount of sequence divergence in promoters of accessory gland proteins. They show that the expression levels of these proteins are more variable among species than randomly selected proteins. Finally, they show that within interspecific hybrids, each copy of the gene maintains its species-specific expression level.

We thank Reviewer 1 for their detailed review and positive feedback on our manuscript, and for their helpful suggestions. We have now fully addressed the points raised by Reviewer 1 and have provided the suggested clarifications and revisions to improve the flow, readability, and presentation of the data, which we believe have improved the manuscript significantly.

The work is done with expected standards of controls and analyses. The claims are supported by the analysis. My main criticism of the manuscript has to do not with the experiments or conclusion themselves but with the presentation. The manuscript is just not very well written, and following the logic of the arguments and results is challenging.

The problem begins with the Abstract, which is representative of the general problems with the manuscript. The Abstract begins with general statements about the evolution of seminal fluid proteins, but then jumps to accessory glands and hybrids, without clarifying what taxon is being studied, and what hybrids they are talking about. Then, the acronym Acp is introduced without explanation. The last two sentences of the Abstract are very cumbersome and one has to reread them to understand how they link to the beginning of the Abstract.

More generally, if this reviewer is to be seen as an "average reader" of the paper, I really struggled through reading it, and did not understand many of the arguments or rationale until the second read-through, after I had already read the bottom line. The paragraph spanning lines 71-83 is another case in point. It is composed of a series of very strongly worded sentences, almost all starting with a modifier (unexpectedly, interestingly, moreover), and supported by citations, but the logical flow doesn't work. Again, reading the paragraph after I knew where the paper was going was clearer, but on a first read, it was just a list of disjointed statements.

Since most of the citations are from the authors' own work, I suspect they are assuming too much prior understanding on the part of the reader. I am sure that if the authors read through the manuscript again, trying to look through the eyes of an external reader, they will easily be able to improve the flow and readability of the text.

We thank the reviewer for their detailed feedback and are glad that they acknowledge our work fully supports the claims of our manuscript. We also appreciate their helpful suggestions for improving the readability of the manuscript and have done our best to re-write the abstract and main text where indicated. In particular, the paragraph between lines 71-83 have been rewritten and we have taken care to write to non-expert readers.

1) In the analysis of expression level differences, it is not clear what specific stage / tissue the levels taken from the literature refer to. Could it be that the source of the data is from a stage or tissue where seminar fluid proteins will be expressed with higher variability in general (not just inter-specifically) and this could be skewing the results? Please add more information on the original source of the data and provide support for their validity for this type of comparison.

These were taken from publicly available adult male Drosophila datasets, listed in the data availability statement and throughout the manuscript. We have provided more detail on the tissue used for analysis of Acp gene expression levels.

2) The sentence spanning lines 155-157 needs more context.

We have added more context to lines 155-157.

3) Line 203-204: What are multi-choice enhancers?

We replaced the sentence with "... such as rapidly evolving enhancers or nested epistasis enhancer networks"

4) Figure 1: The terminology the authors use, comparing the gene of interest to "Genome" is very confusing. They are not comparing to the entire genome but to all genes in the genome, which is not the same.

We have changed the word "genome" to "all genes in the genome" on the reviewer's suggestion.

5) Figure 2: Changes between X vs. Y is redundant (either changes between X and Y or changes in X vs. Y).

We assume that the reviewer is referring to Fig. 2B, which does not measure changes between X and Y, but changes in distribution between Acps and the control group. We have explained this in the figure legend.

The manuscript addresses a general question in evolutionary biology - do control regions diverge more quickly protein coding regions. The answer is that yes, they do, but this is actually not very surprising. The work is probably thus of more interest to people interested in the copulatory proteins or in the evolution of mating systems, than to people interested in broader evolutionary questions.

We appreciate this reviewer's recognition of the significance of our work and would like to point out that there are very few studies looking at promoter evolution as detailed in the introduction. Of particular relevance, our study using Acp genes allows us to directly test the impact of promoter mutations on the expression by comparing two alleles in male accessory glands of Drosophila hybrids. Male accessory glands consist of only two secretory cell types allowing us to study evolution of gene expression in a single cell type (Acps are either expressed in main cells or secondary cells). Amid this unique experimental set up we can conclude that promoter mutations can act dominant, in contrast to mutations in protein coding regions, which are generally recessive. Thus, our study is unique in pointing out a largely overseen aspect of gene evolution.

Reviewer 2

This manuscript explores promoter evolution of genes encoding seminal fluid proteins expressed in the male accessory gland of Drosophila and finds cis-regulatory changes underlie expression differences between species. Although these genes evolve rapidly it appears that the coding regions rarely show signs of positive selection inferring that changes in their expression and hence promoter sequences can underlie the evolution of their roles within and among species.

We thank Reviewer 2 for their thorough review, positive feedback on the importance of our work, and suggestions for improving the manuscript. We have addressed all points raised by the reviewer, including analysis of Acp coding region evolution, additional analyses of hybrid expression data, and improved the clarity of the text.

Figure 1 illustrates evidence that the promoter regions of these gene have accumulated more changes than other sampled genes from the Drosophila genome. While this convinces that the region upstream of the transcription start site has diverged considerably in sequence (grey line compared to black line), Figure 1A also suggests the "Genespan" region which includes the 5'UTR but presumably also part of the coding region is also highly diverged. It would be useful to see how the pattern extends into the coding region further to compare further to the promoter region (although Fig 1H does illustrate this more convincingly).

The reviewer raises an interesting point, and certainly all parts of genes evolve. Fig. 1A shows the evolutionary rates of Acps compared to the genome average from phyloP27way scores calculated from 27 insect species. Since these species are quite distant it is unsurprising that they show divergence in coding regions as well as promoter regions. In fact, we addressed whether promoter regions evolve fast in closely related Drosophila species in Fig. 1H compared to coding regions. We have included an additional analysis of coding region evolution in Figure 1B.

Figure 2 presents evidence for significant changes in (presumably levels of) expression of male accessory gland protein (AcP) genes and ribosomal proteins genes between pairs of species, which is reflected in the skew of expression compared to randomly selected genes.

Correct, we have rephrased the statement for clarity.

Figure 3 shows detailed analysis for 3 selected AcP genes with significantly diverged expression. The authors claim this shows 'substitution' hotspots in the promoter regions of all 3 genes but this could be better illustrated by extending the plots in B-D further upstream and downstream to compare to these regions.

We picked the 300-nucleotide promoter region for this analysis as it accumulated significant changes as shown in Fig. 1E-H, and extending the G plots (Fig. 3B-D) to regions with lower numbers of sequence changes would not substantially change the conclusion. Specifically, this analysis identifies sequence change hotspots within fast-evolving promoter regions, rather than comparing promoter regions to other genomic regions, as we previously addressed. The plot is based on a cumulative distribution function and the significant positive slope in the upstream region where promoters are located identifies a hotspot for accumulation of substitutions. There could be other hotspots, but the point being made is that significant hotspots consistently appear in the promoter region of these three genes.

Figure 4 shows the results of expression analysis in parental lines of each pair of species and F1 hybrids. However the results are very difficult to follow in the figure and in the relevant text. While the schemes in A, C. E and G are helpful, the gel images are not the best quality and interpretations confusing. An additional scheme is needed to illustrate hypothetical outcomes of trans change, cis change and transvection to help interpret the gels. On line 169 (presumably referring to panels D and F although C and D are cited on the next line) the authors claim that Obp56f and CG11598 'were more expressed in D. melanogaster compared to D. simulans' but in the gel image the D. sim band is stronger for both genes (like D. sechellia) compared to the D. mel band. The authors also claim that the patterns of expression seen in the F1s are dominant for one allele and that this must be because of transvection. I agree this experiment is evidence for cis-regulatory change. However the interpretation that it is caused by transvection needs more explanation/justification and how do the authors rule out that it is not a cis X trans interaction between the species promoter differences and differences in the transcription factors of each species in the F1? Also my understanding is that transvection is relatively rare and yet the authors claim this is the explanation for 2/4 genes tested.

We appreciate the reviewer's comments on Figure 4 and the opportunity to improve its clarity. To address these concerns, we have carefully checked the figure citations and corrected any inconsistencies.

The reviewer raises an important point about our interpretation of transvection. We have expanded our discussion of this result to consider why transvection is a plausible explanation for the observed dominance patterns and also consider cis x trans interactions between species-specific promoters and transcription factor binding. While rare, transvection likely has more relevance in hybrid regulatory contexts involving homologous chromosome pairing which we discuss this in the revised text.

Line 112 states that the melanogaster subgroup contains 5 species - this is incorrect - while this study looked at 5 species there are more species in this subgroup such as mauritiana and santomea.

We have corrected the statement about the number of species in the melanogaster subgroup.

Lines 131-134 could explain better what the conservation scores and their groupings mean and the rationale for this approach.

We have clarified what the conservation scores and their groupings mean and the rationale for this approach.

Line 162 - the meaning of the sentence starting on this line is unclear - it sounds very circular.

We have rephrased the statement for more clarity.

Line 168 should cite Fig 4 H instead of F.

We have amended citation of Fig 4F to H.

Reviewer 3

In this study, McQuarrie et al. investigate the evolution of promoters of genes encoding accessory gland proteins (Acps) in species within the D. melanogaster subgroup. Using computational analyses and available genomic and transcriptomic datasets, they demonstrate that promoter regions of Acp genes are highly diverse compared to the promoters of other genes in the genome. They further show that this diversification correlates with changes in gene expression levels between closely related species. Complementing these computational analyses, the authors conduct experiments to test whether differences in expression levels of four Acp genes with highly diverged promoter regions are maintained in hybrids of closely related species. They find that while two Acp genes maintain their expression level differences in hybrids, the other two exhibit dominance of one allele. The authors attribute these findings to transvection. Based on their data, they conclude that rapid evolution of Acp gene promoters, rather than changes in trans, drives changes in Acp gene expression that contribute to speciation.

We thank Reviewer 3 for their thorough review and suggestions. We further thank the reviewer for acknowledging the importance of our findings and for pointing out that it contributes to our understanding of speciation. We have thoroughly addressed all comments from the reviewer and significantly revised the manuscript. We believe that this has greatly improved the manuscript.

Unfortunately, the presented data are not sufficient to fully support the conclusions. While many of the concerns can be addressed by revising the text to moderate the claims and acknowledge the methodological limitations, some key experiments require repetition with more controls, biological replicates, and statistical analyses to validate the findings.

Specifically, some of the main conclusions heavily rely on the RT-PCR experiments presented in Figure 4, which analyze the expression of four Acp genes in hybrid flies. The authors use PCR and RFLP to distinguish species-specific alleles but draw quantitative conclusions from what is essentially a qualitative experiment. There are several issues with this approach. First, the experiment includes only two biological replicates per sample, which is inadequate for robust statistical analysis. Second, the authors did not measure the intensity of the gel fragments, making it impossible to quantify allele-specific expression accurately. Third, no control genes were used as standards to ensure the comparability of samples.

The gold standard for quantifying allele-specific expression is using real-time PCR methods such as TaqMan assays, which allow precise SNP genotyping. To address this major limitation, the authors should ideally repeat the experiments using allele-specific real-time PCR assays. This would provide a reliable and quantitative measurement of allele-specific expression.

If the authors cannot implement real-time PCR, an alternative (though less rigorous) approach would be to continue using their current method with the following adjustments:

• Include a housekeeping gene in the analysis as an internal control (this would require identifying a region distinguishable by RFLP in the control).

• Quantify the intensity of the PCR products on the gel relative to the internal standard, ensuring proper normalization.

• Increase the sample size to allow for robust statistical analysis.

These experiments could be conducted relatively quickly and would significantly enhance the validity of the study's conclusions.

We thank the reviewer for their detailed suggestions for improving the conclusions in Fig. 4. Indeed, incorporating a housekeeping gene as a control supports our results for qualitative analysis of gene expression in hybrids assessing each allele individually (Fig 4), and improves interpretation for non-experts. We have also quantified differential gene expression in hybrids between species alleles and the log2 fold change from D. melanogaster. In addition, we have included an additional analysis in the new Fig. 5 which analyses RNA-seq expression changes in D. melanogaster x D. simulans hybrid male accessory glands. We believe these additions have significantly improved the manuscript and its conclusions.

While the following comments are not necessarily minor, they can be addressed through revisions to the text without requiring additional experimental work. Some comments are more conceptual in nature, while others concern the interpretation and presentation of the experimental results. They are provided in no particular order.

1. A key limitation of this study is the use of RNA-seq datasets from whole adult flies for interspecies gene expression comparisons. Whole-body RNA-seq inherently averages gene expression across all tissues, potentially masking tissue-specific expression differences. While Acp genes are likely restricted to accessory glands, the non-Acp genes and the random gene sets used in the analysis may have broader expression profiles. As a result, their expression might be conserved in certain tissues while diverging in others- an aspect that whole-body RNA-seq cannot capture. The authors should acknowledge that tissue-specific RNA-seq analyses could provide a more precise understanding of expression divergence and potentially reveal reduced conservation when considering specific tissues independently.

We have added a section discussing the limitations in gene expression analysis in the discussion. In addition, we have included an additional Figure analysing gene expression in hybrid male accessory glands (Fig. 5).

1. The statement in line 128, "Consistent with this model," does not accurately reflect the findings presented in Figures 2A and B. Specifically, the data in Figure 2A show that Acp gene expression divergence is significantly different from the divergence of non-Acp genes or a random sample only in the comparison between D. melanogaster and D. simulans. However, when these species are compared to D. yakuba, Acp gene expression divergence aligns with the divergence patterns of non-Acp genes or random samples. In contrast, Figure 2B shows that the distribution of expression changes is skewed for Acp genes compared to random control samples when D. melanogaster or D. simulans are compared to D. yakuba. However, this skew is absent when the two D. melanogaster and D. simulans are compared. Therefore, the statement in line 128 should be revised to accurately reflect these nuanced results and the trends shown in Figure 2A and B.

We have updated the statement for clarity. Here, the percentage of Acps showing significant gene expression changes is greater between more closely related species, but the distribution of expression changes increases between more distantly related species.

1. The statement in lines 136-138, "Acps were enriched for significant expression changes in the faster evolving group across all species," while accurate, overlooks a key observation. This trend was also observed in other groups, including those with slower evolving promoters, in some of the species' comparisons. Therefore, the enrichment is not unique to Acps with rapidly evolving promoters, and this should be explicitly acknowledged in the text.

This is a valid point, and we have updated this statement as suggested.

1. It would be helpful for the authors to explain the meaning of the d score at the beginning of the paragraph starting in line 131, to ensure clarity for readers unfamiliar with this metric.

This scoring method is described in the methods sections, and we have now included reference to thorough explanation of how d was calculated at the indicated section.

1. In Figure 2C-E - the title of the Y-axis does not match the text. If it represents the percentage of genes with significant expression changes, as in Figure 2A, the discrepancies between the percentages in this figure and those in Figure 2A need to be addressed.

We have updated the method used to categorise significant changes in gene expression in the text and the figure legend for clarity.

1. The experiment in Figure 3 needs a better explanation in the text. What is the analysis presented in Figure 3B-D. How many species were compared?

We have added additional details in the results section and an explanation of how sequence change hotspots were calculated in the results section is available.

1. The concept of transvection should be omitted from this manuscript. First, the definition provided by the authors is inaccurate. Second, even if additional experiments were to convincingly show that one allele in hybrid animals is dominant over the other, there are alternative explanations for this phenomenon that do not involve transvection. The authors may propose transvection as a potential model in the discussion, but they should do so cautiously and explicitly acknowledge the possibility of other mechanisms.

We have updated the text to more conservatively discuss transvection, moving this to the discussion section with additional possibilities discussed.

1. The statement at the end of the introduction is overly strong and would benefit from more cautious phrasing. For instance, it could be reworded as: "These findings suggest that promoter changes, rather than genomic background, play a significant role in driving expression changes, indicating that promoter evolution may contribute to the rise of new species."

We have reworded this line following the reviewer's suggestion.

1. Line 32 of the abstract: The term "Acp" is introduced without explaining what it stands for. Please define it as "Accessory gland proteins (Acp)" when it first appears.

We have updated the manuscript to define Acp where it is first mentioned.

1. Line 61: The phrase "...through relaxed,..." is unclear. Specify what is relaxed (e.g., "relaxed selective pressures").

We have included description of relaxed selective pressures.

1. The sentence in lines 74-76, starting in "Interestingly,...." Needs revision for clarity.

We have removed the word interestingly.

1. Line 112: Revise "we focused on the melanogaster subgroup which is made up of five species" to: "we focused on the melanogaster subgroup, which includes five species."

We have made this change in the text.

1. In line 144 use the phrase "promoter conservation" instead of "promoter evolution"

We have updated the phrasing.

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:

In this study, McQuarrie et al. investigate the evolution of promoters of genes encoding accessory gland proteins (Acps) in species within the D. melanogaster subgroup. Using computational analyses and available genomic and transcriptomic datasets, they demonstrate that promoter regions of Acp genes are highly diverse compared to the promoters of other genes in the genome. They further show that this diversification correlates with changes in gene expression levels between closely related species. Complementing these computational analyses, the authors conduct experiments to test whether differences in expression levels of four Acp genes with highly diverged promoter regions are maintained in hybrids of closely related species. They find that while two Acp genes maintain their expression level differences in hybrids, the other two exhibit dominance of one allele. The authors attribute these findings to transvection. Based on their data, they conclude that rapid evolution of Acp gene promoters, rather than changes in trans, drives changes in Acp gene expression that contribute to speciation.

Major comments:

Unfortunately, the presented data are not sufficient to fully support the conclusions. While many of the concerns can be addressed by revising the text to moderate the claims and acknowledge the methodological limitations, some key experiments require repetition with more controls, biological replicates, and statistical analyses to validate the findings.

Specifically, some of the main conclusions heavily rely on the RT-PCR experiments presented in Figure 4, which analyze the expression of four Acp genes in hybrid flies. The authors use PCR and RFLP to distinguish species-specific alleles but draw quantitative conclusions from what is essentially a qualitative experiment. There are several issues with this approach. First, the experiment includes only two biological replicates per sample, which is inadequate for robust statistical analysis. Second, the authors did not measure the intensity of the gel fragments, making it impossible to quantify allele-specific expression accurately. Third, no control genes were used as standards to ensure the comparability of samples.

The gold standard for quantifying allele-specific expression is using real-time PCR methods such as TaqMan assays, which allow precise SNP genotyping. To address this major limitation, the authors should ideally repeat the experiments using allele-specific real-time PCR assays. This would provide a reliable and quantitative measurement of allele-specific expression.

If the authors cannot implement real-time PCR, an alternative (though less rigorous) approach would be to continue using their current method with the following adjustments:

• Include a housekeeping gene in the analysis as an internal control (this would require identifying a region distinguishable by RFLP in the control).
• Quantify the intensity of the PCR products on the gel relative to the internal standard, ensuring proper normalization.
• Increase the sample size to allow for robust statistical analysis. These experiments could be conducted relatively quickly and would significantly enhance the validity of the study's conclusions.

Minor comments

While the following comments are not necessarily minor, they can be addressed through revisions to the text without requiring additional experimental work. Some comments are more conceptual in nature, while others concern the interpretation and presentation of the experimental results. They are provided in no particular order. 1. A key limitation of this study is the use of RNA-seq datasets from whole adult flies for interspecies gene expression comparisons. Whole-body RNA-seq inherently averages gene expression across all tissues, potentially masking tissue-specific expression differences. While Acp genes are likely restricted to accessory glands, the non-Acp genes and the random gene sets used in the analysis may have broader expression profiles. As a result, their expression might be conserved in certain tissues while diverging in others- an aspect that whole-body RNA-seq cannot capture. The authors should acknowledge that tissue-specific RNA-seq analyses could provide a more precise understanding of expression divergence and potentially reveal reduced conservation when considering specific tissues independently. 2. The statement in line 128, "Consistent with this model," does not accurately reflect the findings presented in Figures 2A and B. Specifically, the data in Figure 2A show that Acp gene expression divergence is significantly different from the divergence of non-Acp genes or a random sample only in the comparison between D. melanogaster and D. simulans. However, when these species are compared to D. yakuba, Acp gene expression divergence aligns with the divergence patterns of non-Acp genes or random samples. In contrast, Figure 2B shows that the distribution of expression changes is skewed for Acp genes compared to random control samples when D. melanogaster or D. simulans are compared to D. yakuba. However, this skew is absent when the two D. melanogaster and D. simulans are compared. Therefore, the statement in line 128 should be revised to accurately reflect these nuanced results and the trends shown in Figure 2A and B. 3. The statement in lines 136-138, "Acps were enriched for significant expression changes in the faster evolving group across all species," while accurate, overlooks a key observation. This trend was also observed in other groups, including those with slower evolving promoters, in some of the species' comparisons. Therefore, the enrichment is not unique to Acps with rapidly evolving promoters, and this should be explicitly acknowledged in the text. 4. It would be helpful for the authors to explain the meaning of the d score at the beginning of the paragraph starting in line 131, to ensure clarity for readers unfamiliar with this metric. 5. In Figure 2C-E - the title of the Y-axis does not match the text. If it represents the percentage of genes with significant expression changes, as in Figure 2A, the discrepancies between the percentages in this figure and those in Figure 2A need to be addressed. 6. The experiment in Figure 3 needs a better explanation in the text. What is the analysis presented in Figure 3B-D. How many species were compared? 7. The concept of transvection should be omitted from this manuscript. First, the definition provided by the authors is inaccurate. Second, even if additional experiments were to convincingly show that one allele in hybrid animals is dominant over the other, there are alternative explanations for this phenomenon that do not involve transvection. The authors may propose transvection as a potential model in the discussion, but they should do so cautiously and explicitly acknowledge the possibility of other mechanisms. 8. The statement at the end of the introduction is overly strong and would benefit from more cautious phrasing. For instance, it could be reworded as: "These findings suggest that promoter changes, rather than genomic background, play a significant role in driving expression changes, indicating that promoter evolution may contribute to the rise of new species."

Text edits:

Throughout the manuscripts there are incomplete sentences and sentences that are not clear. Below is a list of corrections:

1. Line 32 of the abstract: The term "Acp" is introduced without explaining what it stands for. Please define it as "Accessory gland proteins (Acp)" when it first appears.
2. Line 61: The phrase "...through relaxed,..." is unclear. Specify what is relaxed (e.g., "relaxed selective pressures").
3. The sentence in lines 74-76, starting in "Interestingly,...." Needs revision for clarity.
4. Line 112: Revise "we focused on the melanogaster subgroup which is made up of five species" to: "we focused on the melanogaster subgroup, which includes five species."
5. In line 144 use the phrase "promoter conservation" instead of "promoter evolution"

Significance

This study addresses an important question in evolutionary biology: how seminal fluid proteins achieve rapid evolution despite showing limited adaptive changes in their coding regions. By focusing on accessory gland proteins (Acps) and examining their promoter regions, the authors suggest promoter-driven evolution as a potential mechanism for rapid seminal fluid protein diversification. While this hypothesis is intriguing and can contribute to our understanding of speciation, more rigorous analysis and experimental validation would be needed to support the conclusions. The revised manuscript can be of interest to fly geneticists and to scientists in the fields of gene regulation and evolution.

Keywords for my expertise: Enhancers, transcriptional regulation, development, evolution, Drosophila.

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

The authors test the hypothesis that promoters of genes involved in insect accessory glands evolved more rapidly than other genes in the genome. They test this using a number of computational and experimental approaches, looking at different species within the Drosophila melanogaster complex. The authors find an increased amount of sequence divergence in promoters of accessory gland proteins. They show that the expression levels of these proteins are more variable among species than randomly selected proteins. Finally, they show that within interspecific hybrids, each copy of the gene maintains its species-specific expression level.

The work is done with expected standards of controls and analyses. The claims are supported by the analysis. My main criticism of the manuscript has to do not with the experiments or conclusion themselves but with the presentation. The manuscript is just not very well written, and following the logic of the arguments and results is challenging. The problem begins with the Abstract, which is representative of the general problems with the manuscript. The Abstract begins with general statements about the evolution of seminal fluid proteins, but then jumps to accessory glands and hybrids, without clarifying what taxon is being studied, and what hybrids they are talking about. Then, the acronym Acp is introduced without explanation. The last two sentences of the Abstract are very cumbersome and one has to reread them to understand how they link to the beginning of the Abstract.

More generally, if this reviewer is to be seen as an "average reader" of the paper, I really struggled through reading it, and did not understand many of the arguments or rationale until the second read-through, after I had already read the bottom line. The paragraph spanning lines 71-83 is another case in point. It is composed of a series of very strongly worded sentences, almost all starting with a modifier (unexpectedly, interestingly, moreover), and supported by citations, but the logical flow doesn't work. Again, reading the paragraph after I knew where the paper was going was clearer, but on a first read, it was just a list of disjointed statements.

Since most of the citations are from the authors' own work, I suspect they are assuming too much prior understanding on the part of the reader. I am sure that if the authors read through the manuscript again, trying to look through the eyes of an external reader, they will easily be able to improve the flow and readability of the text.

More specific comments:

1. In the analysis of expression level differences, it is not clear what specific stage / tissue the levels taken from the literature refer to. Could it be that the source of the data is from a stage or tissue where seminar fluid proteins will be expressed with higher variability in general (not just inter-specifically) and this could be skewing the results? Please add more information on the original source of the data and provide support for their validity for this type of comparison.
2. The sentence spanning lines 155-157 needs more context.
3. Line 203-204: What are multi-choice enhancers?
4. Figure 1: The terminology the authors use, comparing the gene of interest to "Genome" is very confusing. They are not comparing to the entire genome but to all genes in the genome, which is not the same.
5. Figure 2: Changes between X vs. Y is redundant (either changes between X and Y or changes in X vs. Y).

Significance

The manuscript addresses a general question in evolutionary biology - do control regions diverge more quickly protein coding regions. The answer is that yes, they do, but this is actually not very surprising. The work is probably thus of more interest to people interested in the copulatory proteins or in the evolution of mating systems, than to people interested in broader evolutionary questions.

DOI: 10.1038/s41467-025-62205-y

Resource: Mutant Mouse Regional Resource Center (RRID:SCR_002953)

Curator: @AleksanderDrozdz

SciCrunch record: RRID:SCR_002953

What is this?

DOI: 10.1038/s41467-025-62205-y

Resource: None

Curator: @AleksanderDrozdz

SciCrunch record: RRID:MMRRC_068495-MU

What is this?

DOI: 10.1038/s41467-025-62205-y

Resource: (Bethyl Cat# A300-587A, RRID:AB_495510)

Curator: @scibot

SciCrunch record: RRID:AB_495510

What is this?

DOI: 10.1038/s41467-025-62205-y

Resource: (Thermo Fisher Scientific Cat# A301-038A, RRID:AB_817992)

Curator: @scibot

SciCrunch record: RRID:AB_817992

What is this?

DOI: 10.1016/j.celrep.2025.116045

Resource: (IMSR Cat# JAX_006852,RRID:IMSR_JAX:006852)

Curator: @dhovakimyan1

SciCrunch record: RRID:IMSR_JAX:006852

What is this?

DOI: 10.1038/s44319-025-00472-y

Resource: (RRID:CVCL_6996)

Curator: @scibot

SciCrunch record: RRID:CVCL_6996

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678586

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678585

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_2861257

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678584

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678582

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678540

What is this?

DOI: 10.1038/s41596-025-01215-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_3678583

What is this?

DOI: 10.1016/j.celrep.2025.116061

Resource: RRID:Addgene_48138

Curator: @scibot

SciCrunch record: RRID:Addgene_48138

What is this?

A nice way to describe linear functions.

Author response:

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

Reviewer #1 (Public review):

This study investigates how ant group demographics influence nest structures and group behaviors of Camponotus fellah ants, a ground-dwelling carpenter ant species (found locally in Israel) that build subterranean nest structures. Using a quasi-2D cell filled with artificial sand, the authors perform two complementary sets of experiments to try to link group behavior and nest structure: first, the authors place a mated queen and several pupae into their cell and observe the structures that emerge both before and after the pupae eclose (i.e., "colony maturation" experiments); second, the authors create small groups (of 5,10, or 15 ants, each including a queen) within a narrow age range (i.e., "fixed demographic" experiments) to explore the dependence of age on construction. Some of the fixed demographic instantiations included a manually induced catastrophic collapse event; the authors then compared emergency repair behavior to natural nest creation. Finally, the authors introduce a modified logistic growth model to describe the time-dependent nest area. The modification introduces parameters that allow for age-dependent behavior, and the authors use their fixed demographic experiments to set these parameters, and then apply the model to interpret the behavior of the colony maturation experiments. The main results of this paper are that for natural nest construction, nest areas, and morphologies depend on the age demographics of ants in the experiments: younger ants create larger nests and angled tunnels, while older ants tend to dig less and build predominantly vertical tunnels; in contrast, emergency response seems to elicit digging in ants of all ages to repair the nest.

We sincerely thank Reviewer #1 for the time and effort dedicated to our manuscript's detailed review and assessment. The revision suggestions were constructive, and we have provided a point-by-point response to address them.

Reviewer #2 (Public review):

I enjoyed this paper and the approach to examining an accepted wisdom of ants determining overall density by employing age polyethism that would reduce the computational complexity required to match nest size with population (although I have some questions about the requirement that growth is infinite in such a solution). Moreover, the realization that models of collective behaviour may be inappropriate in many systems in which agents (or individuals) differ in the behavioural rules they employ, according to age, location, or information state. This is especially important in a system like social insects, typically held as a classic example of individual-as-subservient to whole, and therefore most likely to employ universal rules of behaviour. The current paper demonstrates a potentially continuous age-related change in target behaviour (excavation), and suggests an elegant and minimal solution to the requirement for building according to need in ants, avoiding the invocation of potentially complex cognitive mechanisms, or information states that all individuals must have access to in order to have an adaptive excavation output.

We sincerely thank reviewer #2 for the time and effort dedicated to our manuscript's detailed review and assessment. We have provided a point-by-point response to the reviewer's comments, which we have incorporated into the revised version of the manuscript.

The only real reservation I have is in the question of how this relationship could hold in properly mature colonies in which there is (presumably) a balance between the birth and death of older workers. Would the prediction be that the young ants still dig, or would there be a cessation of digging by young ants because the area is already sufficient? Another way of asking this is to ask whether the innate amount of digging that young ants do is in any way affected by the overall spatial size of the colony. If it is, then we are back to a problem of perfect information - how do the young ants know how big the overall colony is? Perhaps using density as a proxy? Alternatively, if the young ants do not modify their digging, wouldn't the colony become continuously larger? As a non-expert in social insects, I may be misunderstanding and it may be already addressed in the citations used.

We thank the reviewer for this interesting question. We find that the nest excavation is predominantly performed by the younger ants in the nest, and the nest area increase is followed by an increase in the population. However, if the young ants dig unrestricted, this could result in unnecessary nest growth as suggested by reviewer #2. Therefore, we believe that the innate digging behavior of ants could potentially be regulated by various cues such as;

(a) Density-based: If the colony becomes less dense as its area expands, this could serve as a feedback signal for young ants to reduce or stop digging, as described in references (25, 29, 30).

(b) Pheromone depositions: If the colony reaches a certain population density, pheromone signals could inhibit further digging by young ants, references (25, 29), or space usage as a proxy for the nest area. 

Thus, rather than perfect information, decentralized control, and digging-based local cues probably regulate the level of age-dependent digging, without the ants needing to estimate the overall colony size or nest area.

In any case, this is an excellent paper. The modelling approach is excellent and compelling, also allowing extrapolation to other group sizes and even other species. This to me is the main strength of the paper, as the answer to the question of whether it is younger or older ants that primarily excavate nests could have been answered by an individual tracking approach (albeit there are practical limitations to this, especially in the observation nest setup, as the authors point out). The analysis of the tunnel structure is also an important piece of the puzzle, and I really like the overall study.

We thank the reviewer for the comments. We completely agree that individual tracking of ants within our experimental setup would have been the ideal approach, but we were limited by technical and practical limitations of the setup, as pointed out by the reviewer, such as; 

(a) Continuous tracking of ants in our nests would have required a camera to be positioned at all times in front of the nest, which necessitates a light background. Since Camponotus fellah ants are subterranean, we aimed to allow them to perform nest excavation in conditions as close to their natural dark environment as possible. Additionally, implementing such a system in front of each nest would have reduced the sample sizes for our treatments.

(b) The experimental duration of our colony maturation and fixed demographics experiments extended for up to six months (unprecedented durations in these kinds of measurements). These naturally limited our ability to conduct individual tracking while maintaining the identity of each ant based on the current design.

These details are described in detail within the revised version of the manuscript.

Reviewer #3 (Public review):

Summary:

In this study, Harikrishnan Rajendran, Roi Weinberger, Ehud Fonio, and Ofer Feinerman measured the digging behaviours of queens and workers for the first 6 months of colony development, as well as groups of young or old ants. They also provide a quantitative model describing the digging behaviours and allowing predictions. They found that young ants dig more slanted tunnels, while older ants dig more vertically (straight down). This finding is important, as it describes a new form of age polyethism (a division of labour based on age). Age polyethism is described as a "yes or no" mechanism, where individuals perform or not a task according to their age (usually young individuals perform in-nest tasks, and older ones foraging). Here, the way of performing the task is modified, not only the propensity to carry it or not. This data therefore adds in an interesting way to the field of collective behaviours and division of labour.

The conclusions of the paper are well supported by the data. Measurements of the same individuals over time would have strengthened the claims.

We sincerely thank reviewer #3 for the time and effort dedicated to our manuscript's detailed review and assessment. We completely agree with the reviewer’s comments on the measurements of the same individuals over time, however, we were limited by the technical and experimental limitations as described above and pointed out by reviewer #2.

Strengths:

I find that the measure of behaviour through development is of great value, as those studies are usually done at a specific time point with mature colonies. The description of a behaviour that is modified with age is a notable finding in the world of social insects. The sample sizes are adequate and all the information clearly provided either in the methods or supplementary.

We thank reviewer #3  for this assessment.

Weaknesses:

I think the paper is failing to take into consideration or at least discuss the role of inter-individual variabilities. Tasks have been known to be undertaken by only a few hyper-active individuals for example. Comments on the choice to use averages and the potential roles of variations between individuals are in my opinion lacking. Throughout the paper wording should be modified to refer to the group and not the individuals, as it was the collective digging that was measured. Another issue I had was the use of "mature colony" for colonies with very few individuals and only 6 months of age. Comments on the low number of workers used compared to natural mature colonies would be welcome.

Regarding the main comment 1

We completely agree with the reviewer’s comment on considering inter-individual variability based on activity levels. We have discussed how individual morphological variability could influence digging behavior (references: 28, 31), and we will elaborate further on this aspect in future revisions.

Regarding the main comment 2:

The term ‘colony maturation’ in our study refers to the progressive development of colonies from a single queen, distinguishing it from experiments that begin with pre-established, demographically stable colonies. We provide a detailed explanation for this terminology in the revised version of the manuscript. We were practically limited by the continuation of the experiments for more than 6 months of age, predominantly due to the stability of nests, as they were made with a sand-soil mix. We also acknowledge that the colony sizes attained in our maturation experiments may be smaller than those of naturally matured colonies. This trend was observed generally in lab-reared colonies and could be attributed to differences in microclimatic conditions, foraging opportunities, space availability, and other factors. We have explicitly described these details in the revised version of the manuscript.

Reviewer #1 (Recommendations for the authors):

The experimental design is fantastic. The large quasi-2D should allow for the direct visualization of the movements of individuals and the creation of the nest, and the inclusion of non-workers (specifically, a mated queen and pupae) is new and important. However, I have some questions and concerns about the results, as outlined below. Also, I found the paper difficult to read, and the connections between the various experiments and the model were not always clear. 

We thank the reviewer for the time and effort dedicated to reviewing our manuscript. We have modified the manuscript substantially to address the comments and readability. 

The assumption that the digging rate is constant across ants may be a strong one. Previous work (see, for instance, Aguilar, et al, Science 2018) has demonstrated a very heterogeneous workload distribution among ants. I am not sure what implications that may have for the results here, but the authors should comment on this choice. Related to the point above, given a constant digging rate, the variation in digging is attributed to an age-dependent "desired target area". Can the authors comment on the implications of this, specifically in contrast to a variable digging rate? The distinction between digging rate differences and target area differences seems to be important for the authors. However, the way this is presented, it is difficult to fully understand or appreciate this importance and its implications. What is the consequence of this difference, and why is this important?

We apologize to the reviewer for the confusion.

Our model does not assume that the digging rate (da/dt, Equation 1) remains constant throughout the experiment. Instead, we only treat the basal digging rate (r) as a constant.

The variable digging rate (da/dt, Equation 1) is derived by multiplying the basal rate constant (r) by the term (1 - a/a_age), which accounts for deviations from the age-dependent target area that the ants aim to achieve. This makes the actual digging rate dynamic, as it responds to changes in excavated area (e.g., expansion or rapid collapse)

For example, according to our model (Equation 1), two ants with the same basal digging rate (r) may exhibit markedly different actual digging rates at a given time if they differ in age. This occurs because the variable digging rate (da/dt) depends not only on ‘r’ but also on the age-dependent term (1 - a/a_age). Also, we emphasize that the use of a basal digging rate constant aligns with prior studies (refs. 24, 29, 30).

In our work, we demonstrate that after a collapse event, ants of all ages dig at rates comparable to those observed in the initial (pre-collapse) phase of the experiment. This occurs because the ants are far from their age-dependent target area, effectively resetting their digging behavior. By comparing maximum digging rates pre- and post-collapse, we provide strong empirical evidence that this rate is age-independent (SI Fig. 6A, 6B), supporting the conclusion that the basal digging rate constant (r) is a fundamental property of the ants' behavior, unaffected by age.

We agree with the reviewer that individual tracking of ants within our experimental setup would have been the ideal approach. Then, we could have taken the inter-individual variability of the digging activity into account. However, we were limited to doing so by the technical and practical limitations of the setup, such as; 

(a) Continuous tracking of ants in our nests would have required a camera to be positioned at all times in front of the nest, which necessitates a light background. Since Camponotus fellah ants are subterranean, we aimed to allow them to perform nest excavation in conditions as close to their natural dark environment as possible. Additionally, implementing such a system in front of each nest would have reduced the sample sizes for our treatments.

(b) The experimental duration of our colony maturation experiments extended for up to six months (unprecedented durations in these kinds of measurements). These naturally limited our ability to conduct individual tracking while maintaining the identity of each ant based on the current design.

In light of these points, the following lines are added to the discussion (line numbers: 283-295), signifying the above points:

“Our age-dependent model demonstrates that the digging behavior in Camponotus fellah is governed by a basal digging rate constant (r) modulated by the age-dependent feedback (1 − a/aage). Crucially, we show that after a collapse, the maximum digging rates return to their pre-collapse levels, suggesting that this basal rate ’r’ represents an age-independent ceiling on how fast ants can dig, regardless of age or context (SI Fig. 6 A, B). Previous studies have demonstrated both homogeneous and heterogeneous workload distribution, with varying digging rates among ants (24, 29, 30, 35). Studies showing heterogeneous workload distribution relied on continuous individual tracking of ants to quantify digging rates (35). However, this approach was not feasible in our current design due to the experimental durations of both our colony maturation and fixed demographics experiments. Additionally, sample size requirements naturally limited our ability to conduct continuous individual tracking during nest construction in our study. Thus, based on empirical measurements from our fixed-demographics experiments and supported by the age-independent post-collapse digging rates, we adopted a constant basal digging rate for simulating our age-dependent model—an assumption aligned with both prior literature and the collective dynamics observed in our system (24,29,30)”.

Model: as presented, the model seems to lack independent validation. The model seems to have built-in that there is an age-dependent target area, and this is what is recovered from the model. I am failing to see what is learned from the model that the experiments do not already show. Also, the model has no ant interactions, though ants are eusocial and group size is known to have a large effect on behavior (this is acknowledged by the authors at the beginning of the discussion). Can the authors comment on this?My recommendation would be to remove the model from this paper or improve the text to address the above comments.

We did not draw the conclusion of the age-dependent target area from our model. We used the fixed demographics experiments to quantify the age-dependent area target as a function of the age of individuals. We then used this age-dependent area target in our model to quantify the excavation dynamics of the colony maturation experiments, where ants span a variety of ages, as the nest population changes over time, resulting in natural variation in the ages of individuals within the nest.  These results could not have been obtained by performing any of the individual experiments, whether colony maturation or the fixed demographics, young or old, on their own. The need for different age demographics was crucial to quantify the age-dependent effects in nest excavation, which were lacking in previous studies. 

First, the age-dependent model provides a very good estimate for the natural growth of the nest.  More importantly, after fixing an age threshold of 56 days (mean + standard deviation of the young ant age), the model provides an estimate of which ants are doing the majority of the digging during natural nest expansion. This teaches us that during natural expansion, the older ants are far from their density target and therefore do not engage in any substantial digging, which is shown in Figure 4. C. 

On the other hand, the younger ants are close to their area targets and induced to dig. Indeed, the target area fitted for the age-independent model closely approximates the empirically measured age-dependent target when extrapolated to very young ants. This provides further support for the idea that, in the colony maturation experiments, the youngest ants are responsible for most of the digging.

Our model is a simple analytical model, inspired by earlier models that used a fixed area target (such as density models) for nest construction. However, because we knew the precise age of workers in our experiments, we were able to obtain age-dependent area targets, thereby challenging the use of a constant area target (as employed in prior studies) in light of our findings from the fixed demographics of young and old colonies.

Empirically Quantifiable Parameters: We wanted our model to have empirically quantifiable parameters. Since we did not continuously record the experiment, we could not quantify agent-agent interactions, pheromonal depositions, or similar factors.

Minimal Model Design: We aimed to keep the model as minimal as possible, which is why we did not include complex interactions such as those found in continuous tracking experiments.

However, the model does set up some interesting hypotheses that could easily be tested with the experimental setup (e.g., marking the ants / tracking individual activity levels). For instance, it is hypothesized that older ants dig less often, but when they do dig, they do so at the same rate. Given the 2D setup, the authors could track individual ants and test this hypothesis. Also, if the desired target area does decrease with age, the authors could verify this hypothesis by placing older ants into arenas with different-sized pre-formed nests to observe how structure is changed to achieve the desired area/ant.

We thank the reviewer for this comment.

We believe that the confusion with the usage of a constant basal digging rate is resolved now. To briefly reiterate, ants dig at variable rates that can be decomposed to a (constant on short time scales but age-dependent) basal rate times the (variable) distance from the density target. The suggested experiments are beyond the scope of our current study, and further studies could utilize the suggested experimental design with better time-resolved imaging for individual ant tracking that could verify the predictions from our model. 

Specific comments:

Title:

The title suggests a broad result, yet the study focuses on one ant species. Please modify the title to more accurately reflect the scope of the work.

We thank the reviewer for the comment.

The title is modified as “Colony demographics shape nest construction in Camponotus fellah ants.”

Introduction:

Important information and context are missing about this ant species. For instance, please add the following about this species in the introduction:

What is their natural habitat and substrate? How does the artificial soil compare?

What is their (rough) colony size? [later, discuss experiment group size choice and potential insights/limitations of results when applied to the natural system].

The details have been added to the introduction (line numbers : 49-55) and the materials and methods section (Study species).

“Camponotus fellah ants are native to the Near East and North Africa, particularly found in countries like Israel, Egypt, and surrounding arid and semi-arid regions, where they prefer to nest in moist, decaying wood, including tree trunks, branches, or stumps (49,50). The species lives in monogynous colonies with tens to thousands of individuals. Nests are commonly found in a sand-loamy mix, which is a combination of sand, soil, clay, or gravel, providing structural stability and moisture retention (51). They are typically found under rocks, in the crevices of dried vegetation, or dry, sandy soils, sometimes in areas with loose gravel, with a colony size ranging from tens to thousands of workers”.

What is the natural life expectancy of a worker? A queen? [later, discuss fixed demographic age choices in this context and/or why were age ranges chosen for experiments?].

The lifespan of ants, including both queens and workers, varies significantly based on caste, species, and environmental conditions.

(1) Queen Longevity: From the literature, Camponotus fellah queens can live up to 20 years, with one documented case reaching 26 years (50). 

(2) Worker Longevity: In contrast to queens, the lifespan of workers is much shorter. Lab studies on Camponotus fellah (82) and other Camponotus species (83) suggest that workers can live for several months depending on environmental conditions, colony health, and caste-specific roles (e.g., minor vs. major workers)

(3) Laboratory vs. Natural Conditions: Worker longevity is highly variable between laboratory and natural conditions

Therefore, in the context of the old worker lifespan in our experiments, ~200 days (roughly 6–7 months), we strongly believe that the worker lifespan used in our experiments represents a substantial portion of a worker's expected life. While exact figures for C. fellah workers are unavailable, inferences from related species suggest that workers nearing 200 days are approaching the latter stages of their lifespan, making them meaningfully "old". 

The details are added to the main text (line numbers: 124-127) and discussion (line numbers: 278-282).

Why was this species chosen? Convenience, or is there something special about this species that the readers should know? Specifically, is there something that might make the results more general or of broader interest?

Camponotus fellah was chosen for this study because it is native to Israel, making it convenient to collect and maintain in the lab. Additionally, its nuptial flights occur close to the study location, ensuring a steady supply of colonies. We were able to provide them with a nesting substrate similar to what they naturally use, as their nests are typically found in a sand-loamy mix, similar to the sand-soil mix in our artificial nests. This was possible because we had the opportunity to observe their habitat and nesting behavior in the wild, allowing us to gather preliminary information on their natural nesting conditions.

Results:

Line 60: "several brood items" - how many exactly? Was this consistent across experiments? Do mated queens ever produce more pupae during the experiments?

Yes, the number of brood items (5) was added consistently across the experiments. Additionally, the mated queen did produce pupae during the course of the experiments, which was evident from the noticeable increase in the number of workers in the nest. This was significantly higher than the number of brood items present at the start of the study.

The above points are added to the section (line numbers : 68-69).

Figure 1: Panel A - The food ports are never mentioned in the text. Are the ants fed during the experiments? If so, what? With what frequency? Is the water column replenished/maintained? If so, how and how often? panel C - how long did this experiment last?

We thank the reviewer for pointing this out. We have now updated the nest maintenance section in the Materials and Methods (line numbers : 349-354) part to include all the necessary details and clarifications.

“We provided food to the ants ad libitum through three separate tubes containing water, 20 % sucrose water, and protein food. The protein mixture included egg powder, tuna, prawns, honey, agar, and vitamins. Each of the three tubes was filled with 5 ml of their respective contents and sealed with a cotton stopper to prevent overflow. The tubes were positioned at a slight angle and connected using a custom-made plexiglass adapter to facilitate the flow of liquids. These tubes were replenished once depleted, and regularly replaced once the nest maintenance was carried out bi-weekly.”

Line 76: "...excavation was commenced by the founding queen". How were the queen and pupae introduced into the system?

We initiated colony maturation experiments by introducing a single mated queen and several brood items (pupae) at random positions on the soil layer of the nest (line numbers : 68-69)

Line 87: Please provide bounds for 11cm2/ant value. Is there any biological or physical justification for this number?

We thank the reviewer for the suggestion. We have now provided the bounds as requested (line numbers : 97-101). 

We were unable to pinpoint a specific biological justification based solely on this treatment. However, on extrapolating the age-dependent area fit we derived from the fixed demographics experiment, we found that at the age of 1 day, an ant has a target area of approximately 11.17 cm², which is the largest age-dependent area target possible within our experimental setup.

From the colony maturation experiment, we obtained the value of  11.6 (±1.15) cm² as the area per ant. The consistency between the area per ant obtained from two completely different treatments across different colonies yielded similar results. We propose that under standardized conditions, a 1-day-old ant has a theoretical maximum target area of 11.17 cm²—the highest value observed in our experimental framework.

Lines 98-99: "one straightforward possibility would be that newborn ants are the ones that dig". This statement contradicts the results presented in Figures 1 and S1 - the population increase seems to occur at least a few days before increased excavation in nearly all cases.

We apologize for any confusion caused by our initial phrasing. To clarify, we proposed that a lag likely exists between population growth and nest area expansion. This lag could arise from two sequential processes: (1) newborn ants require time to mature and become active (first delay), and (2) digging to expand the nest takes additional time (second delay; estimated at ~10 days from the cross-correlation analysis). Thus, our results suggest that it is not the population that lags behind the area, but rather the area that lags behind the population, as demonstrated in Figures 2D and SI. Figure. S1.

The sentence “one straightforward possibility would be that newborn ants are the ones that dig” is modified as below (line numbers : 112-119) to prevent further confusion.

“One possible explanation is that, although all ants are capable of digging, it is primarily the newly emerged ants who perform this task. In this case, nest expansion would lag behind colony growth due to two delays: first, the time needed for young ants to mature enough to begin digging, and second, the physical time required to excavate additional space (e.g., around 10 days). This mechanism could eliminate the need for ants to assess overall colony density, as each new group of active workers simply enlarges the nest as they become ready. An alternative possibility is that all ants, regardless of age, respond to increased density by initiating excavation. In that scenario, nest expansion would follow more immediately after the emergence of new individuals, making delays less prominent (24, 29, 30)”.

Line 105: How do group sizes compare to natural colony size? Line 106: How do "young" and "old" classifications compare to natural life expectancy?

We have already addressed this question in an earlier comment. The details are added to the main text (line numbers: 124-127) and discussion (line numbers: 278-282).

Line 118-119: How are nests artificially collapsed?

We have added a new section in the Materials and Methods section that describes the nest collapsing procedure (Nest artificial collapse - line numbers : 386-399).

Figure 2 Panel A: The white dotted line is nearly impossible to see. Please use a more visible color.

We thank the reviewer for the comment.

We changed the solid circles to violet and the dotted line color to continuous white.

Figure 3: The use of circle markers as post-collapse recovery in young and old as well as old pre-collapse is confusing. Use different symbols for old pre-collapse vs young and old post-collapse.

We thank the reviewer for pointing out the confusion. We have revised the figure markers as suggested and modified the main text accordingly.

• Young; pre-collapse : star

• Young; post-collapse : diamond

• Old; pre-collapse : circle

• Old; post-collapse: triangle.

Figure 3 Panel C: Indicate that fixed demographic values here are pre-collapse. Also, as presented, it appears that there is a large group-size dependence that is not commented on. Previous results (Line 87 and Figure 2C) suggest a constant excavation area per ant of 11cm2/ant. Figure 3, panel C appears to suggest a group-size dependence. If these values are divided by group size, is excavated area per ant nearly constant across groups? How does the numerical value compare to the slope from Figure 2C?

We thank the reviewer for their insightful comments.

First, we would like to clarify that the area target of 11.1 (±1) cm²/ant, as described in Line 87, was obtained from the colony maturation experiments. In these experiments, we were unable to track the age of each individual ant, so the area target was calculated by normalizing the total excavated area by the number of ants.

We normalized the excavated area by the group size for both young and old colonies as suggested, and found that the area per ant was not significantly different across the group sizes (see new SI Fig. 5A). This indicates that the excavated area per ant remains relatively constant within each demographic group. Moreover, this shows that the total excavated area is proportional to group size, in agreement with previous works (24, 29, and 30). 

We have explicitly described the above information in the line numbers: 142-146

Regarding the slope comparisons, the slope of Figure 2C (10.71), from the colony maturation experiments, is the largest, followed by the area per ant from the short-term young (8.79 ± 0.98) cm²/ant, and short-term old experiments (5.16 ± 0.44) cm²/ant.

Lines 128-129: "...younger ants aim to approach a higher target area". Seems hard to know what they "aim" to do... rephrase to report what they are observed to do.

We thank the reviewer for the comment. The sentence is rephrased as suggested (line numbers : 158-161).

“In the previous sections, we showed that in fixed-demographics experiments, younger ants excavated a significantly larger nest area compared to older ants (Fig. 3. C).  This difference emerged despite similar temporal patterns in digging rates across age groups, with excavation activity peaking within the first 7 days before asymptotically decaying as nest expansion approached saturation (SI Fig. 8).”

Lines 133-141: The model description is not clear. Specifically, what parameters are ant-dependent? How does A relate to a?

We appreciate the reviewer's request for clarification. In our model:

(1) Equation 1 describes the change in the excavated area due to the digging activity of a single ant. Here, the variable 'a' represents the area excavated by one ant. This formulation allows us to capture the individual digging behavior and its impact on the excavation process.

(2) Equation 2 extends this concept to the total area excavated in the nest, denoted by 'A'. Specifically, 'A' is the sum of the areas excavated by all ants present in the nest. In other words, it aggregates the individual contributions of each ant, linking the microscopic digging behavior to the macroscopic excavation dynamics.

Therefore, the relationship between 'a' and 'A' is as follows:

●     'a' = Area excavated by a single ant.

●     'A' = ∑ 'a' (Summed over all ants in the nest).

We have explicitly mentioned this in the line numbers “ 161-179”, and describe the model assumptions and parameters in detail.

Figure 4:

Figure 4, Panel A: The equation quoted in the caption does not match the data in the figure. The equation has a positive slope and negative intercept, while the figure has a negative slope and a positive intercept. Please provide the correct equation and bounds on fit parameters.

We thank the reviewer for spotting this typing mistake.

The equation was already updated in the reviewed preprint published online. The correct equation and the fit bound are provided in the figure caption.

“Target areas decrease linearly with the ant age (y = −0.032x + 11.22 , 95 % CI (Intercept : (-0.035,-0.027), Slope : (10.53,11.91)), R2 = 0.96 ).”

Figure 4, Panel A: There seem to be three "fixed target area per ant values" in the paper: around 11cm2/ant (line 87), 11.6 cm2/ant (SI Figure 2), and linearly dependent value from fit to Figure 4A. The distinctions between these values and their significance are hard to keep track of. Can the authors add a discussion somewhere that helps the reader better understand? Is there a way to connect/rationalize/explain these different values in terms of demographics?

We thank the reviewer for the suggestion.We have added a paragraph in the discussion (line numbers : 270-277) describing the area targets.

“In our colony maturation experiments, we found that area per ant was highest when the workers were youngest, with values around 11.1–11.6 (±1–1.15). This aligns with observations from naturally growing nests, where newly eclosed ants dominate the population and nest volumes are relatively large. Supporting this, fixed-demographics experiments showed that the area excavated per ant declines linearly with worker age, indicating that the youngest ants contribute most to excavation. Notably, the target area we fit for the age-independent model (11.6 ± 1.15) closely matches the extrapolated value for very young workers (Fig. 4. A), reinforcing the idea that young ants are the primary excavators during early colony growth. In contrast, during events like collapses or displacement, when space is urgently needed, ants of all ages participate in excavation.”

Figure 4, Panel A: What are various symbols and colors for data with error bars? If consistent with Figure 3, then this panel and subsequent model confound two factors: (1) the age dependence and (2) the behavioral differences pre- and post-collapse (structures are different pre-and post-collapse, according to SI Figure 6; line 120: "...colonies ceased digging when they recovered 93{plus minus}3% of the area lost by the manual collapse..."; lines 201-202: "We find significant quantitative and qualitative differences between nests constructed within this natural context and nests constructed in the context of an emergency") and behavior is different (according to SI Figure 7 and line 119: "...all ants dig after collapse...")). Therefore, without further supporting evidence, it does not seem that these data should be used to fit a single line that defines a model parameter a_age for each ant in equation 2.

The symbols are the area per ant quantified from the fixed demographics of young, and old experiments. The symbols show the following;

A.  Star - Young, pre-collapse

B.  Diamond - Young, post-collapse 

C.  Circle - Old, pre-collapse

D.  Triangle - Old, post-collapse.

The details are clearly described in the figure caption. 

We apologize to the reviewer for the confusion. We argue that the data can be fit by a single line to quantify the parameter ‘a_age’ as follows. 

A. All data presented in Figure 4A were obtained from the same fixed-demographics experiments (containing only young and old ants) under experimental collapse conditions, pre- and post-collapse. These results, therefore, exclusively reflect emergency nest-building behaviors during emergency scenarios and do not include any observations from natural colony maturation processes.

B. Age-dependent excavation differences: As correctly noted by the reviewer, the observed difference in excavated area before versus after collapse reflects the natural aging of ants in our experimental colonies. While colonies recovered >90% of lost area post-collapse, the residual variation was not negligible—instead, it systematically correlated with colony age structure. By tracking colonies across this demographic transition, we obtained additional data points spanning a broader developmental spectrum. This extended range strengthened our ability to detect and quantify the linear relationship between worker age and excavation output.

C.The quoted sentence (lines 201-202, submitted version) refers to comparisons across all three experimental cases: (1) fixed-demographics young ants, (2) fixed-demographics old ants, and (3) the natural scenario (mixed-age colonies). Importantly, these comparisons are based on pre-collapse steady-state excavation areas, ensuring a consistent baseline across treatments. We highlight quantitative and qualitative differences between these distinct experimental groups, not between pre- and post-collapse phases within the same treatment. The pre- and post-collapse data within fixed-demographics groups were analyzed separately to avoid conflating aging effects with emergency responses.

To avoid confusion, the whole paragraph in the discussion (line numbers : 253-260) is rephrased.

In lines 201-202; “We find significant quantitative and qualitative differences between nests constructed within this natural context and nests constructed in the context of an emergency”. 

Here, by natural context, we mean the nests excavated in the colony maturation experiments. We believe that it could have been confusing, and the sentence is modified as answered for the previous question. 

Figure 4, Panel B: This uses the model with a_age determined by from Figure 4A and the life table (as shown in the supplemental), whereas the supplemental Figure SI 8 uses the fixed blue line a_age value for the model, which comes from the colony maturation experiments. The age-independent model in the supplemental fits the data better, yet the authors claim the supplemental model cannot be applied to the data because of their experimentally determined age-dependent target area. Given the age-independent target area model fits better, additional evidence/justification is needed to support the choice of the model.

We agree with the reviewer that the age-independent model fits the data well. However, we believe that the fixed area target cannot be used to explain the excavation dynamics for the following reasons.

We make an important assumption in our model: that the ants rely on local cues and that individual ants can not distinguish between the fixed demographics and colony maturation experiments (line numbers : 161-166). Given this assumption, the ants cannot change their behavior between experiments, meaning the same model should fit all of our results. However, the fixed demographics experiments revealed a significant difference in the areas excavated by young vs. old cohorts, despite having the same group size. If the ants regulated the excavated area based on an age-independent constant density target model, then the excavated area in the fixed demographics of young and old colonies would have been similar. This discrepancy indicates that the target area per ant is not constant, as assumed in the age-independent density model (SI. Fig. 8). We emphasize that while the age-independent model provides a better fit for the excavated area in colony maturation experiments, the age-dependence of excavation is empirically supported by fixed-demographics experiments. Therefore, we implemented this age-dependence through a variable target area within the age-dependent model framework to explain excavation dynamics in the colony maturation experiments.

These details are explicitly mentioned in the main text (line numbers : 187 - 198)

Figure 4, Panel C: Is this plot entirely from the model, or are the data points measured from experiments? Please label this more clearly.

We apologize to the reviewer for the confusion.

The Figure 4C is based on the age-dependent digging model. We applied the model to population data from the long-term experiments (n = 22). By setting an age threshold of 56 days (since ants used in the short-term young experiment had an average age of 40 ± 16 days), we categorized the ants into young and old groups. We then quantified the area dug by the young ants, the queen, and the old ants in terms of the percentage of the total area excavated. We hypothesized that, because young ants have a lower digging threshold, they would perform the majority of the digging. We indeed confirm this in Figure 4C.

This information is added to the main text and described in detail (line numbers: 200 - 208).

Lines 162-165: "...Furthermore, we quantified the area dug by each ant in the normal colony growth experiment as estimated from the age-dependent model and found that all ants excavated more or less the same amount...". Figure 4D shows a distribution with significant values ranges from 1-16 cm2... how is this interpreted as "more or less the same amount" and what is the significance of this?

We apologise to the reviewer for the confusion.

We quantified the percentage contribution to the excavated area of each histogram bin (provided in the new SI table: 4), and found that the area excavated between 5 cm² and 13 cm² accounts for 73.76% of the total excavated area. This indicates that most ants dug within this range rather than exhibiting extreme variations. Additionally, the mean excavation amount is 7.84 cm², with a standard deviation of 3.44 cm², meaning that most values fall between 4.4 cm² and 11.28 cm², which aligns well with the 5–13 cm² range. Since the majority of the excavation is concentrated within this narrow interval, and the mean is well centered within it, this suggests that ants excavated more or less the same amount, rather than forming distinct groups with highly different excavation behaviors.

We have modified the main text (line numbers: 209-216) to include these points.

The biological significance of this finding is that since all ants in the colony maturation experiments are born inside the nest, we hypothesize that they should excavate similar amounts. To test this, we quantified the area contribution of each ant over the entire duration of the experiment using the age-dependent digging model as described above and found that they indeed excavated more or less the same amount. From our analysis of fixed demographics experiments, we showed that the youngest ants excavate the largest area. Since the majority of the youngest ants participated in the colony maturation experiments, this further supports our hypothesis.

Figure 5.

Figure 5, Panels A-C: Please provide a scale bar. 

The scale bar is provided in the figure as suggested. The algorithm for the cutoffs for tunnel vs wide tunnels is described in detail in the section “Nest skeletonization, segmentation, and orientation.”

Figure 5, Panel E: Why does the chamber error bar for 5 ants go to zero?

In Figure 5, E, we plot the standard error, as described in the figure caption. In the experiments, the chamber area contributions were (0,0,39.94,0) respectively. The mean of the 4 numbers is 9.985, the standard deviation is 19.97, and the standard error is 9.985. So, the mean and the standard error are the same, so the lower error bar goes to zero, and the upper error bar goes to 19.97. This implies that in these experiments, the chamber area is often zero.

Figure 5, Panel I: Why are there no chambers for young colonies in I when they are in the histogram in E?

We apologize to the reviewer for the confusion. We initially missed adding the chamber orientation data of the young colonies to Panel I, but it has now been included.

Line 212: "...densities of ants never become too high...". What is too high? Is there some connection to biological or physical constraints?

Under normal growth conditions, nest volume is kept proportional to the number of ants, ensuring that the density remains within a specific range. This prevents overcrowding, which could otherwise lead to excessively high densities.

Yes, we believe there is likely a connection to both biological and physical constraints. The proportional relationship between nest volume and the number of ants is likely driven by factors such as:

(1) Biological Constraints:

Ant Colony Size: Ants typically adjust their behavior and social structure to maintain an optimal population size relative to available resources and space.Overcrowding could lead to potentially a breakdown in colony function.

Colony Health: High densities can lead to faster epidemic spread, leading to negative effects on reproduction, foraging efficiency, and overall colony health. By maintaining density within a specific range, the colony can thrive without these adverse effects.

(2) Physical Constraints:

Spatial Limitations: The physical space within the nest limits how many ants can occupy it before space becomes constrained. The nest’s structure and size must physically accommodate the ants, and the volume must be large enough to prevent overcrowding, and efficient resource distribution.

Lines 272 and 302: How often were photos taken? These two statements seem to suggest different data collection rates.

As stated in line 272, photos were taken every 1 to 3 days. During each photo session, four photos were taken, with each photo separated by 2 seconds, as mentioned in line 302. To avoid confusion, we rephrased the sentence (line numbers: 359-361).

“We photographed the nest development every 1-3 days. During each photography session, four pictures of the nest were taken, with a 2-second interval between each.”

Reviewer #2 (Recommendations for the authors):

Some more minor points/questions/clarifications:

This might be pedantic, but I don't think the nest serves as the skeleton of the superorganism, while it does change and grow, the analogy becomes weak beyond that point. The skeleton serves to protect the internal organs of the organism, facilitates movement and muscle attachment, and creates new blood cells. I would be more comfortable with a statement that the nest can grow or shrink according to need.

We sincerely thank the reviewer for their time and effort in providing a detailed review and assessment of our manuscript. A point-by-point response to the comments is provided below.

The analogy of treating a nest structure to the skeleton of a superorganism was based on the following points;

(a) Protection: A nest protects the colony on a collective scale. This is analogous to protecting "organs" by a skeletal framework.

(b) Organization and Division of Space: The skeletal structure organizes the body's internal layout, just as nest structures are organized into various spatial compartments for various colony functions, with specific regions designated for brood chambers, food storage, and waste disposal.

Thus, we believe that the analogy can still be valid in a metaphorical way.

Does this statement need justification with a citation, or is that information contained in the subsequent clause? "However, for more complex structures where ants congregate in specific chambers, workers are less likely to assess the overall nest density." The idea that workers do (or do not) assess overall density touches on many issues, including that of perfect information and adaptive responses, that it seems it needs to be well founded in previous work to be stated in such unequivocal terms.

We thank the reviewer for this comment. The references for this argument are provided in the next sentence. We have now moved these references to the relevant sentence (reference number: 24, 29,30; line number : 30-31 ) 

Can you give some more information on this statement? "Experiments were terminated either when the queen died or when she became irreversibly trapped after a structural collapse." Why was this collapse irreversible and therefore unlike treatment 2? Did the queen die in these instances? Was this event more likely than in natural colonies? And if so, was there something inherently different about your experiments that limit interpretation under natural conditions (e.g. the narrow nature of the observation setup? The consistency of the sand?)

Our nest excavation experiments were terminated under two primary scenarios: (1) the queen died of natural causes, reflecting the baseline mortality expected when queens are brought into laboratory conditions, or (2) the nest experienced a structural collapse that left the queen irreversibly trapped. The second scenario is further elaborated below:

Irreversible Collapses: These collapses were classified as irreversible because the queen could not be rescued alive. This occurred when the structural stability of the nest failed, burying the queen in a manner that prevented recovery. In some cases, the collapse resulted in the queen's immediate death, while in others, she was trapped beyond reach, and any rescue attempt risked further structural damage.

Collapse and Experimental Context: These collapses were not uniquely associated with natural colonies or fixed-demographic experiments; rather, they occurred across various experimental setups.

The sentence is modified as below to improve clarity (line numbers : 70-72 ).

“In all instances where a collapse resulted in the queen's death or her being irreversibly trapped in the nest, the experiment was excluded from analysis starting from the point of the collapse, as such events did not reflect normal colony dynamics.”

I want to make sure I understand the following statement: "Moreover, the area excavated by the young cohorts was similar to that excavated by naturally maturing colonies at the point in which they reached the same population size (Tukey's HSD; group size: 5; p = 0.61, group size: 10; p = 0.46, group size: 15; p = 0.20)." Do I have it right that this means a group of (e.g. 10) young ants excavates an area similar to that of a group of 10 naturally maturing ants at the same age as the young ants?

Yes, the interpretation provided is correct. We apologize to the reviewer for the confusion. We have rephrased the sentence for better readability (line numbers : 146-148).

“Furthermore, the area excavated by the young cohorts was comparable to that excavated by naturally maturing colonies when they reached the same population size (Tukey's HSD; group size: 5, p = 0.61; group size: 10, p = 0.46; group size: 15, p = 0.20)”

How old do ants get? Is the 'old' demographic (~200 days) meaningfully old in the context of the overall worker lifespan? While the results certainly demonstrate there is an age effect, I would like to understand how rapid this is in terms of overall lifespan.

The lifespan of ants, including both queens and workers, varies significantly based on caste, species, and environmental conditions.

(1) Queen Longevity: From the literature, Camponotus fellah queens can live up to 20 years, with one documented case reaching 26 years. This remarkable longevity underscores the queen's central role in maintaining the colony.

(2) Worker Longevity: In contrast to queens, the lifespan of workers is much shorter.

However, specific data on worker longevity in Camponotus fellah colonies are lacking. Studies on other Camponotus species (50, 82) suggest that workers can live for several months depending on environmental conditions, colony health, and caste-specific roles (e.g., minor vs. major workers).

(3) Laboratory vs. Natural Conditions: Worker longevity is highly variable between laboratory and natural conditions

Therefore, in the context of the old worker lifespan in our experiments of, ~200 days (roughly 6–7 months) we strongly believe that the worker lifespan used in our experiments represents a substantial portion of a worker's expected life. While exact figures for C. fellah workers are unavailable, inferences from related species suggest that workers nearing 200 days are approaching the latter stages of their lifespan, making them meaningfully "old."

These details are added to the main text (line numbers : 124 - 127) and to the discussion (line numbers : 278-282)

Reviewer #3 (Recommendations for the authors):

We sincerely thank the reviewer for their time and effort in providing a detailed review and assessment of our manuscript. A point-by-point response to the comments is provided below.

L10: "fixed demographics": I find this term unclear, what does it mean, it should specify if the groups are with or without a queen.

We thank the reviewer for the comment. The sentence is modified in the abstract, and definitions are later added in detail in the introduction (line numbers : 8-10) and the Materials and Methods section (Fixed demographics colonies). 

“We experimentally compared nest excavation in colonies seeded from a single mated queen and allowed to grow for six months to excavation triggered by a catastrophic event in colonies with fixed demographics, where the age of each individual worker, including the queen, is known”.

The details of the “fixed demographics” treatments were explained in the later portion of the text (line numbers: 58-61).

L36: I think it is documented that younger individuals are the ones who involved in nest construction in many species.

Previous studies on nest construction were predominantly performed on mature colonies of specific age demographics or rather mixed demographics, where age was not considered as a factor influencing nest construction. Some studies have speculated that young ants could be the most probable ones to dig, but this has not been experimentally verified to the best of our knowledge.

L50: I do not think the colony should be called mature after only 6 months, given that colonies reach thousands of workers.

The sentence is changed as suggested (line numbers : 56-57).

“The "Colony-Maturation" experiment observed the development of colonies up to six months, starting from a single fertile queen and progressing to colonies with established worker populations.” 

L60: Where was the queen introduced? It is specified in the Methods but a word here would be helpful.

The detail is added as suggested (line numbers : 68-69).

“We initiated colony maturation experiments by introducing a single mated queen and several brood items (n = 5, across all experiments) at random positions on the soil layer of the nest.”

L106: Young vs Old workers 40 vs 171 days. Maybe cite a reference or provide a reason for the selection of those ages?

Previous studies have shown that the Camponotus fellah queens can live up to 20 years, with one documented case reaching 26 years (50). To the best of our knowledge, specific data on worker longevity in Camponotus fellah colonies in natural conditions are lacking. Lab studies on Camponotus fellah (82) and other Camponotus species (50) suggest that workers can live for several months depending on environmental conditions, colony health, and caste-specific roles (e.g., minor vs. major workers). 

We intentionally selected workers from two distinct age groups: younger ants (40 ± 16 days old) and older ants (171.56 ± 20 days old). These ages represent functionally different life stages - the younger group had completed about 25% of their expected lifespan at the start of the experiment, while the older group had lived through most of theirs (50, 82). This 4-fold age difference allowed us to compare excavation behaviors across fundamentally different phases of adult life.

Our experiments lasted for 60-90 days, during which all participating workers continued to age. To ensure all ants remained alive throughout the experiments, and given the constraints of the experimental timeline, we selected young and old workers within the specified age range. 

These details are added to the main text (line numbers :  124 -127), and the discussion (line numbers  : 278-282)

L122-123: But usually ants can vary highly in their behaviours. Can the authors comment on their choice to consider an average, implying that all ants of the same age had the same digging rates?

We thank the reviewer for the comment.

In our experiments, we could not track each worker's activity over time. As described in the methods, we took snapshots of the nest structure over days and recorded the population size of the nest. Thus, we could not capture the activity of single ants in the nest as described in the response to major comments in the reviewed preprint.

We agree that individual tracking of ants within our experimental setup would have been the ideal approach. Then, we could have taken the inter-individual variability of the digging activity into account. However, we were limited to doing so by the technical and practical limitations of the setup, such as; 

(a) Continuous tracking of ants in our nests would have required a camera to be positioned at all times in front of the nest, which necessitates a light background. Since Camponotus fellah ants are subterranean, we aimed to allow them to perform nest excavation in conditions as close to their natural dark environment as possible. Additionally, implementing such a system in front of each nest would have reduced the sample sizes for our treatments.

(b)The experimental duration of our colony maturation and fixed demographics experiments extended for up to six months (unprecedented durations in these kinds of measurements). These naturally limited our ability to conduct individual tracking while maintaining the identity of each ant based on the current design.

To clarify this, we have added the following to the discussion (line numbers: 286-292).

“Previous studies have demonstrated both homogeneous and heterogeneous workload distribution, with varying digging rates among ants (24,29,30,35). Studies showing heterogeneous workload distribution relied on continuous individual tracking of ants to quantify digging rates (35). However, this approach was not feasible in our current design due to the experimental durations of both our colony maturation and fixed demographics experiments. Additionally, sample size requirements naturally limited our ability to conduct continuous individual tracking during nest construction in our study.”

L171: A line on how the nest structure was acquired and data extracted would be welcome here.

The algorithm for the nest structure segmentation, data extraction, and analysis is added in detail to the SI section: Nest skeletonization, segmentation, and orientation. The line is modified (line numbers : 221-224) in the main text as suggested.

“We compared nest architectures by segmenting raw nest images into chambers and tunnels (see SI Section: Nest Skeletonization, Segmentation, and Orientation). Chambers were identified as flat, horizontal structures, while tunnels were narrower and more vertical in orientation (see SI Fig. 9, SI Section: Nest Skeletonization, Segmentation, and Orientation)”.  

Figure 3: Where does the data of the mean in panel C come from: is it the mean of the first 30 days, before the collapse? How is it comparable with the rest?

We apologize to the reviewer for the confusion.

In panel C, the mean values (solid stars and circles) for fixed-demography colonies (young/old groups) represent pre-collapse excavation areas. For colony maturation experiments (where no collapses were induced), we instead plot the mean saturated excavation area for each group size. This allows direct comparison of mean excavated areas across experimental conditions at equivalent colony sizes.

To improve readability, the following sentences are added to the main text (line numbers : 139 - 146 ) 

“We compared the saturated excavation areas (pre-collapse) from fixed-demographics experiments (young and old groups) with those from colony maturation experiments of the same colony sizes (Fig. 3C). We find that, for a given age cohort (young or old), the saturation areas increase linearly with the colony size (GLMM, F(35,37); p < 0.0001) (Fig. 3 C, SI. Fig 7 A). The observed proportional scaling between excavated area and group size aligns with previous studies, even though those studies did not explicitly account for age demographics (24, 29, 30). After normalizing the pre-collapse excavated area by group size for both young and old colonies, we found no significant difference in area per ant across group sizes (SI Fig. 5. A). This indicates that the excavated area per ant remains relatively constant within each demographic group”.

L209-210: I would be more parsimonious in saying that the results presented prove that the target area decreases with age, as the individual behaviour of the ants was not monitored. Suggestion: rephrase to "the target of the group decreases with age".

The sentence is rephrased as suggested (line numbers : 265-266).

“Our results reveal that this target area of the group decreases linearly with age, such that young ants are more sensitive to shortages in space.”

L246: Are C.fellah colonies really found with such few workers?

Previous studies have speculated that mature Camponotus fellah colonies are a monogynous species typically founded by a single queen following nuptial flights (50,51,82), and can range from tens to thousands of workers. However, during the founding stage (as in our experiments), colonies naturally pass through smaller developmental sizes comparable to the matured colonies.

Author response:

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

Reviewer #1 (Public review): 

Summary:  

The Szczupak lab published a very interesting paper in 2012 (Rodriquez et al. J Neurophysiol 107:1917-1924) on the effects of the segmentally-distributed non-spiking (NS) cell on crawl-related motoneurons. As far as I can tell, the working model presented in 2012, for how the non-spiking (NS) cell impacts the crawling motor pattern, is the same functional model presented in this new paper. Unfortunately, the Discussion does not address any of the findings in the previous paper or cite them in the context of NS alterations of fictive crawling. Aside from different-looking figures and some new analyses, the results and conclusions are the same. 

Reviewers #1 and #2 called our attention to our failure to cite the Rodriguez et al. 2012 article in the context of the main goal of the present work. We do now explain how the present study is framed by the published work. See lines 74-79.

In Rodriguez et al. 2012, we hypothesized that the inhibitory signals onto NS were originated in the motoneuron firing. We now cite this reference in line 104. In the current manuscript we further investigated the connection between the inhibitory signals onto NS and the motoneuron activity (Figure 2) and proved that the hypothesis was wrong. Thus, the model presented here differs from the one proposed in Rodriguez et al. 2012.

In Rodriguez et al. 2012, we speculated that the inhibitory signals received by NS were transmitted to the motoneurons, but an important control was missing in that study. In the current study depolarization of NS during crawling is tested against a control series that allows to properly examine the hypothesis (lines 138-147). But, most important, because NS is so widely connected with the layer of motoneurons it was necessary to test the effect on other motoneurons during the fictive crawling cycle. We now explain this rationale in lines 249-257.

Strengths: 

The figures are well illustrated. 

Weaknesses:  

The paper is a mix of what appears to be two different studies and abruptly switches gears to examine how closely the crawl patterning is in the intact animal as compared to the fictive crawl patterning in the intact animal. Unfortunately, previous studies in other labs are not cited even though identical results have been obtained and similar conclusions were made. Thus, the novelty of the results is missing for those who are familiar with the leech preparation. The lack of appropriate citations and discussion of previous studies also deprives the scientific community of fully comprehending the impact of the data presented and the science it was built upon.  

The main aim of the manuscript is to learn the role of premotor NS neurons in the crawling motor pattern studied using spike sorting in extracellular nerve recordings. This readout allows to  simultaneously monitor a larger number of units  than in any previous study. This approach aims to determine whether and how a recurrent inhibitory peripheral circuit is involved in coordinating or modulating the rhythmic motor pattern.

Our rationale was that the known effect of NS on one particular motoneuron (DE-3) may have overlooked a more general effect on crawling (lines 253-257). Moreover, we wanted to investigate whether this effect was due to the recurrent inhibitory circuit or if other elements were involved, and to study whether the modulation was mediated by the recurrent synapse between NS and the motoneurons.

In the context of this aim we studied the rhythmic activity of cell DE-3, together with motoneurons that fire in-phase and anti-phase, in isolated ganglia (Figure 4). To reveal the effect of NS manipulation we applied a quantitative analysis that showed the phase-specific effect of NS (Figure 6). 

Given that this is the first study using a spike sorting algorithm to detect and describe the activity of motoneurons in nerve recordings we found it reasonable to compare these results with an in vivo study; thus, providing information to the general reader, that supports the correspondence between the ex vivo and the in vivo patterns.

(1) Results, Lines 167-170: "While multiple extracellular recordings have been performed previously (Eisenhart et al., 2000), these results present the first quantitative analysis of motor units activated throughout the crawling cycle. The In-Phase units are expected to control the contraction stage by exciting or inhibiting the longitudinal or circular muscles, respectively, and the Anti-Phase units to control the elongation stage by exciting or inhibiting the circular or longitudinal muscles, respectively."  

Reviewer: The first line above is misleading. The study by Puhl and Mesce (2008, J. Neurosci, 28:4192- 420) contains a comprehensive analysis of the motoneurons active during fictive crawling with the aim of characterizing their roles and phase relationships and solidifying the idea that the oscillator for crawling resides in a single ganglion. Intracellular recordings from a number of key crawl-related motoneurons were made in combination with extracellular recordings of motoneuron DE-3, a key monitor of crawling. In their paper, it was shown that motoneurons AE, VE-4, DI-1, VI-2, and CV were all correlated with crawl activity, and fired repeatedly either in phase or out-of-phase with DE-3. They were shown to be either excitatory or inhibitory. At a minimum, the above paper should be cited. 

The sentence in the submitted manuscript explicitly refers to the quantitative analysis of extracellular recordings, but we recognize that it may lead to confusion. We have now added a clarification (lines 197-199). 

The article by Puhl and Mesce 2008 shows very nice intracellular recordings of the AE, CV, VE-4, DE-3, DI-1, and Vi-2, accompanied by extracellular recordings of DE-3 in the DP nerve. In all cases, there is only one intracellular recording paired with the DP nerve recording.

While it is possible to perform up to 3-4 simultaneous intracellular recordings, these are technically challenging, and more so when the recordings have to last 10-20 minutes. Due to this difficulty, and because our objective was to record multiple units simultaneously in order to comprehensively describe the different crawling stages, we implemented the spike sorting analysis on multiple extracellular recordings. This approach enabled us to reliably obtain multiple units per experiment and thus execute a quantitative analysis of the activity of each identified unit.

The article by Puhl and Mesce 2008 mentions several quantitative aspects of the neurons that fire in-phase or out-of-phase with DE-3, but, as far as we understand, there is no figure that summarizes activity levels and span in the way Figures 4 and 6 do in the current manuscript. To the best of our knowledge, no previous work renders this information.

It is very important for us to emphasize that the work by Puhl and Mesce was seminal for our research. We cited it four times in the original manuscript and 10 times in the present version. But, like any important discovery, it sets the ground for further work that can refine certain measurements that in the original discovery were not central.

This is why we believe that the cited sentence in our manuscript is not misleading.  However, to comply with the requirement of Reviewer #1, we added a sentence preceding the mentioned paragraph (lines 185-187) that acknowledges the description made using intracellular recordings, and explains the need for implementing the approach we chose.

The submitted paper would be strengthened if some of these previously identified motoneurons were again recorded with intracellular electrodes and concomitant NS cell stimulation. The power of the leech preparation is that cells can be identified as individuals with dual somatic (intracellular) and axonal recordings (extracellular). 

Most of the motoneurons mentioned by Reviewer #1 are located on the opposite side (dorsal) of the ganglion to NS (ventral), and therefore, simultaneous intracellular recordings in the context of fictive crawling are challenging.

In the publication of Rodriguez et al. 2009, Mariano Rodriguez did manage to record NS from the dorsal side together with DE-3 and MN-L (!) and this led to the discovery that these motoneurons are electrically coupled, but the recurrent inhibitory circuit masks this interaction. Repeating this type of experiments during crawling, which requires stable recordings for around 15 minutes, is not a reasonable experimental setting.

Rodriguez et al. 2012 shows intracellular recordings of motoneurons AE and CV during crawling in conjunction with NS, and their activity presented the expected correlation. 

The shortfall of this aspect of the study (Figure 5) is that the extracellular units have not been identified here. 

The Reviewer is right in that the extracellular units have not been identified in terms of cell identity. As we explained earlier, most motoneurons are on the opposite side (ventral/dorsal) of the ganglion relative to NS. 

However, we do characterize the units in terms of the nerve through which they project to the periphery and their activity phase. In lines 345-349 we use this information and, based on published work, we propose possible cellular identities of the different units.

In xfact, these units might not even be motoneurons. 

We are surprised by this comment. The classical work of Ort and collaborators (1974) showed that spikes detected in extracellular nerve recordings were emitted by specific motoneurons, and several previous publications have validated extracellular nerve recordings as a means to study fictive motor patterns (Wittenberg & Kristan 1992, Shaw & Kristan 1997, Eisenhart et al. 2000).

For further reassurance, we only took in consideration units whose activity was locked to DE3; any non-rhythmical activity was filtered out (see lines 433-435). 

They could represent activity from the centrally located sensory neurons, dopamine-modulated afferent neurons or peripherally projecting modulatory neurons. 

Peripheral nerves also contain axons from sensory neurons. However, in a previous article, we studied the activity of mechanosensory neurons (Alonso et al. 2020) and showed that they remain silent during crawling. Moreover, the low-threshold T sensory neurons are inhibited in phase with DE-3 bursts and NS IPSPs (Kearney et al. 2022). Alonso et al. 2000 showed that spiking activity of T cells affects the crawling motor pattern, revealing the relevance of keeping them silent.

What does the Reviewer mean by “dopamine-modulated afferents”? We are not aware of this category of leech neurons.

The neuromodulatory Rz neurons project peripherally through the recorded nerves, but intracellular recordings of these neurons from our lab show no rhythmic activity in those cells during dopamine-induced crawling.

Essentially, they may not have much to do with the crawl motor pattern at all.

Does the Reviewer consider that neurons engaged in a coherent rhythmic firing could be unrelated to the pattern? As indicated above, the units reported in our manuscript were selected because dopamine evoked their rhythmic activity, locked to DE-3. 

Does the Reviewer consider that dopamine could evoke spurious neuronal activity?

(2) Results Lines 206-210: "with the elongation and contraction stages of in vivo behavior. However the isometric stages displayed in vivo have no obvious counterpart in the electrophysiological recordings. It is important to consider that the rhythmic movement of successive segments along the antero-posterior axis of the animal requires a delay signal that allows the appropriate propagation of the metachronal wave, and this signal is probably absent in the isolated ganglion." 

Reviewer: The so-called isometric stages, indeed, have an electrophysiological counterpart due in part to the overlapping activities across segments. This submitted paper would be considerably strengthened if it referred to the body of work that has examined how the individual crawl oscillators operate in a fully intact nerve cord, excised from the body but with all the ganglia (and cephalic ganglion) attached. Puhl and Mesce 2010 (J. Neurosci 30: 2373-2383) and Puhl et al. 2012 (J. Neurosci, 32:17646 -17657) have shown that "appropriate propagation of the metachronal wave" requires the brain, especially cell R3b-1. They also show that the long-distance projecting cell R3b-1 synapses with the CV motoneuron, providing rhythmic excitatory input to it.  

We would like to draw the Reviewer’s attention to the fact that Puhl and Mesce 2008, 2010 and Puhl et al. 2012 characterized crawling in intact (or nearly intact) animals considering the whole body. In our in vivo analysis, we studied the changes in length of the whole animal and of sections demarcated by the drawn points, as described in the Materials and Methods/Behavioral

Experiments. Because of this different analysis, we defined “isometric” stages as those in which a given section of the animal does not change its length. We now clarify this (line 230).

In the paragraph cited by the Reviewer, we intended to state that, in the context of our study, the intersegmental lag caused by the coordinating mechanisms has no counterpart “in the electrophysiological recordings of motoneurons in the isolated ganglia”. We have now completed this idea with the expression underlined in the previous sentence (line 231).

As the Reviewer indicates, in the intact nerve cord the behavioral isometric stages correspond to the “waiting time” between segments. We did refer to the metachronal order but did not cite the articles by Puhl and Mesce 2010 and Puhl et al. 2012; we now do so (lines 234).

For this and other reasons, the paper would be much more informative and exciting if the impacts of the NS cell were studied in a fully intact nerve cord. Those studies have never been done, and it would be exciting to see how and if the effects of NS cell manipulation deviated from those in the single ganglion.  

The Reviewer may consider that a systematic analysis of multiple nerves in several ganglia along the whole nerve cord would have been a different enterprise than the one we carried out. The Reviewer is right in recognizing the interest of such study, but in our opinion, the value of the present work lies in presenting a thorough quantitative analysis of multiple nerves to demonstrate its usefulness for the study of the network underlying leech crawling. In this manuscript, we used it to analyze the role of the premotor NS neuron. Without the recording of units firing in-phase and out-ofphase with DE-3, we would have been unable to assess the span of NS effects.

(3) Discussion Lines 322-324. "The absence of descending brain signals and/or peripheral signals are assumed as important factors in determining the cycle period and the sequence at which the different behavioral stages take place." 

Reviewer: The authors could strengthen their paper by including a more complete picture of what is known about the control of crawling. For example, Puhl et al. 2012 (J Neurosci, 32:17646-17657) demonstrated that the descending brain neuron R3b-1 plays a major role in establishing the crawlcycle frequency. With increased R3b-1 cell stimulation, DE-3 periods substantially shortened throughout the entire nerve cord. Thus, the importance of descending brain inputs should not be merely assumed; empirical evidence exists.  

We now strengthen the concept using “known descending brain signals” (line 358) and cite Puhl et al. 2012. We believe that extending the discussion to cell R3b-1 does not contribute meaningfully to the focus of this manuscript.

(4) Discussion Lines 325-327: "the sequence of events, and the proportion of the active cycle dedicated to elongation and contraction were remarkably similar in both experimental settings. This suggests that the network activated in the isolated ganglion is the one underlying the motor behavior." 

Reviewer: The results and conclusions drawn in the current manuscript mirror those previously reported by Puhl and Mesce (2008, J. Neurosci, 28:4192- 420) who first demonstrated that the essential pattern-generating elements for leech crawling were contained in each of the segmental ganglia comprising the nerve cord. Furthermore, the authors showed that the duty cycle of DE-3, in a single ganglion treated with dopamine, was statistically indistinguishable from the DE-3 duty cycle measured in an intact nerve cord showing spontaneous fictive crawling, in an intact nerve cord induced to crawl via dopamine, and in the intact behaving animal. What was statistically significant, however, was that the DE-3 burst period was greatly reduced in the intact animal (i.e., a higher crawl frequency), which was replicated in the submitted paper.  

There is no doubt that the article by Puhl and Mesce 2008 is seminal to the work we present here. The Reviewer seems to suggest that we do not recognize the value of this work. The contrary is true, all our related papers cite this important breakthrough. We cite the paper very early in the article in the Introduction (see lines 51 and 52-53). Likely, we would like the Reviewer to recognize the novelty of the current report. To clarify what has been shown and what is new in our manuscript, considerer the following:

i. Figures 1-6 in Puhl and Mesce 2008 provide representative intracellular recordings that describe neurons that fire in phase and out of phase relative to DE-3. Some general measurements are given in the text, but none of these figures quantify the relative activity of neurons that fire in different stages; only DE-3 activity was quantified. A quantitative description of multiple units active in phase and out of phase with DE-3 is presented here for the first time, are we wrong? This quantification is particularly relevant when assessing how a treatment affects the function of the circuit.

ii. Regarding the cycle period, we referred to the work from the Kristan lab, which reported this value long before the requested reference. We now cite Puhl and Mesce 2008 in lines 222 regarding in vivo measurements, and in line 221 regarding isolated ganglia.

iii. Regarding the duty cycle: 

Puhl and Mesce 2008 measured the duty cycle of DE-3 in three configurations: a. spontaneous whole cord, b. DA-mediated whole cord and c. DA mediated single ganglion crawling. However, it does not report the duty cycle of neurons out-of-phase with DE-3. Our current manuscript carried out this analysis. One could argue that the silence between DE-3 bursts captures that value, but this is a speculation that needed a proper measure.

Puhl and Mesce 2008 does not indicate the duty cycle of the contraction and elongation stages in vivo. Our current manuscript does. 

Therefore, the sentence cited by the Reviewer refers to data presented in this manuscript, and not in any prior manuscript. It is true that Puhl and Mesce 2008 inspire the intuition that the sentence is true, but does not present the data that the current manuscript does.

Finally, our study focused only on the body sections corresponding to the same segmental range used in the ex vivo experiments, rather than the whole animal. The comparison was made only to validate that the duty cycles of neurons firing in phase and out of phase with DE-3 matched the dynamic stages in the studied sections of the leech (line 364).

In my opinion, the novelty of the results reported in the submitted manuscript is diminished in the light of previously published studies. At a minimum, the previous studies should be cited, and the authors should provide additional rationale for conducting their studies. They need to explain in the discussion how their approach provided additional insights into what has already been reported.  

Throughout our reply, we have provided a detailed explanation of the rationale and necessity behind each experiment. Following the Reviewer’s suggestion, we have rephrased the research objectives, included what is known from our previously published work, and highlighted the substantial new data contributed by the present study. See lines 80-85. 

Additionally, we further cite our published article in lines 93, 104, 138, 146 and 250. 

Reviewer #2 (Public review):  

The paper is well-written overall. The findings are clearly presented, and the data seems solid overall. I do have, however, a few major and some minor comments representing some concerns.

My major comments are below. 

(1) This may seem somewhat semantic, yet, it has implications on the way the data is presented and moreover on the conclusions drawn - a single ganglion cannot show fictive crawling. It can demonstrate rhythmic patterns of activity that may serve in the (fictive) crawling motor pattern. The latter is a result of the intrinsic within single-ganglion connectivity AND the inter-ganglia connections and interactions (coupling) among the sequential ganglia. It may be affected by both short-range and long-range connections (e.g., descending inputs) along the ganglia chain. 

Semantics is not a trivial issue in science communication. It entails metaphors that enter the bibliography as commonly used “shortcuts” to a complex concept that are adopted by a community of researchers. And yes, indeed, they can be misleading.

However, if recording the activity in an isolated ganglion shows that a wide group of motoneurons, that control known muscle movements, presents a rhythmic output that maintains the appropriate cycle period and phase relationships, the “shortcut” is incomplete but could be valid (Puhl and Mesce 2008). If we were to include the phase lag component, a single ganglion cannot generate the fictive motor output.

Because any new study builds knowledge on the basis of the cited bibliography, the way we name concepts is a sensitive point. Adopting the terminology used by previous publications (Puhl and Mesce 2008) seems important to allow readers to follow the development of knowledge. However, attending the observation made by Reviewer #2, we included a sentence clarifying that the concept “fictive crawling” does not include intersegmental connectivity (lines 54-57)

(2) The point above is even more critical where the authors set to compare the motor pattern in single ganglia with the intact animals. It would have made much more sense to add a description of the motor pattern of a chain of interconnected ganglia. The latter would be expected to better resemble the intact animal. Furthermore, this project would have benefitted from a three-way comparison (isolated ganglion-interconnected ganglia-intact animal.  

As we answered to Reviewer #1, the present manuscript does not intend to present a thorough study on how the activity in the isolated nervous system compares with the animal behavior. To do so we would have needed to perform a completely different set of experiments. To better define the relevance of our comparison with the in vivo experiments we rephrased the objective of the behavioral analysis (lines 197-199).

The main aim of the manuscript is to learn the role of premotor NS neurons in the crawling motor pattern studied using a readout (spike sorting in extracellular nerve recordings) that allows simultaneous screening of a larger number of units than in any previous study, in order to determine whether and how a recurrent inhibitory peripheral circuit is involved in coordinating or modulating the rhythmic motor pattern.

Our rationale was that the known effect of NS on one particular motoneuron (DE-3) may have overlooked a more general effect on crawling (lines 253-257). Moreover, we wanted to investigate whether this effect was due to the recurrent inhibitory circuit or if other elements were involved, and to study whether the modulation was mediated by the recurrent synapse between NS and the motoneurons.

In the context of this aim we studied the rhythmic activity of cell DE-3, together with motoneurons that fire in-phase and anti-phase, in isolated ganglia (Figure 4). To reveal the effect of NS manipulation we applied a quantitative analysis that showed the phase-specific effect of NS (Figure 6). 

Given that this is the first study using a spike sorting algorithm to detect and describe the activity of motoneurons in nerve recordings we found it reasonable to compare these results with an in vivo study; thus, providing information to the general reader, that supports the correspondence between the ex vivo and the in vivo patterns.

(3) Two previous studies by the same group are repeatedly mentioned (Rela and Szczupak, 2003; Rodriguez et al., 2009) and serve as a basis for the current work. The aim of one of these previous studies was to assess the role of the NS neurons in regulating the function of motor networks. The other (Rodriguez et al., 2009) reported on a neuron (the NS) that can regulate the crawling motor pattern. LL 71-74 of the current report presents the aim of this study as evaluating the role of the known connectivity of the premotor NS neuron in shaping the crawling motor pattern. The authors should make it very clear what indeed served as background knowledge, what exactly was known about the circuitry beforehand, and what is different and new in the current study. 

Rela and Szczupak 2003 and Rodriguez et al. 2009 analyze the interactions of motoneurons with NS. We believe that Reviewer #2 refers here to Rodriguez et al. 2012. A similar observation was made by Reviewer #1. Below, we copy the answer previously stated:

Following the Reviewer’s suggestion, we have rephrased the research objectives, included what is known from our previously published work, and highlighted the substantial new data contributed by the present study. See lines 80-85. 

Additionally, we further cite our published article in lines 93, 104, 138, 146 and 250. 

Reviewer #1 (Recommendations for the authors):  

Please edit for correct word usage. 

Reviewer #2 (Recommendations for the authors):  

Minor Concerns 

(1) LL33-36: These lines are somewhat vague and non-informative. Why is the functional organization of motor systems an open question? What are the mechanisms at the level of the nerve cord that are an open question? Maybe be more explicit? 

We did as suggested (lines 30-32).

(2) L62: The homology between the NS neurons and the vertebrate Renshaw cells is mentioned already in the Abstract and here again. While a reference is provided (citing the lead author of this current work), the reader would benefit from some further short words of explanation regarding the alleged homology. 

We included a description of Renshaw cell connectivity (lines 64-65).

(3) LL90-92: The NS recording in Figure 1 (similar to Figure 3 in Rodriguez et al.) demonstrates clear distinct IPSPs. Could these be correlated with DE-3 spikes? 

We investigated this correlation in detail and the answer is that there is no strictly a 1:1 DE-3 spike to IPSP correlation. NS receives inputs from other dorsal and ventral excitors of longitudinal muscles, and the NS trace is too “noisy” to reflect any short-term correlation. Originally we proposed that the NS IPSPs were due to the polysynaptic interaction between the MN and NS (Rodríguez et al. 2012). However, the present work demonstrates that the IPSPs in NS are caused by a source upstream from the MNs. 

(4) LL145-145: Do you mean - inhibitory signals FROM NS premotor neurons? Not clear. 

We see the confusion, and we rewrote the sentence (lines 164). We hope it is clearer now: “…inhibitory signals onto NS premotor neurons were transmitted to DE-3 motoneurons via rectifying electrical synapses and counteracted their excitatory drive during crawling, limiting their firing frequency.”

(5) LL153-154: Why isn't AA included in Figure 4A? 

Reading our original text, the Reviewer #1 is right in expecting to see the AA recording. We changed the sentence: “we performed extracellular recordings of DP along with AA and/or PP root nerves” (lines 171-172).

We dissected the three nerves but, unfortunately, we did not always obtain good recordings from the three of them.

(6) LL237-238: The statistical significance (B- antiphase) is not clear. Furthermore, with N of 7-8, I'm not sure the parametric tests utilized are appropriate. 

Regarding the Reviewer's concern about the tests, please note that all the assumptions made for each model were tested (see now Materials and Methods lines 466-467).The information on each model is provided in Supplementary Table 2 under the column 'Model, random effect,' which specifies whether a Linear Mixed Model (LMM) or a Generalized Linear Mixed Model (GLMM) was implemented. For GLMMs, the corresponding distribution and link function are also specified. For the analysis of Max bFF of Anti-Phase motor units, we found a significant interaction between epoch and treatment, indicating a difference between treatments. This is indicated on the left of the y-axis (##). In control experiments, all three comparisons (pre-test, pre-post, test-post) show significant differences in Max bFF: this variable decreased (slightly but significantly) along the subsequent epochs, suggesting a change over time. We now corrected the text to indicate that these changes were small (line 268). In contrast, Max bFF in depo experiments remained stable between pre-test and pre-post, but significantly decreased between the depo and post epochs. Thus, in our view the comparison between control and the test supports the conclusion that NS depolarization was limited to counteracting this decrease (lines 270-273). Supplementary Table 2 provides the significance and modeled estimated ratio for each comparison in the column for pairwise simple contrasts.

Thanks to this question, we realized that the nomenclature used in the table for the epochs (pre - depo - post) needed to be changed to pre - test - post, and we have now corrected it.

(7) LL240-241: I fail to see a difference from Control. 

For the Relative HW of In-Phase units, we also found a significant interaction between epoch and treatment, indicating a difference between treatments, as denoted to the left of the y-axis (#). Then, the significance of the comparisons across epochs within each treatment are shown in the figure (*). What is important to notice is that obtaining the same significance for each treatment does not imply identical results, but we failed to describe this in our original text and we do now in lines 275-279.

(8) LL244-245: I must admit that Table 2 is beyond me. Maybe add some detail or point out to the reader what is important (if at all). 

We have now clarified what each column of the tables indicates in the corresponding legends. 

Here, we also share an insight into how the experiments were designed and analyzed:

To account for possible temporal drifts of the variables during the recordings that could mask or confuse the results, we compared two experimental series: one in which NS was subjected to depolarizing current pulses (depo), and another series (ctrl) in which the neurons were not depolarized.

The statistical analysis was made using Linear Mixed Models (LMMs) or Generalized Linear Mixed Models (GLMMs). In these analyses treatments and epochs are used as explanatory variables to evaluate the interaction between these factors. These models allow us to determine whether changes in each variable across epochs differ depending on the treatment. For example, whether the variation in firing frequency from pre to test to post differs between control experiments and those in which NS was depolarized.

A significant interaction between treatment and epoch indicates that NS depolarization affected the variable. In such cases, we performed pairwise comparisons between epochs (pre-test, test-post, pre-post) within each treatment. In contrast, the absence of a significant interaction can result from two possibilities: either the variable did not change across epoch in either treatment, or a similar temporal drift occurred in both cases.

(9) LL245-256: Move this paragraph to the discussion. 

Because we introduced a rationale for the experiments described in Figure 6 (lines 282-284) the paragraph was mostly removed, but the part that supports the methodological approach was left.

(10)  LL259-260: see my second minor point above. This is explained in LL270-272 for the first time. 

We amended according to comment (2).

(11) Figures: The quantitative analysis shown in Figure 3B is very useful. Why isn't this type of analysis utilized for the comparisons shown in Figures 4 and 6? 

We chose different ways of plotting the data based on their nature. In Figure 3B, we present data from an identified neuron (DE-3) recorded in different experiments. In contrast, in Figure 6 we analyze data from neurons classified into the same group based on their activity during the fictive crawling cycle, but their individual identity was not ascertained. Therefore, we consider it important to plot the results for each unit individually, to assess the effect of temporal drift and NS depolarization.

(12) Figures: Figure 7 is meant to be compared to Figure 1C; the point being the addition of an inhibitory connection onto the NS neuron. Why are other details of the figure also different (different colored M)? 

While Figure 1C illustrates the known connection between NS and both DE-3 and CV motoneurons, Figure 7 shows the connections between NS and the different groups of motor units described in this study. The units are represented in the circuit using the same colors that identify them in Figures 4 and 6. Since the CV motoneuron was not recorded in this study, the circuit represents the AntiPhase neurons but does not identify them with CV. Figure 7 legend now clarifies what the colors represent, and Figure 1C has been updated to match the same color scheme.

En la resolución tendrá que realizarse una valoración racional de cada prueba, sin que sea válido el argumento que se valoró conforme a la íntima convicción del juzgador.

DOI: 10.1038/s41598-019-42630-y

Resource: RRID:BDSC_33613

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_33613

What is this?

DOI: 10.1038/s41598-019-42630-y

Resource: RRID:BDSC_58498

Curator: @scibot

SciCrunch record: RRID:BDSC_58498

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_1560

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_1560

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_992

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_992

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_27495

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_27495

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_35789

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_35789

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_58993

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_58993

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_458

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_458

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_25855

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_25855

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_7415

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_7415

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_55850

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_55850

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_57183

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_57183

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_7013

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_7013

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_58986

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_58986

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_26207

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_26207

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_7011

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_7011

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_36823

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_36823

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_8164

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_8164

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_51753

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_51753

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_14640

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_14640

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_41803

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_41803

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_62315

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_62315

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_48030

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_48030

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_58768

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_58768

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:BDSC_27390

Curator: @bdscstockkeepers

SciCrunch record: RRID:BDSC_27390

What is this?

DOI: 10.1038/s41467-019-09943-y

Resource: RRID:Addgene_73915

Curator: @scibot

SciCrunch record: RRID:Addgene_73915

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (BCRJ Cat# 0278, RRID:CVCL_0132)

Curator: @dhovakimyan1

SciCrunch record: RRID:CVCL_0132

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (Santa Cruz Biotechnology Cat# sc-17829, RRID:AB_626633)

Curator: @scibot

SciCrunch record: RRID:AB_626633

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (Thermo Fisher Scientific Cat# MA5-14520, RRID:AB_10979488)

Curator: @scibot

SciCrunch record: RRID:AB_10979488

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (MP Bio Cat# 0855484, RRID:AB_2334495)

Curator: @scibot

SciCrunch record: RRID:AB_2334495

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_2809831

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (Santa Cruz Biotechnology Cat# sc-74439, RRID:AB_1130818)

Curator: @scibot

SciCrunch record: RRID:AB_1130818

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (Cell Signaling Technology Cat# 7076, RRID:AB_330924)

Curator: @scibot

SciCrunch record: RRID:AB_330924

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_2848937

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (RRID:CVCL_0063)

Curator: @scibot

SciCrunch record: RRID:CVCL_0063

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_10862910

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: None

Curator: @scibot

SciCrunch record: RRID:AB_2809903

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: (ATCC Cat# HTB-63, RRID:CVCL_0318)

Curator: @scibot

SciCrunch record: RRID:CVCL_0318

What is this?

DOI: 10.1038/s41419-025-07863-y

Resource: RRID:Addgene_97329

Curator: @scibot

SciCrunch record: RRID:Addgene_97329

What is this?