Note: This rebuttal 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
Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.
Reply to the Reviewers:
We would like to thank the reviewers for taking the time to provide us with insightful and constructive comments, which helped us in improving the manuscript. We have performed additional experiments and data analysis while also improving the presentation of our analysis. In addition, we have noted one experiment (bulk RNA-sequencing of pax2a-Low and pax2a-High thyrocyte population) for further revision to address the concerns raised by Reviewer #1.
Please find below our point-by-point response. Please note that the figure references in the response refer to the ‘Revised Manuscript’, ‘Revised Figures’ and ‘Revised Supplementary Data’ file.
Reviewer #1:
Gillotay et al. performed an unbiased profiling of the zebra fish thyroid gland and captured different cell populations by single-cell gene expression analysis. Using bioinformatic tools, they identified seven clusters corresponding to expected, but also poorly characterized, sub-populations, such as non-follicular epithelial cells. They also found two transcriptionally distinct types of thyrocytes and validated this heterogeneity using a new transgenic Pax2a reporter line. Using this tool, they identified and located Pax2a-low and -high thyrocytes within thyroid follicles. Finally, they highlighted a dense intercellular signaling network based on ligands expressed by the diverse sub-populations present in the thyroid and receptors expressed by the thyrocytes. This is a descriptive work calling for more in-depth analyses.
The authors thank the reviewer #1 for acknowledging the strengths of the manuscript and for stating that it could be of interest to a larger audience. We appreciate the reviewer’s advice on supportive experiments and for improving the clarity of the analysis and presentation. We have tried to experimentally address the concerns of the reviewer and hope this provides in-depth substantiation of our observations.
**Major comments on main conclusions :**
Conclusion 1 : Identification of 7 clusters
- This reviewer was particularly surprised by the relative abundance of the different sub-populations identified. For example, 267 thyrocytes out of 6249 cells from the thyroid gland is less than 5% of the total thyroid cell number. In comparison, authors identified three times more immune cells or non-follicular epithelial cells. Authors should comment on these numbers, on the dissection (contaminants?) and dissociation procedure.
Do these relative abundance reflect the proportion of thyrocytes, immune, stromal cells they normally observe in adult thyroid sections ? It would be interesting to have a better resolution for figure 6A in order to evaluate the number of nuclei stained only with DAPI as compared to nuclei stained with DAPI and surrounded by E-cadherin staining. Based on the image, this reviewer seriously doubts that follicular cells represent less than 5% of the total cell number in this organ, considering that the colloid is a cell-free zone.
We fully agree with the reviewer statement. We have clarified the nature of the dissociated tissue in Results, Methods and Fig. 1 C- G. In this, we have performed 3D confocal imaging of the region utilized for single-cell RNA-Seq. (Fig. 1C). Further, we quantified percentage of thyrocytes in transverse sections across the dissociated region (Fig. 1 D – E). Our results demonstrate presence of 5.9 ± 1.9 % thyrocytes in the region. Lastly, we provide the FACS plots of the cells utilized for single-cell RNA-Seq. In this, we obtained around 4% thyrocytes among the live cells. Taken together, our quantifications of the thyrocyte proportion in the tissue matches well with the percentage of thyrocytes obtained in the single-cell atlas, suggesting lack of thyrocyte loss during the procedure.
We agree with the reviewer’s observation that the follicle lumen represents a cell-free zone. However, it is worth noting that the cells surrounding the follicles, in particular gills and stroma, have a higher cell density. Additionally, the zebrafish thyroid follicles sit loosely on the ventral aorta, thus making it difficult to manually separate the follicles from the surrounding tissue without destroying the organ. We avoided injuring the organ to minimize cell-death associated with manual dissection.
- When using the webtool developed on the thyrocyte population, one can notice that only a small fraction of the thyrocyte population expresses common thyroid-specific genes such as Tpo or duox. Was this expected ? It would be interesting to comment on this observation and confirm using standard localization technique to demonstrate that this is real and not due to the sequencing. Is it specific to the zebrafish ? On the other hand, Tg is expressed in most thyrocytes, but surprisingly also in all the clusters at a fairly good level. This should be commented... Is it normal ? due to the sequencing quality ? or clustering ?
We believe these are technical issues related to single-cell sequencing and thank the reviewer for their insight in seeing this.
The non-uniform expression of thyroid-specific marker genes (Tpo and duox) likely represents dropout effects, possibly due to the low expression levels of these genes. To address this issue, we propose to perform bulk RNA-Seq. of thyrocytes (segregated by pax2a expression levels). Bulk RNA-Seq. is more sensitive than single-cell RNA-Seq. and should provide the expression levels of these genes in the thyrocytes.
The expression of Tg in non-thyrocyte population likely represents cross-contamination of free RNAs released from ruptured cells. Since Tg mRNA is highly expressed in the thyroid follicular cells, release of the mRNA from a few injured cells would contaminate the cell suspension, leading to its detection in non-thyroid cells. However, the expression represents background noise signal. To test this, we utilized DecontX, a recently developed approach for background correction (S. Yang et al. 2020). In this, the expression of a gene is modelled as a mixture of expression in the expected population plus background expression. With this, we could robustly reduce Tg mRNA expression in non-thyrocytes (Supp. Fig. 3). This supports our hypothesis that the Tg expression in non-thyrocytes likely represents cross-contamination of mRNA from ruptured cells.
- In order to validate and locate the different populations identified in the thyroid, this reviewer suggests to perform in situ hybridization or immunostaining, based on the specific marker genes identified in each cluster. This experiment could lead to the precise identification of the different sub-populations and their respective localization. These experiments would also help in the interpretation of the cellular interaction network.
We have characterized the different cell types surrounding the thyroid follicles using various reporter lines. The data is presented in Fig. 4.
Conclusion 2 : two distinct types of thyrocytes
- This is an interesting observation. However, from a non-expert it is difficult to understand why the authors propose two populations. Based on the points distribution (Figure 4A), this reviewer would rather identify 3 or 4 clusters but not the two shown in red and blue.... Did the authors impose two populations for the clustering ? Did they perform a permutation test to confirm the pax2a significant fold change seen between clusters is not a false positive generated by the clustering ? Could they show, as a supplementary file, the same graph with points colored based on Pax2a expression ?
We concur with the reviewer that the number of potential clusters in thyrocyte population might be more than two. In-fact, there is no upper limit to the diversity present in the thyrocyte population. However, our message in the manuscript is that the population is not homogenous, and there are at-least two populations based on pax2a expression level. In this regard, we do believe that we have only scratched the tip of the iceberg and further investigation is needed to completely answer this issue. Nonetheless, we are the first to demonstrate genetic heterogeneity among the population.
The clustering of thyrocytes followed the guidelines suggested by Seurat package. For thyrocytes, we utilized the principal components displaying significant deviation from uniform distribution. For cluster identification, we utilized a resolution of 0.3, which was same as the one utilized for clustering the entire organ (further details on this provided as a response to Minor Concern #3 by Reviewer #1). A plot of pax2a expression, along with tg expression, is provided as Supp. Fig. 8.
Finally, to strengthen our observation, we conducted independent analysis of transcriptional diversity in the population using a recently developed method called ROGUE (Ratio of Global Unshifted Entropy) (Liu et al. 2019). The method provides genes that display transcriptional heterogeneity within the cell population. Assessment is based on expression entropy, a measure of the degree of uncertainty, or promiscuity, in the expression of a gene (Teschendorff and Enver 2017). For this, we utilized raw counts of thyrocytes so as to provide an alternative analysis of our data. The analysis, presented as Fig. 6D, demonstrates significant entropy, or transcriptional heterogeneity, for pax2a and cathepsin B (ctsba). The full list of genes displaying transcriptional heterogeneity in thyrocytes is provided as Supp. Table 4.
- This reviewer was also surprised by the relatively "heavy" approach (generation of the Pax2a reporter line) used to demonstrate the existence of two types of thyrocytes. Knowing that the reporter line was validated with a very good Pax2a antibody... The use of the reporter line is a bit short. The authors could for example validate the two populations of Figure 4A using the Pax2a FACS-sorted cells and RT-qPCR.
We completely agree with the suggestion of the author and plan to perform bulk RNA-sequencing using the pax2a reporter line to corroborate our results. This is the advantage provided by the generation of knock-in line. In addition, we have performed antibody staining against endogenous pax2a protein (Fig. 8 E – F), which validates the transcriptional heterogeneity observed in our single-cell RNA-Seq. data.
- In addition, data available in the sequencing dataset could be used to prove that the two populations are really active thyrocytes. This reviewer would suggest to present a table with the expression level of common and thyroid-specific genes such as TshR, Nis, Tpo, Duox, Tg, Pax2, TTF1 and other known transcription factors in the two populations to demonstrate that these two types of cells are indeed thyrocytes. Finally, image quality (Figure 6) could be improved and high-magnification images with several thyrocyte marker could be shown to convince the readers.
We strongly agree with the reviewer that this is a very important concern to address. To address this, we have taken three steps:
We have included the expression level of tg in the two populations in Fig. 6C and Supp. Fig. 8.
We performed antibody staining against pax2a on thin section obtained from Tg(tg:nls-EGFP) animals (Fig. 8 E-F). In this, we observed pax2a-Low cells with tg reporter expression, suggesting that they are indeed differentiated thyrocytes.
We plan to perform bulk RNA-sequencing of cells from pax2a-Low and pax2a-High population. This will allow us to validate the transcriptional differences observed by single-cell RNA-Seq., and allow us to demonstrate expression of thyroid-specific markers genes that are missing from our dataset (for instance, duox).
Conclusion 3 : cellular interaction network
- Most of the interactions revealed by the analysis seem to belong to the extracellular matrix and not to classical ligands such a Wnt, TGFb, FGF, PDGF,..... could the authors comment on this ? Considering that both endothelial cells and epithelial cells assemble their own basement membrane, the analysis will obviously reveal interactions between endothelial cells and epithelial cells....
We appreciate that the reviewer pointed this out. The enrichment of ligands related to extracellular matrix, and not growth factors, likely represents the homeostatic nature of the organ. Growth at 2 and 8 mpf is low (if not absent). Correspondingly, gene expression related to development and cell-cycle might be reduced. As stated in response to the next concern, extending the atlas to juvenile stage (1 mpf) would be beneficial to understand the regulators of cell-cycle.
However, to improve the cellular interaction network, we have incorporated physical information from the characterization of cell-populations surrounding the thyroid follicles (Fig. 4). Our experiments suggested that stromal, gills and NFE do not physically contact the thyrocytes. Thus, interactions based on ligands incorporated into the cell membrane were removed for these cell-populations.
**Minor comments to improve the Ms :**
- Could you explain how from 2 x 12 000 FACS-sorted live-cells (from six animals at each stage; 2 mpf and 8 mpf) you obtain 6249 cells (pooled of 2 mpf and 8 mpf), and why the two stages were first sorted separately and then pooled (?), as no differential analysis is carried out for the two stages.
The number of cells obtained for analysis represents cells that were successfully incorporated into droplets during library preparation and generated high-quality data that passed quality control (Supp. Fig. 2). FACS sorted cells were utilized for droplet generation using 10X Chromium that encapsulates cells with single-Poisson distribution (Zheng et al. 2017). This leads to approximately 50% cell capture rate, which is the ratio of the number of cells detected by sequencing and the number of cells loaded. Thus, we obtained 13,106 sequenced cells from 24,000 input cells (54.6 % cell capture rate). Further, the quality control criteria removed 6,857 low-quality cells (52.3 % dropout rate). We chose a stringent cut-off for quality control so as to remove technical artefacts from the analysis. We have added these detail to the Result section.
We pooled the two samples as the stages represent the range of homeostasis in zebrafish. We decided not to include differential expression between the two stages as the number of cells in multiple clusters were too low for individual stages (less than 100), and thus not trustworthy. It future, it would be of interest to extend the analysis for rapidly-growing juvenile (less than 1 mpf) and old-age (greater than 1.5 ypf animals) animals and to perform single-cell or bulk RNA-Seq. with high cell numbers. We have mentioned this drawback in the discussion section.
- Which method was used for the graph-based clustering ? KNN ? Louvain ? Random walk ?
The details have been added to the Method section. Specifically, the clustering was performed using graph-based method, shared nearest neighbour (SNN), which is default for Seurat 2.3 package.
- How did you define the numbers of clusters ?
The number of clusters were defined by using the first five principal components as they displayed significant deviation from uniform distribution as accessed by JackStraw analysis. Further, a resolution of 0.3 was used in Seurat as the clusters generated by this parameter could be annotated using a cell-type specific marker from literature.
- Figure 4B, the color-code for the expression level would help the reader.
The color code has been added (Fig. 6B in revision).
- Figure 4C, violin plot for Pax2a: why do we find cells that do not express this gene in the two populations ? The same is true for tbx2a and ahnak ... is the clustering optimal ?
The detected expression of pax2a depends on its biological expression and technical dropout rate. Thus, the pax2a-High cluster also contain cells with no detectable expression of pax2a. Similar detection dynamics can be expected for other genes.
We have experimentally validated the variability in pax2a expression using antibody staining for endogenous pax2a protein in tg:nls-EGFP transgenic line (Fig. 8 E-F). With this, we can validate the presence of pax2a heterogeneity within thyrocytes.
- Figure 4C, blue violin plot for ptp4a3 does not seem to fit with the distribution of the points.
Due to the high number of cells that do not express ptp4a3, the cells collapse on each other at the bottom of the graph, thus making the violin plot seem different from the distribution. However, the plots were made with Seurat without changing any parameters.
- what is the function of tbx2a, ahnak, ptp4a3 and dusp5, which are not mentioned in the text.
The genes have been implicated in regulation of cell proliferation. However, we acknowledge that the evidence based on a handful of genes needs to be strengthened. For this, we have removed the figure panels from the revision, and will instead identify genetic markers based on bulk RNA-sequencing analysis of pax2a-High and pax2a-Low population.
- Line 195: "pax2amKO2 reporter expression perfectly overlapped with PAX2A antibody staining". This reviewer would be more cautious as the images (Fig. 5 C, D and F) do not show a perfect colocalization: one can observe only blue or only red staining.
We have edited the text to “The pax2apax2a-T2A-mKO2 (abbreviated as pax2amKO2) reporter expression overlapped with PAX2A antibody staining in a majority of regions at 9.5 hours post-fertilization.” The regions with single colors in Fig. 5C (Fig. 7C in revision) are due to differences in staining intensity between different regions.
- Line 246: it is proposed to "study the functional and replicative differences among the two sub-populations of thyrocytes". This reviewer completely agrees and the suggestion made (vide supra) to use the datasets to assemble a table with the expression level of common and thyroid-specific genes such as TshR, Nis, Tpo, Duox, Tg, Pax2, TTF1 and other known transcription factors in the two populations could already give some indications on the functionality of these two types of cells. Expression of genes involved in cell-cycle control and/or apoptosis would be another possibility to better characterize the two populations. Lastly, the authors could perform the comparative analysis of ligand-receptor pairs between these two sub-populations to examine if they differentially interact with their environment.
We agree with the reviewer for the need to characterize the two populations in detail. Using the current dataset, we are limited to the 265 genes differentially expressed between the two thyrocyte states. Thus, we propose to perform bulk RNA-sequencing of the two populations to obtain a better picture of the cellular identity. In this, we will perform sequencing of each population in triplicates. With this, we will avoid the dropout effects suffered by the single-cell analysis. The experiment would demonstrate the differentiation status of each cell population, and provide insights into other pathways active within each population. Further, we will identify ligands and receptors that are enriched in a particular population.
Text improvement:
Intro: thyroglobulin (TG) appears twice (line 46 and 49)
Results: Fig. 5 (not 8) (line 203 and 205)
Figure 3: stromal? (not skeletal)
Figure 4: fold change scale is missing
Figure 5: Thyroid gland (Gland)
Figure sup 2: Fluorescence-activated cell sorting (FACS) of zebrafish thyroid gland
Figure sup 3: number of unique molecular identifier
Figure sup 4: "are present in the zebrafish"
We thank the reviewer for pointing these errors. We have edited them.
Reviewer #1 (Significance (Required)):
The work performed by Gillotay et al. is clearly novel but descriptive. It provides a useful database to propose hypotheses to further study the thyroid gland. The single-cell RNA-seq analysis brings a molecular appreciation of the various thyroid cell populations, thyrocyte heterogeneity and intercellular signaling network. Although focused on the thyroid, results will be of interest to a larger audience than the thyroid community, especially the demonstration (if further and better validated) of thyrocyte heterogeneity and the intercellular communication possibilities.
In response to comments by Reviewer 1, we plan to perform bulk RNA-Sequencing of pax2a-Low and pax2a-High thyrocytes. We believe that this will help address all the remaining concerns of the reviewer.
Reviewer #2:__
__
**Summary**
In this manuscript the authors present a single-cell transcriptome atlas of the zebrafish thyroid gland (possibly also including some adjacent tissues depending on how the dissection was performed, see comments below). The atlas includes cell clusters that are expected to be found in the thyroid of any higher vertebrate species (thyrocytes, blood vessels, lymphatic vessels, immune cells and fibroblasts), but also musculature/gills and a less well-defined population of non-follicular epithelium. The data will be made publicly available as a resource, by what seems to be a user-friendly web-interface (more accessible to a broader audience than customary raw sequencing data deposition, that I suppose will also be provided). The results are used to describe putative autocrine or paracrine signaling networks. The authors are able to further subdivide the thyrocyte cell cluster into two sub-populations with different transcriptomic features. Interestingly, these populations differ in their expression of for instance the key transcription factor pax2a, which is further demonstrated by the use of a novel zebrafish reporter strain.
In general, this is a clearly interesting, novel, nicely structured and well written manuscript and the data presented seems to be of high quality.
We would like to thank reviewer 2 for the constructive comments. We are glad that the reviewer finds our work interesting, novel and of high quality. We appreciate the reviewer’s advice on additional experiments, analysis and on improving the clarity of the text. We have addressed all the concerns raised, and hope that our revised manuscript satisfies the reviewer.
**Major comments**
Key conclusions are convincing and performed with scientific rigor. As will be further discussed below the seemingly superficial description of the extent of tissue that was subjected to transcriptome analysis makes it a bit difficult to assess reproducibility outside the authors' lab.
We acknowledge the lack of clarity in the description of the tissue utilized for single-cell analysis. We have corrected this providing a detailed step-by-step protocol for dissecting the organ in Methods Section, titled ‘Dissection of the zebrafish thyroid gland’. Additionally, using immunofluorescence based imaging of the region and FACS, we estimate the proportion of thyroid follicular cells within the region. The results are presented in Fig. 1 C – G.
As far as I can see very little methodological detail or information is provided about how the dissection of the thyroid region was performed, more than that "the thyroid gland was collected" or that "we dissected out the entire thyroid gland". This is essential to understand the significance of the cell populations that are described based on the transcriptomic data. The section "Tissue collection" describes dissection of the thyroid for whole-mount imaging. From Fig. 6A it seems that substantial amounts of non-thyroid tissue are included in this dissection. Is it the same kind of dissection that was performed for transcriptomics? Was the string of thyroid follicles shelled out from their surroundings or was some kind of en bloc dissection, including other neighboring tissues, performed (as suggested from the transcriptome cell populations data including e.g. gill transcripts)? In the latter case it would be good if the authors discuss in more detail what other tissues or structures that are expected to also be included in the dissected tissue and transcriptomic data.
In response to this concern of the reviewer, we have improved the clarity of the text in Results and Methods section. We have added the following text to the Result section, “The zebrafish thyroid gland is composed of follicles scattered in the soft tissue surrounding the ventral aorta (Fig. 1 A, B). Ventral aorta extends from the outflow tract of the zebrafish heart and carries blood from the ventricle to the gills. Dissection of the ventral aorta associated region (detailed in Methods section) provided us with tissue that included the thyroid follicles and parts of zebrafish gills (Fig. 1C). Using Tg(tg:nls-EGFP) transgenic line, which labels thyrocytes with nuclear green fluorescence (Fig. 1D), we estimated presence of 5.9 ± 1.9 % thyrocytes within the dissociated region (Fig. 1E).”
Further, the Methods section now defines a step-by-step protocol for dissociation (‘Dissection of the zebrafish thyroid gland’).
In addition, we have improved the characterization of the dissected region, as stated in response to the next comment.
It would facilitate understanding if the thyroid is outlined in Fig. 1A as well as the region that was dissected for downstream single cell sequencing.
A whole-mount 3D reconstruction of the region is presented in Fig. 1C. A transverse section from the region is presented as Fig. 1D, while quantification of the percentage of thyrocytes in the transverse sections is provided in Fig. 1E.
The clusters seem logical given what cell types that could be expected in the region (but depending on how dissection was performed). The only exception is cluster 7 (non-follicular epithelium; NFE). I do understand that relative sizes of the clusters do not necessarily reflect the endogenous relative abundance of different cell types, as I guess they may be more or less prone to enzymatic dissociation, survival etc. Nevertheless, the number of cells in the NFE cluster (831 cells) seems sizeable relative to the number of thyrocytes (267 cells). In my opinion, a major weakness of the current manuscript is that little detail is provided about this cell population and that no attempt to at least spatially localize these cells is presented.
Detailed characterization of the cell-populations surrounding the thyroid follicles is now presented in Fig. 4. In addition, we have quantified the percentage of thyrocytes in the region (Fig. 1 E), which demonstrates that thyrocytes represent 5.9 ± 1.9 % of the cell-population. Additionally, we have presented FACS analysis of the dissociated region as Fig. 1 F - G, which corroborates the imaging based quantification. Both quantifications are in the same range as the proportion of thyrocytes identified in the single-cell analysis (4.27 % - 267 out of 6249 cells). Thus, we do not believe that cell-loss had a big impact on the relative abundance of thyrocytes in the single-cell atlas.
A detailed characterization of NFE cells is provided in response to the next three comments, which includes visualization of the population using TP63 / p63 antibody staining in Fig. 4D. Particularly, Fig. 4D contains 72 thyrocytes and 302 TP63+ nuclei, thereby pointing to the higher relative abundance of NFE in the region.
The NFE cells are characterized by TP63 expression and the authors speculate that they might show homology to main cells of solid cell nests. From previous zebrafish literature it seems like what is supposedly ultimobranchial bodies (or ultimobranchial glands more similar to those in avian species) are located rather distant from the thyroid follicles (Alt et al 2006). Is it possible that these structures are included in the dissection that has been performed? As solid cell nests are supposed to be ultimobranchial body remnants (with calcitonin positive and calcitonin negative epithelial cells) it would be good if the authors discuss in more detail what is known about the ultimobranchial bodies in zebrafish, if they are located inside the zebrafish thyroid, in an anatomical region that is included in the dissected tissue of this study or in a region that is likely not included.
As stated by the reviewer, the ultimobranchial bodies lie distant to the thyroid gland. They lie as a pair of follicles on top of the sinus venous (Alt et al. 2006), which is a blood vessel that delivers blood to the atrium. In contrast the thyroid follicles sit loosely on top of ventral aorta that connects to the ventricle via the outflow tract (Fig. 1B). Thus, the collection of the ultimobranchial bodies is not expected. This is also corroborated by the absence of calcitonin (calca) expression in the NFE (Supp. Table 1). This has been added to Discussion section.
In higher vertebrates, P63 expression is typically seen in basal cells of stratified epithelia (as for instance in the esophagus), in myoepithelial cells, in the urothelium and in the thymus. Is it possible that the TP63 expressing NFE population corresponds to cells of the zebrafish thymus (that might perhaps explain the seemingly large immune cell population in cluster 4)? Could TP63 expressing NFE cells represent the esophagus (if included in the dissection)? As so little detail is provided about the dissection procedure this opens up for speculation and it would be good if the authors discuss these possibilities, as some of them might perhaps be unlikely or even impossible given regional anatomy of the zebrafish and how the dissection was performed.
The dissection region is now characterised in detail in Fig. 1 C – E. The presence of immune cells (macrophages) in the proximity of thyroid follicles is validated in Fig. 4B. The presence of NFE is presented as Fig. 4D, and explained in detail in response to the next comment.
To gain better understanding of the sizeable TP63 expressing NFE population the authors briefly mention the possibility of future studies utilizing a TP63 reporter. If a reporter line is not available, immunofluorescent detection of P63 (as presented for PAX2A in Fig. 5 and E-cadherin in Fig. 6) would be desirable to provide more insight into the location of the NFE population. Given the proficiency the authors demonstrate in this manuscript when it comes to zebrafish imaging, at least whole-mount immunostaining of P63 in the region that was dissected for transcriptome analysis seems clearly feasible, both with respect to resources and time needed (perhaps in the range of 1-3 months).
To clarify the presence of NFE cells, we have followed the suggestion of the reviewer and performed immunostaining against TP63. The result is presented as Fig. 4D. The staining was performed on thin (8 µm) sections, allowing us to ensure uniform antibody penetration. As depicted in the image, TP63+ cells are part of the gills. This population possibly represents a progenitor population, similar to the TP63+ basal layer in the zebrafish (Guzman et al. 2013) and mammalian (A. Yang et al. 1999) epithelium. Additionally, a sub-set of TP63+ cells were observed in the region between the thyroid follicles and gills. Our data provides an exciting opportunity for an in-depth study of these cells in future, particularly using tp63-regulatory region driven transgenic reporter and Cre lines.
**Minor comments**
It is a little bit confusing that different color coding for the various cell populations is used in Fig. 3B as compared to Figs. 1 and 2.
The color coding for Fig. 3B (Fig. 5B in revised manuscript) has been modified in accordance to the once used in Fig. 1 and 2.
In the discussion of the intercellular interaction network (Fig. 3B) the authors clearly point out that anatomical barriers are not modelled. Nevertheless, it would be more informative if this description was able to sort out ligands that are secreted, from ligands that are not secreted and would require physical interaction between thyrocytes and a different cell population for signaling to occur.
We have now improved the intercellular network to resemble the putative physical interactions. By characterizing the different cell-populations present in the atlas (Fig. 4), we recognized that gills, stromal and NFE are not in physical proximity of the thyrocytes. Thus, these three cell-populations would not be able to communicate via ligands attached to the cell membrane. Hence, we have pruned the network to remove cell-membrane attached ligands for these three cell-populations. Only secreted ligands were considered for the mentioned cell-populations. In accordance, the figure and Supplementary Table 2 has been updated.
The authors describe thyroid solid cell nests as "... lumen containing irregular structures located between the thyroid lobes in mammals...". Solid cell nests of the thyroid in higher vertebrates (e.g. humans and dogs) are located within the thyroid lobes and not between the lobes (i.e. the right and left thyroid lobes). Moreover, the authors write that "Recently, epithelial cells have been reported in a structure called the Solid Cell Nests (SCN) of the thyroid..." and give reference to a paper from 1988. If that is recent or not might be a matter of opinion, but to the best of my knowledge, solid cell nests were describe already by Getzowa in 1907 and I suppose that the epithelial identity (or at least epithelioid morphology) has been appreciated for long.
We thank the reviewer for pointing this out. We have added the reference to the original study by Dr. Sophia Getzowa identifying SCN (Getzowa 1907). As the original study is in German, we have also added a reference to a recent article attributing the discovery of SCN to Dr. Getzowa and performing immunohistochemistry analysis of the tissue (Ríos Moreno et al. 2011). Notably, the authors note the presence of TP63 staining, along with the absence of TG and Calcitonin staining, in SCN main cells – similar to the expression profile of NFE in our atlas. Finally, we have edited the text to accurately describe their location in the mammalian thyroid gland.
Reviewer #2 (Significance (Required)):
The authors provide a single-cell transcriptomic atlas of the zebrafish thyroid gland. To the best of my knowledge this is certainly a unique resource. In that sense it provides novel and significant information that will surely facilitate our further understanding of thyroid biology. It will surely be of great interest and value to the thyroid community, but probably also to a wider audience interested in e.g. zebrafish biology, endodermal biology and the biology of endocrine glands. Their finding and direct demonstration of transcriptional heterogeneity within the thyrocyte population is very interesting, also in different contexts of thyroid disease. However, I leave it to other referees to comment on the conceptual uniqueness of the current manuscript (i.e. a single-cell transcriptomic atlas of a zebrafish organ). Does it provide conceptually unique information, or does it add to an expanding collection of single-cell transcriptomic atlases of zebrafish organs?
References:
Alt, Burkhard, Saskia Reibe, Natalia M. Feitosa, Osama A. Elsalini, Thomas Wendl, and Klaus B. Rohr. 2006. “Analysis of Origin and Growth of the Thyroid Gland in Zebrafish.” Developmental Dynamics 235 (7): 1872–83. https://doi.org/10.1002/dvdy.20831.
Getzowa, Sophia. 1907. “Über Die Glandula Parathyreodeaa, Intrathyreoideale Zellhaufen Derselben Und Reste Des Postbranchialen Körpers.” Virchows Archiv Für Pathologische Anatomie Und Physiologie Und Für Klinische Medizin 188 (2): 181–235. https://doi.org/10.1007/BF01945893.
Guzman, A., J. L. Ramos-Balderas, S. Carrillo-Rosas, and E. Maldonado. 2013. “A Stem Cell Proliferation Burst Forms New Layers of P63 Expressing Suprabasal Cells during Zebrafish Postembryonic Epidermal Development.” Biology Open 2 (11): 1179–86. https://doi.org/10.1242/bio.20136023.
Liu, Baolin, Chenwei Li, Ziyi Li, Xianwen Ren, and Zemin Zhang. 2019. “ROGUE: An Entropy-Based Universal Metric for Assessing the Purity of Single Cell Population.” BioRxiv, January, 819581. https://doi.org/10.1101/819581.
Ríos Moreno, María José, Hugo Galera-Ruiz, Manuel De Miguel, María Inés Carmona López, Matilde Illanes, and Hugo Galera-Davidson. 2011. “Inmunohistochemical Profile of Solid Cell Nest of Thyroid Gland.” Endocrine Pathology 22 (1): 35–39. https://doi.org/10.1007/s12022-010-9145-4.
Teschendorff, Andrew E., and Tariq Enver. 2017. “Single-Cell Entropy for Accurate Estimation of Differentiation Potency from a Cell’s Transcriptome.” Nature Communications 8 (1): 15599. https://doi.org/10.1038/ncomms15599.
Yang, A, R Schweitzer, D Sun, M Kaghad, N Walker, R T Bronson, C Tabin, et al. 1999. “P63 Is Essential for Regenerative Proliferation in Limb, Craniofacial and Epithelial Development.” Nature 398 (6729): 714–18. https://doi.org/10.1038/19539.
Yang, Shiyi, Sean E. Corbett, Yusuke Koga, Zhe Wang, W Evan Johnson, Masanao Yajima, and Joshua D. Campbell. 2020. “Decontamination of Ambient RNA in Single-Cell RNA-Seq with DecontX.” Genome Biology 21 (1): 57. https://doi.org/10.1186/s13059-020-1950-6.
Zheng, Grace X. Y., Jessica M. Terry, Phillip Belgrader, Paul Ryvkin, Zachary W. Bent, Ryan Wilson, Solongo B. Ziraldo, et al. 2017. “Massively Parallel Digital Transcriptional Profiling of Single Cells.” Nature Communications 8 (1): 14049. https://doi.org/10.1038/ncomms14049.
