Idiosyncratic and generic single nuclei and spatial transcriptional patterns in papillary and anaplastic thyroid cancers

  1. Adrien Tourneur
  2. Joel Rodrigues Vitória
  3. Manuel Saiselet
  4. Ligia Craciun
  5. Denis Larsimont
  6. Anne Lefort
  7. Frederick Libert
  8. Carine Maenhaut
  9. Sabine Costagliola
  10. Maxime Tarabichi
  11. Mirian Romitti
  12. Vincent Detours

Reviewed by

Read the full article See related articles

Discuss this preprint

Start a discussion What are Sciety discussions?

Listed in

Log in to save this article

Abstract

Sixty percent of papillary thyroid cancers (PTCs) are driven by BRAF v600E , a mutation associated with high inter- and intra-tumoral heterogeneity. PTCs may become highly aggressive anaplastic thyroid cancers (ATC). While single cell transcriptomics may resolve this heterogeneity, it is potentially confounded by technical effects whose correction may dampen inter-tumor variations. Here we profiled ATCs and BRAF v600E PTCs with single nuclei RNA-seq and spatial transcriptomics, and an experimental design disentangling biological and technical variations. It reveals that much transcriptional variation in cancer cells and several immune cell types is idiosyncratic, i.e. tumor-specific, a phenomenon obscured by batch integration in a number of single cell studies. Idiosyncrasies are associated in some cases with genomic aberrations and global tissue states like hypoxia. Beyond idiosyncrasies, differentiation markers SLC5A5 (NIS), TPO, TG and TSHR are lost in a sequence mirrored by their gain during human thyroid organoids maturation, suggesting a new classification of cancer cell states. PTC cells retain TSHR expression and show features of partial EMT with a massive expression of FN1 , which promotes proliferation via an autocrine loop. In contrast, ATCs undergo full blown EMT, with expression of mesenchymal extracellular components and loss of TSHR . Finally, we show that the microenvironment of cancer cells is driven by inflammation. These findings may help future stratifications of BRAF v600E PTCs.

Article activity feed

  1. Review Commons

    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

    This is a reply/revision plan, not definitive. Planned and already implemented revisions are underlined.

    First of all, we wish to express our gratitude to the reviewers: they helped to improve the paper.

    Reviewer #1:* **

    Reviewer #1 wrote: Major Comments: 1.Differential gene/pathway analysis across epithelial clusters: What are the differential genes or pathways among the epithelial clusters? Without CCA/Harmony integration, do the tumor subgroups show distinct differences? In addition, I suggest applying NMF or hdWGCNA to identify shared modules and test whether ATC and PTC harbor overlapping regulatory modules.

    Reply plan: Both reviewers suggested some regulatory network analysis. We proposed to run SCENIC+ (Nature Methods, 2023, https://doi.org/10.1038/s41592-023-01938-4) on our data__.__

    __Reviewer #1 wrote: __2.Validation of TSHR/TPO-based subgrouping: While the TSHR/TPO grouping appears appropriate for stratification at the single-cell level, it is necessary to exclude sequencing depth as a confounding factor. Should validate the existence of these subpopulations using mIHC/IF on corresponding samples. *

    __Reply plan: __We made claims about RNA expression, not protein expression. Thus, validation should be at the RNA level:

    • We already replicated part of our analysis on the dataset published by Lu et al. (JCI 2023, https://doi.org/10.1172/JCI169653), see Figs. 3 and 4. This effort will be extended to all single cell analysis results from our study in the revised paper.
    • We will also present plots demonstrating that the sequencing depth is similar in the different cancer cell subgroups-further excluding it as a confounding factor. __Reviewer #1 wrote: __*3.Impact of mutational differences on conclusions: According to Supplementary Table 1, almost all PTC cases carried BRAF mutations, whereas four ATC patients harbored no BRAF mutation. Could this difference influence the conclusions of the study? Although the authors briefly mention this in the Discussion, a more thorough clarification is warranted. *

    Reply plan: The dataset of Lu et al. includes BRAF-mutated ATCs along with BRAF-mutated PTCs. Therefore, the replication mentioned earlier will also address those concerns. In fact, Fig. 4E-I already confirm in Lu et al. data the ordered loss of markers. Replication will be extended to other results of the study and be more emphasized in the paper.

    __Reviewer #1 wrote: __4. The statement "Myeloid and T cells also grouped in specific clusters" seems descriptive. Is this clustering biologically meaningful? Please elaborate.

    __Reply plan: __This is an important point, and accordingly, a cell mixing experiment was specifically designed to sort apart technical effects from biological effects. We therefore know with certainty that the myeloid and T cell patients-specific clusters are the result of biological variation (Fig. 1). We further demonstrate that part of this variation is associated with hypoxia (Supp. Fig 4). So yes, the clustering is biologically meaningful.

    __Reviewer #1 wrote: __Minor Comments: In Figure 2C, the "Epith TSHR-" population resembles myeloid cells. Could the authors clarify why this is the case? For the correlation analysis in Figure 2C, were highly variable genes or all genes used?

    Reply plan: There is a simple explanation: The Epith TSHR- population the reviewer is referring to are cells from anaplastic thyroid cancers (ATC), which are tumors notoriously infiltrated by macrophages (Supp. Fig. 4). A high correlation of Epith TSHR- and macrophages proportion across our panel of ATC and papillary cancer (PTC) is therefore expected. Among other things, Fig. 2C shows that high correlation, but it is not meant to and does not show that Epith TSHR- and macrophages "resemble" one another. It shows that their proportions are highly positively correlated. This correlation analysis does not rely on gene expression but on cell type proportions. It measures co-occurrence rather than resemblance. The text has been clarified in order to prevent any confusion.

    __Reviewer #2: __

    __Reviewer #2 wrote: __1. This study largely confirms established facts that 1) PTC due to BRAF driver mutation is a heterogeneous tumour entity and 2) ATC is the most dedifferentiated of all thyroid cancers. Although interesting, observations of a highly variable tissue cell composition including immune cells and the gradual loss of thyroid differentiation markers, in part linked to tumor subclone development featured by altered chromosomal copy numbers, are thus not surprising.

    __Reply plan: __We wish to respectfully express our take on this perception of the work:

    • There is a difference between conjecturing a high heterogeneity in the cell composition of thyroid cancers and establishing it with the level of accuracy and quantitative rigor our analysis provides. The extreme amplitude of that variation was surprising to us: the size of the microenvironment makes from 8.4 to 80% of the cells in PTCs driven by the same BRAF mutation.
    • We don't simply show that a subclone characterized by a large number of copy number events is less differentiated. We go all the way proving that those copy number alterations are associated with specific cell states that produce specific histology (Fig. 5). It required a combination of single cell transcriptomics, spatial transcriptomics and sophisticated computational analysis to establish that connection between genomic changes and histology. The fragmentation of epithelial sheets uncovered from CNV analysis had escaped the attention of pathologist colleagues and ours at first, this is not a parameter typically assessed in diagnostic, to our knowledge.
    • We don't simply show that there is a gradual loss of differentiation markers: this loss is ordered in a very specific way that mirrors the gain of markers during thyroid organoid differentiation.

    __Reviewer #2 wrote: __2. Considering tumor progression, comparison of PTC and ATC should preferably include specimens with the same driver mutation (BRAF or RAS), which is not the case here. This notion should be more clearly explained to readers. An optional improvement would be to conduct similar analyses on an ATC specimen that contains more differentiated PTC tumor portions arguably suggesting that PTC progresses to ATC (by mechanisms that are still largely unexplored).

    __Reply plan: __This is clearly a limitation of our study. As already proposed in our reply to reviewer number one, we will extend to all our single cell results the replication of our analysis in the dataset of Lu e al., which includes ATCs and PTCs harboring the BRAF-mutation.

    __Reviewer #2 wrote: __3. Comments on findings of lymphocytic infiltration need to be balanced. Although autoimmune thyroid disease in infered a risk factor of developing malignancy it is unlikely that the majority of TCGA samples of PTC is associated with thyroiditis as indicated in Fig. 3 and Suppl Fig. 3. Immune cell abundance may rather reflect the tumor immune microenvironment (TIME).

    __Reply plan: __The figure the reviewer is referring to demonstrates that PTC occurring in a background of thyroiditis also has a higher proportion of B cells. We did not claim, and the figure did not show, that "the majority of TCGA samples of PTC is associated with thyroiditis", because they don't. This point has been clarified.

    __Reviewer #2 wrote: __4. Some tissue sections seem of quite poor quality either shape-wise of containing rifts e.g. PTC7 in Fig. 3 and PTC2 in Fig. 5. The authors should explain whether and how this might influence analysis.

    __Reply plan: __Spatial transcriptomics is typically performed on frozen sections. Frozen sections, which are obviously of lower visual quality than slice from FFPE preserved samples. Since no computational analyses were performed on the image, this lower quality has no impact on our results. Regarding RNA quality, the RINs were >7 for all tumors. RINs are now presented in Supp table S1.

    __Reviewer #2 wrote: __The experiment on mouse ESC/organoids (Fig. 6H-J) does not show much of an expected enhanced thyroid progenitor cell proliferation after induction of the mutant Braf allele by tamoxifen, which raises doupt whether the subsequent promoted growth by fibronectin at all is oncogene-related. This differs from the impact of BrafCA activation along with mouse thyroid development in vivo (Schoultz et al iScience 2023 PMID: 37534159). In the same experimental setup, it appears that mutant Braf prevents follicle formation (Fig. 6I). A control experiment investigating the influence of fibronectin in the absence of oncogene activation should be conclusive. The effect of Braf and fibronectin on thyroid organoid structure and function should be better explained, if necessary based on complementary experiments, and discussed in relation to the claimed association of fibronectin expression to "...low amounts of thyroid differentiation markers...) and "...loss of epithelial structure (PTC7, Figure 6E)." in the previous section of Results.

    __Reply plan: __The induction of the mutant Braf allele for 7 days increases the percentage of BrdU+ cells by 1.43 fold (p-value for Wilcoxon test = 0.035). The effects observed by Schoultz et al. are certainly more dramatic, but they result from an oncogenic activity spanning 1 to 6 months (4 to 26 times longer) in an in vivo model. Most importantly, oncogenic activity is initiated in Nkx2.1+ cells and not Tg+ cells, thus much earlier during development. These two models are thus not comparable. As for the effects of fibronectin on thyroid structure, we do not claim that our organoid model recapitulates the complex interactions between cancer cells and their microenvironment that shapes tissue morphology in vivo. This is now clarified in the text.

    We presented controls with no oncogene expression and no Fn1, controls with oncogene induction and no FN1 and organoids with oncogene induction and Fn1 treatment. This alone establishes the effect of Fn1 on induced organoids, which was our goal. We regard it as a novel and interesting but non-essential development in our paper.

    As the reviewer points out, while our results show an increased proliferation in Braf-mutated organoids treated with Fn1, they do not allow us to conclude on any potential interaction between Fn1 and the oncogenic process. The suggested experiment with Fn1 in absence of oncogene activation would add information, but we cannot follow up for practical lab management reasons detailed in Section 4 below.

    __Reviewer #2 wrote: __6. Concerning EMT profiling (Supplementary Fig. 7B) , there is a great similarity of ATC tumor cells and fibroblasts, and as stated in the text the malignant status of the former is indicated by chromosomal aberrations (refering to Suppl fig. 6). However, looking at Suppl. Fig. 7B it is evident that fractions of cells identified as fibroblasts express TG and TSHR suggesting mismatch. How was this comparison done in order to exclude mismatch? Is there no other profiled markers that distinguish cancer cells from stromal cells that can support conclusions?*

    Reply plan: TG-a thyrocyte marker-seems expressed by fibroblasts in Supplementary Figure 7B. The reviewer suggests this could be caused by an incomplete distinction between bona fide fibroblasts and thyrocytes in advanced EMT state. We argue that

    • Ambient TG RNA leaking out of thyrocytes nuclei contaminates the transcriptomes of all cell types. It is a well-known technical problem, with dedicated software packages to mitigate it. We preprocessed our data with one of them, SoupX, which corrected for most, but not all, ambient RNA contamination.
    • The plot below shows that there is nothing special about fibroblasts in that respect. For example, B and T cells are contaminated by TG at levels comparable to fibroblasts, endothelial cells and pericyte to higher levels.
    • In addition, the UMAP of Fig. 2A shows that EMT cells and fibroblast form very distinct clusters. Furthermore, the fibroblast cluster but not the two EMT clusters contain cells from PTC, and the PTC cluster do not contain cells with DNA copy number aberration. Thus, although both EMT cells and fibroblasts express the typical mesenchymal marker of Supplementary Fig. 7B, they are easy to distinguish on the basis of their overall transcriptomes.
    • The panel below has been added to the Supplementary Figure 7B. [Panel cannot be displayed here]

    __Reviewer #2 wrote: __*In the same figure, it appears there are no clear differences in EMT marker expression among PTC samples regardless of differentiation state, suggesting that the gradual loss of thyroid differentiation in PTC tumor cells and EMT are not parallel and potentially linked phenomena? Please clarify this dissociation of results. Is possible that refocusing on other EMT markers than the top 10, of which almost all concerns various collagen genes, might better reveal partial EMT in PTCs?

    __Reply plan: __The technical basis of this comment is related to the previous point. Our perception is that the mesenchymal markers in Supplementary Fig. 7B show a binary effect, i.e. strong expression in ATC and no expression in PTC (beyond ambient RNA noise)-not a gradual effect. Thus, there is no correlation of COL1A1 and other mesenchymal markers with dedifferentiation in PTC as these markers are not expressed beyond the noise level of the experiments. A lot has been written about EMT in PTC, but one of the findings of our study is that while ATC undergo full EMT, EMT in PTC is very limited. PTC express FN1 but no other major mesenchymal markers such as collagens I and III, for example.

    __Reviewer #2 wrote: __*7. According to Suppl. Table 1, the ATC2 tumor does not harbor any mutations. What about chromosomal aberrations, was that included in analysis? Considering previous consistent reports of a high mutation burden in ATC, if not supported by other data (clinical, pathological) the diagnosis might be questioned for this particular case included in multiple analyses of the present study.

    Reply: There is little doubt about the diagnostic of ATC2 by our pathologist collaborators

    • The histology of this tumor is strikingly anaplastic, i.e. without structure, as shown in the image below.
    • This tumor has a high level of macrophages infiltration typical of ATCs (Supplementary Fig. 4).

    __Reviewer #2 wrote: __Minor comments: -The logical order of presentation of Results might benefit from first presenting specific PTC data following by ATC dito. I´m thinking of swapping the section of EMT in ATC to end of Results.*

    *Reply plan: We miss why the reviewer thinks that way. We believe that discussing the microenvironment, then tumor cells bring conciseness and clarity about how we propose to stratify the latter. By contrast, the suggested tumor type-centered structure entails going back and forth between the microenvironment and tumor cells, diluting the messages about both.

    __Reviewer #2 wrote: __-Methods paragraph "Mouse ESC-derived thyroid organoids experiments" (starting with "ccc") seems to be missing some essential information.

    *Reply plan: A sentence was missing, indeed, and has been re-introduced in the manuscript. We thank the reviewer for catching that error.

  2. Review Commons

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

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    Summary:

    Expression pattern profiling of human thyroid cancer tissues by combining single cell/nuclei RNAseq analysis, spatial transcriptomics and immunofluorescence on corresponding tumor histologic sections. Papillary and anaplastic thyroid carcinomas (PTC n=10 and ATC n=4) were compared; some data were extracted from TCGA. The results indicate that ATCs consists of completely dedifferentiated tumor cells whereas PTCs show variable levels of dedifferentiation, which in a sense mimics the the reverse process of thyroid differentiation as observed in stem cell-based organoids. Moreover, PTC and ATC tumors show different levels of epithelial-mesenchymal transition. Fibronectin is inferred a role in promoting tumor growth, supported by functional studies on organoids. Authors suggest that global profiling of differentiation state is a promising technique to stratifiy tumor heterogeneity, with potentially might be useful distinguishing thyroid malignancies suitable or not to adjuvant treatment e.g. with radioiodine (RAI) therapy.

    Major comments:

    1. This study largely confirms established facts that 1) PTC due to BRAF driver mutation is a heterogeneous tumour entity and 2) ATC is the most dedifferentiated of all thyroid cancers. Although interesting, observations of a highly variable tissue cell composition including immune cells and the gradual loss of thyroid differentiation markers, in part linked to tumor subclone development featured by altered chromosomal copy numbers, are thus not surprising.
    2. Considering tumor progression, comparison of PTC and ATC should preferably include specimens with the same driver mutation (BRAF or RAS), which is not the case here. This notion should be more clearly explained to readers. An optional improvement would be to conduct similar analyses on an ATC specimen that contains more differentiated PTC tumor portions arguably suggesting that PTC progresses to ATC (by mechanisms that are still largely unexplored).
    3. Comments on findings of lymphocytic infiltration need to be balanced. Although autoimmune thyroid disease in infered a risk factor of developing malignancy it is unlikely that the majority of TCGA samples of PTC is associated with thyroiditis as indicated in Fig. 3 and Suppl Fig. 3. Immune cell abundance may rather reflect the tumor immune microenvironment (TIME).
    4. Some tissue sections seem of quite poor quality either shape-wise of containing rifts e.g. PTC7 in Fig. 3 and PTC2 in Fig. 5. The authors should explain whether and how this might influence analysis.
    5. The experiment on mouse ESC/organoids (Fig. 6H-J) does not show much of an expected enhanced thyroid progenitor cell proliferation after induction of the mutant Braf allele by tamoxifen, which raises doupt whether the subsequent promoted growth by fibronectin at all is oncogene-related. This differs from the impact of BrafCA activation along with mouse thyroid development in vivo (Schoultz et al iScience 2023 PMID: 37534159). In the same experimental setup, it appears that mutant Braf prevents follicle formation (Fig. 6I). A control experiment investigating the influence of fibronectin in the absence of oncogene activation should be conclusive. The effect of Braf and fibronectin on thyroid organoid structure and function should be better explained, if necessary based on complementary experiments, and discussed in relation to the claimed association of fibronectin expression to "...low amounts of thyroid differentiation markers...) and "...loss of epithelial structure (PTC7, Figure 6E)." in the previous section of Results.
    6. Concerning EMT profiling (Supplementary Fig. 7B) , there is a great similarity of ATC tumor cells and fibroblasts, and as stated in the text the malignant status of the former is indicated by chromosomal aberrations (refering to Suppl fig. 6). However, looking at Suppl. Fig. 7B it is evident that fractions of cells identified as fibroblasts express TG and TSHR suggesting mismatch. How was this comparison done in order to exclude mismatch? Is there no other profiled markers that distinguish cancer cells from stromal cells that can support conclusions? In the same figure, it appears there are no clear differences in EMT marker expression among PTC samples regardless of differentiation state, suggesting that the gradual loss of thyroid differentiation in PTC tumor cells and EMT are not parallel and potentially linked phenomena? Please clarify this dissociation of results. Is is possible that refocusing on other EMT markers than the top 10, of which almost all concerns various collagen genes, might better reveal partial EMT in PTCs?
    7. According to Suppl. Table 1, the ATC2 tumor does not harbor any mutations. What about chromosomal aberrations, was that included in analysis? Considering previous consistent reports of a high mutation burden in ATC, if not supported by other data (clinical, pathological) the diagnosis might be questioned for this particular case included in multiple analyses of the present study.

    Minor comments:

    • The logical order of presentation of Results might benefit from first presenting specific PTC data following by ATC dito. I´m thinking of swapping the section of EMT in ATC to end of Results.
    • Methods paragraph "Mouse ESC-derived thyroid organoids experiments" (starting with "ccc") seems to be missing some essential information.

    Significance

    The study confirms at single cell level the fundamental difference of PTC and ATC that is evident clinically and biologically, but does not address the intriguing issue how ATC may progress from PTC.

    Tumor heterogeneity of BRAFV600E-driven PTC in terms of dedifferentiation of functional parameters, which are of potential clinical relevance, is well documented.

    Reviewer expertise: thyroid development, thyroid cell and tumor biology, superficial knowledge in scRNAseq analysis

  3. Review Commons

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

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    This study is well designed with a rational sample collection strategy. The authors collected PTC and ATC tissue samples for snRNA and spRNA sequencing, clearly characterizing tumor heterogeneity. Using representative thyroid differentiation markers (TSHR, TPO, TG, NIS), they distinguished different differentiation states of PTC and ATC and further validated the role of FN1 in organoid models. However, the manuscript is largely descriptive in nature, and several key issues remain to be addressed.

    Major Comments:

    1.Differential gene/pathway analysis across epithelial clusters: What are the differential genes or pathways among the epithelial clusters? Without CCA/Harmony integration, do the tumor subgroups show distinct differences? In addition, I suggest applying NMF or hdWGCNA to identify shared modules and test whether ATC and PTC harbor overlapping regulatory modules. 2.Validation of TSHR/TPO-based subgrouping: While the TSHR/TPO grouping appears appropriate for stratification at the single-cell level, it is necessary to exclude sequencing depth as a confounding factor. Should validate the existence of these subpopulations using mIHC/IF on corresponding samples. 3.Impact of mutational differences on conclusions: According to Supplementary Table 1, almost all PTC cases carried BRAF mutations, whereas four ATC patients harbored no BRAF mutation. Could this difference influence the conclusions of the study? Although the authors briefly mention this in the Discussion, a more thorough clarification is warranted. 4.The statement "Myeloid and T cells also grouped in specific clusters" seems descriptive. Is this clustering biologically meaningful? Please elaborate.

    Minor Comments:

    In Figure 2C, the "Epith TSHR-" population resembles myeloid cells. Could the authors clarify why this is the case? For the correlation analysis in Figure 2C, were highly variable genes or all genes used?

    Significance

    This study provides a comprehensive single-nucleus and spatial transcriptomic atlas of papillary and anaplastic thyroid carcinomas. Its strengths include well-designed sample collection, high-resolution profiling of tumor heterogeneity, and validation of FN1 function. By stratifying malignant cells with thyroid differentiation markers (TSHR, TPO, TG, NIS), the authors delineate differentiation states and highlight mechanisms of progression from PTC to ATC. However, the study remains mainly descriptive, and additional analyses of gene modules, pathway regulation would increase its conceptual depth. The findings will interest researchers in thyroid cancer, tumor heterogeneity, and the single-cell/spatial genomics field, with potential relevance for translational oncology.

    Field of expertise: thyroid cancer biology, single-cell and spatial transcriptomics.

  4. Version published to 10.1101/2024.02.15.580495 on bioRxiv