The intratumor microbiome predicts prognosis across gender and subtypes in papillary thyroid carcinoma.

While the intratumor microbiome has become increasingly implicated in cancer development, the microbial landscape of papillary thyroid carcinoma (PTC) is essentially uninvestigated. PTC is characterized by varied prognosis between gender and cancer subtype, but the cause for gender and subtype-based dissimilarities is unclear. Women are more frequently diagnosed with PTC, while men suffer more advanced-staged PTC. In addition, tall cell variants are more aggressive than classical and follicular variants of PTC. We hypothesized that intratumor microbiome composition distinctly alters the immune landscape and predicts clinical outcome between PTC subtypes and between patient genders. Raw whole-transcriptome RNA-sequencing, Level 3 normalized mRNA expression read counts, and DNA methylation 450 k sequencing data for untreated, nonirradiated tumor, and adjacent normal tissue were downloaded from the Genomic Data Commons (GDC) legacy archive for 563 thyroid carcinoma patients. Microbe counts were extracted using Pathoscope 2.0 software. We correlated microbe abundance to clinical variables and immune-associated gene expression. Gene-set enrichment, mutation, and methylation analyses were conducted to correlate microbe abundance to characterize microbes' roles. Overall, PTC tumor tissue significantly lacked microbes that are populated in adjacent normal tissue, which suggests presence of microbes may be critical in controlling immune cell expression and regulating immune and cancer pathways to mitigate cancer growth. In contrast, we also found that microbes distinctly abundant in tall cell and male patient cohorts were also correlated with higher mutation expression and methylation of tumor suppressors. Microbe dysbiosis in specific PTC types may explain observable differences in PTC progression and pathogenesis. These microbes provide a basis for developing specialized prebiotic and probiotic treatments for varied PTC tumors.

Introduction

Thyroid cancer is the fifth most common cancer in women, and it is estimated that over 12,000 men and 40,000 women will be diagnosed with thyroid cancer in 2021 [1]. Advanced-stage thyroid cancer presents itself more frequently in men than women, but the cause of this disparity is unclear [2], [3]. Despite advancements in diagnostic technology for thyroid cancer, precision medicine and nonsurgical treatment options for papillary thyroid carcinoma (PTC) and its most common subtypes, classical (CPTC), follicular variant (FVPTC), and the more aggressive tall cell papillary thyroid carcinoma (TCPTC) are severely lacking [4], [5], [6], [7], [8], [9].

A recent study identified the presence of microbes in blood and thyroid cancer tissue, suggesting that there exists a predictable and significant abundance of microbes within PTC issues [10]. Furthermore, there is mounting evidence that the intratumor microbiome is implicated in cancer progression and mitigation through their metabolites. For example, the absence of folate-producing Bifidobacterium and Lactobacillus was correlated with hypomethylation at regions of the p53 gene in colorectal cancer [11]. This research suggests that tissue-specific microbes may play a tumor suppressive role, specifically by causing or preventing the repair of DNA damage. Past microbiome studies have used The Cancer Genome Atlas (TCGA) data and machine learning to identify diagnostic microbial signatures in various cancers but have not elucidated how microbe composition varies in PTC and how microbes interact with pathways to influence cancer progression [10]. We hypothesize that the intratumor microbiome in thyroid carcinoma plays a critical role in determining the varying immune landscape across different types of PTC.

In this study, we characterized the abundance of microbes in PTC using RNA sequencing data obtained from TCGA. We identified microbes differentially abundant between males and females and between different PTC subtypes. To determine the mechanisms by which specific microbes may influence the immune landscape, we correlated microbe abundance to clinical variables, immune-associated (IA) gene expression, immune-cell expression, and cancer and immunological pathways. Finally, we correlated microbe abundance with known copy number variations (CNV) and mutations and methylation sites (Fig. 1). With this study, we aim to characterize the thyroid carcinoma intratumor microbiome and investigate microbiome’s unique clinical significance in thyroid carcinoma prognosis.

Fig. 1

Schematic of project workflow and analyses. After obtaining sequencing data from The Cancer Genome Atlas (TCGA), we used the Pathoscope software to extract microbial read counts. We normalized the microbial counts using Aitchison’s log transformation. 1) We performed differential abundance to identify which microbes were significant differentially abundant in tumor vs. normal samples. The overlap in differentially abundant microbial species between various cohorts is demonstrated through venn diagrams: male versus female thyroid carcinoma and classical papillary thyroid carcinoma (CPTC) versus follicular variant papillary thyroid carcinoma (FVPTC) versus tall cell papillary thyroid carcinoma (TCPTC). 2) Three methods of contamination correction were used to identify potential contaminants among the significantly differentially abundant microbes we identified in Step 1. We performed correction by sequencing plate, sequencingdata, and abundance in one individual sample vs. all samples. We also compiled a list of common contaminants from previously published research and eliminated these microbes from further study. 3) MACIS (distant Metastasis, patient Age, Completeness of resection, local Invasion, and tumor Size) score and clinical variable analysis of differentially abundant microbes were performed using the Kruskal-Wallis statistical test. Clinical variable data (pathologic stage, TNM stage) for each sample were obtained from the Broad Institute GDAC Firehose database. 4) Correlation of microbial abundance to immune-associated gene expression and immune cell infiltration: We performed gene differential expression analysis on Level 3 normalized mRNA expression read count data obtained from Genomic Data Commons legacy archive/TCGA. We then correlated microbe abundance with dysregulated immune-associated gene expression using Kruskal-Wallis and visualized and categorized the data using Cytoscape Reactome FIViz software. We used CIBERSORTx software to obtain immune cell expression correlated with microbe abundance. 5) We used Gene Set Enrichment Analysis to identify immunological and oncogenic signatures correlated with microbial abundance. 6) We used the Repeated Evaluation of Variables conditionAL Entropy and Redundancy (REVEALER) to identify chromosomal amplifications, deletions, and mutations correlated with microbial abundance. 7) We used DNA methylation 450k sequencing data and a modified workflow of methylationArrayAnalysis from Bioconductor to identify differentially methylated sites in thyroid carcinoma. We correlated differentially methylated sites with microbial abundance using the Kruskal-Wallis statistical test.

Materials and methods

Data acquisition

Raw whole-transcriptome RNA-sequencing, Level 3 normalized mRNA expression read counts, and DNA methylation 450 k sequencing data for untreated, nonirradiated tumor and adjacent normal tissue were downloaded from the Genomic Data Commons (GDC) legacy archive for 563 thyroid carcinoma (THCA) patients (354 CPTC, 101 FVPTC, 35 TCPTC, 135 male, 366 female tumor samples) [12], [13]. The GDC provides web-based access to data and metadata from cancer genomics studies aggregated in TCGA, and The Broad Institute developed a suite of tools and pipelines to assist in aggregating and analyzing various types of large-scale genomic and proteomic data. Clinical information for all patients was downloaded from Broad GDAC Firehose [14]. Genomic alteration information for each patient was obtained from the last analysis report (2016) of the Broad Institute TCGA Genome Data Analysis Center [15].

Extraction of microbial reads

RNA-sequencing data were filtered for bacterial reads via direct alignment to a reference library of bacterial sequences using the Pathoscope 2.0 software and the NCBI nucleotide database [16]. To account for the inherent compositionality of microbe data, we performed a ratio transformation of the microbe abundance data. Ratio transformations capture the relationships between samples in the data, regardless of whether the data are counts or proportions. Log-ratios make the dataset linearly related, so that log-ratio abundances of features can be compared to other samples in the data [17], [18]. We used the R “compositions” package to perform the variance-stabilizing Aitchison’s log transformation to eliminate negative correlation bias and account for large amounts of variability between samples.

Evaluation of contamination

Contamination of RNA-sequencing data may originate from appliance contamination or sequencing errors [19]. In order to account for contamination, we employed multiple contamination correction methods based on TCGA sampling and sequencing methodology [20]. First, scatterplots of microbial abundance versus date of sequencing were produced; potential contaminants were determined based on unrepresentative overabundance of bacteria, or overabundance solely across three (or less) consecutive dates. Second, boxplots of microbe abundance versus sequencing plates were produced to identify unrepresentative overabundance of bacteria, or overabundance of bacteria solely in samples sequenced on the same plate. Third, total abundance across all microbes was compared against the abundance of each individual microbe, per patient. As the total number of read counts for each sample increases, it is reasonable to assume that the abundance of any specific microbe will also increase if the microbe is not a contaminant. We identified potential contaminants in scatter plots of individual microbe abundance versus total abundance with a slope of zero (margin around 0 of ±0.1)

Differential abundance between patient cohorts and correlation to MACIS (distant metastasis, patient age, completeness of resection, local Invasion, and tumor Size) score and clinical variables

The Kruskal-Wallis test was used to determine differential abundance (p < 0.05) and correlations between microbe abundance and lymph node metastasis, MACIS score, neoplasm cancer status, pathologic stage, pathologic tumor-node-metastasis stage, and residual tumor (p < 0.05).

Correlation of microbial abundance to immune-associated elements, copy number variants and mutations, and methylation sites

Dysregulated immune-associated (IA) genes are defined as genes with altered expression between tumor and normal tissue. The genes were determined using edgeR analysis of mRNA expression data obtained from TCGA (FDR < 0.05 and |log fold change| > 1) [13]. EdgeR offers highly accurate statistical methods for multigroup experiments [21], [22]. Dysregulated IA genes were clustered using the Cytoscape Reactome FIViz software. The Kruskal-Wallis test was used to determine the presence of correlation between differentially abundant microbes and dysregulated IA genes (p < 0.01). The p-value cut-off was reduced to 0.01 to focus on the most statistically significant correlations. The threshold was determined by p-value in order to evaluate effect size and assess the probability of the results supporting the hypothesis.

Using the software tool CIBERSORTx and mRNA expression data, we used gene expression profiles to estimate immune infiltration [23]. The input was a matrix of gene expression read counts by THCA samples, compiled from mRNA sequencing data [13]. CIBERSORTx software returned expression values for 22 immune cell types.

We conducted Gene Set Enrichment Analysis (GSEA) to correlate microbe abundance to known cancer pathways and immunological signatures. The GSEA software and signatures from the C2, C6, and C7 collections were downloaded from the Molecular Signatures Database [24]. Statistically significantly (p < 0.05) enriched signatures were plotted.

The surface-level trends of mutation presence were analyzed by calculating the percentage of patients with each mutation, indicated by a binary value per mutation. The GDAC files were compiled into input files for the REVEALER (repeated evaluation of variables conditional entropy and redundancy) algorithm, which identifies sets of specific CNVs, including chromosomal deletions, amplifications, and alterations, and mutations that are most likely implicated in changes to the target expression profile. The target profile was identified as the expression of a single immune-associated gene. The REVEALER algorithm runs in multiple iterations in order to identify the most prominent genomic alterations. For our study, we set the maximum number of iterations to three. The algorithm also allows for the use of a seed, or a particular mutation of CNV gain or loss event that may account for target activity. Because we did not know the individual genetic alterations that were responsible for each IA genes’ dysregulation, the seed was set to null. Significant correlations were indicated by p < 0.05 and CIC > 0.15 [25].

Using a modified workflow of methylationArrayAnalysis from Bioconductor, we converted B-values to M-values and then performed probe-wise differential methylation analysis on the matrix of M-values in limma to obtain t-statistics and p-values for each CpG site [26], [27]. The Kruskal-Wallis test was used to correlate microbe abundance to the extent of methylation.

Validation with external dataset

77 CPTC and 48 FVPTC samples were obtained from the European Nucleotide Archive (Project ID: PRJEB11591) [28]. Pathoscope 2.0 was used to align and extract microbial reads. Microbe abundance was plotted to compare relative abundance between tumor samples from TCGA and ENA.

Results

Microbial landscapes of tumor tissue between genders are distinct

Principal component analysis (PCA) was conducted to measure heterogeneity in microbe species between tumor and normal samples, tumor samples of varied subtype, and tumor samples from different genders (Fig. 2). Overall, the principal component plots show major clustering and indicate that the microbe populations between the tumor samples are relatively homogenous. However, the third PCA plot contains a clustering of female patient tumor samples in the top right corner which makes clear that the microbial landscapes of tumor tissues from the different genders may be distinct. We found the microbes that contributed most to the top-right cluster to be cyanobacteria species.

Fig. 2

Principal component analysis plots to visualize microbial landscape in the following comparisons: Tumor versus normal tissue, classical papillary thyroid carcinoma versus follicular variant papillary thyroid carcinoma versus tall cell papillary thyroid carcinoma, and males versus females.

Contamination correction

To verify that our measurement of intratumor microbial abundance was consistent across samples and not due to contamination, we conducted three methods of contamination correction. We identified 17 and 21 potential contaminants using correction by sequencing date and plate respectively (Fig. 3A, B). We compared individual microbial abundance to the abundance of all microbes to identify aberrant abundance potential contaminants (Fig. 3C). Finally, we compiled a list of 175 commonly known contaminants from studies of bacterial contamination from reagents and of the hospital microbiome [29], [30]. In all, we identified 202 potential contaminants to be eliminated from further investigation (Supplementary Table 1).

Fig. 3

Contamination correction. (A) Potential microbial contaminants were identified by producing scatterplots of individual microbe abundance versus date of sequencing of sample. Independent clusters of high abundance across a few dates indicates that the microbe is not consistently differentially abundant across all samples. (B) Potential microbial contaminants were identified by producing boxplots of individual microbe abundance versus sequencing plate. Singular boxplots of high abundance in one sequencing plate indicates that the microbe is not consistently differentially abundant across all samples. Instead, the microbe appears to be differentially abundant due to contamination in one sequencing plate. (C) Read counts of total abundance across all microbes was compared against the abundance of each individual microbe, per patient. Potential contaminants were identified in scatter plots with a slope of zero. This method elicited zero potential contaminants, as all plots had a slope with magnitude significantly greater than zero.

Differential abundance of intratumor microbes

We found 45 microbes in CPTC, 34 in FVPTC, and 33 in TCPTC to be differentially abundant between tumor and normal tissue. We found 33 microbes in male samples and 49 microbes in female samples to be differentially abundant between tumor and normal tissue. Of all the statistically significant correlations, we focused on investigating the strongest correlations (p < 0.05) of microbes that have been previously discovered and studied in the context of cancer and inflammatory disease. Micrococcus luteus, Frankia sp., Anabaena sp. K119, and uncultured Gammaproteobacteria bacterium were all similarly overabundant in normal tissue in CPTC, FVPTC, and TCPTC (Fig. 4A). Trueperella pyogenes and Stenotrophomonas maltophilia K279a were similarly dysregulated in CPTC and FVPTC (Fig. 4B). We also identified the strongest correlations of microbes uniquely differentially abundant in each subtype cohort (Fig. 4C). Rhodococcus fascians D188 was overabundant in normal samples in CPTC. Acinetobacter baumannii AB0057 was overabundant in normal samples in FVPTC. Bradyrhizobium sp. BTAi1 was overabundant in normal samples in TCPTC.

Fig. 4

Differential abundance of microbes. (A) Microbes identified to be similarly differentially abundant across all three papillary thyroid carcinoma (PTC) subtypes. Yellow corresponds to classical PTC (CPTC) data, green corresponds to follicular variant PTC (FVPTC) data, and red corresponds to tall cell PTC (TCPTC) data. (B) Microbes identified to be similarly differentially abundant in CPTC and FVPTC. (C) Microbes identified to be uniquely differentially abundant in CPTC, FVPTC, and TCPTC. (D) Synechococcus sp. CC9311 was significantly differentially abundant in both male and female thyroid tumor samples. However, the microbe was absent in male tumor tissue versus normal tissue but overabundant in female tumor tissue versus normal tissue. (E) Venn diagram schematics summarizing microbes found to be differentially abundant in patient cohorts. The number indicates the number of most significant correlations, and the down arrow indicates that the microbe was sparse in tumor samples compared to normal samples. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In males, Synechococcus sp. CC9311 was overabundant in the normal samples while in females, it was overabundant in tumor samples (Fig. 4D). Overall, normal tissue had significantly more differentially abundant microbes compared to tumor tissue samples (Fig. 4E).

Microbe abundance correlation with MACIS (distant metastasis, patient age, completeness of resection, local invasion, and tumor size) score and clinical variables

The majority of significant correlations between microbe abundance and clinical variables were negative, consistent with the finding that tumor tissue contained lower microbe abundance compared to adjacent normal tissue. Frankia sp. and uncultured Gammaproteobacteria bacterium, which are abundant in all PTC normal tissue samples, are correlated with lower MACIS score, while Bradyrhizobium sp. BTAi1, which is uniquely abundant in TCPTC normal tissue, is correlated with higher MACIS score. Frankia sp. and Anabaena sp. K119, both of which are overabundant in normal tissue samples of all PTC subtypes, are negatively correlated with Pathologic M stage, while Stenotrophomonas maltophilia, dysregulated in only CPTC and FVPTC, is positively correlated with pathologic M stage (Fig. 5).

Fig. 5

Correlation to clinical variables. Select boxplots of microbes significantly positively and negatively correlated with MACIS (distant Metastasis, patient Age, Completeness of resection, local Invasion, and tumor Size) prognostic score and pathological M stage. M0 indicates that the cancer has not spread to other parts of the body, while M1 indicates spread of cancer/distant metastases. All boxplots were produced using the Kruskal-Wallis test.

Microbe abundance correlation to immune-associated genes

Dysregulated IA genes found by comparing gender and subtypes cohorts to normal samples were inputted into Reactome FIViz software and clustered into gene modules (Fig. 6A). We found that the greatest number of dysregulated IA genes with the greatest magnitude was correlated to Module 2 pathways in CPTC and FVPTC. However, in TCPTC, IA genes that were correlated to Module 3 pathways were the most dysregulated (Fig. 6B). IA gene dysregulation was comparable between male and female patient cohorts. We determined that CPTC had the greatest absolute number of correlations between microbe dysbiosis and IA gene dysregulation (Fig. 6C).

Fig. 6

Correlation to immune dysregulation (A) Dysregulated immune-associated (IA) genes of each subtype and gender cohort were clustered into modules by Cytoscape Reactome FIViz software. Module 0 corresponded to cytokine expression and signaling-related pathways. Modules 1, 2, and 3 describe signaling and cancer pathways that were grouped together according to strongest pairwise relationships of functional interactions, or interactions in which two proteins are involved in the same reaction as input, catalyst, activator, or inhibitor. (B) Pathway enrichment analysis using Reactome FI was performed to identify pathways most correlated with IA gene modules (p < 0.05 and FDR < 0.05) (C) Barplot summarizing the number of microbes (differentially abundant in each patient cohort) correlated with IA genes in each module. Differentially abundant microbes were correlated with dysregulated IA gene expression using the Kruskal-Wallis test. (D) Radar plots visualizing levels of immune cell infiltration, measured using CIBERSORTx software, in subtype and gender cohorts.

Microbe abundance correlation to immune cell infiltration

We found FVPTC samples to be strongly correlated to the expression of resting CD4 memory T-cells. In tumor versus normal samples, we observed higher expression of immune cells in normal tissue. Notably, follicular helper T-cells, resting dendritic cells, and activated mast cells were highly expressed in normal tissue in male patients, while naive B cells and resting CD4 memory T cells were highly expressed in normal tissue in female patients (Fig. 6D).

Microbe abundance correlation to cancer and immunologic pathways

By using GSEA to correlate microbe abundance to known cancer and immunologic pathways and signatures, we found that decreased abundance of microbes in CPTC and FVPTC tumor tissue was correlated with decreased expression of EIF4E_UP and SRC_UP.V1_UP oncogenic signatures. In TCPTC, we found that abundance of Anabaena sp. K119 in tumor tissue was correlated with the upregulation of the EIF4E_UP signature, which is consistent with the more aggressive nature of TCPTC. Among the gender cohorts, the majority of significantly correlated pathways were positively enriched with microbe abundance (Fig. 7A). Although Synechococcus sp. CC9311 was overabundant in female but not in male tumor tissue, all overlapping pathways were similarly dysregulated in both gender cohorts (Fig. 7B). We found that the majority of most enriched pathways uniquely dysregulated in the different subtypes were related to cell growth (Fig. 7C). In males, the greatest number of uniquely dysregulated pathways belonged to the tumor suppression-related group, and in females, the greatest number of uniquely dysregulated pathways belonged to the DNA checkpoint and damage-related group (Fig. 7D).

Fig. 7

Gene Set Enrichment Analysis (GSEA). (A) Heatmap of GSEA correlations associating microbe abundance to dysregulated pathways in the three subtype cohorts: classical papillary thyroid carcinoma, follicular variant papillary thyroid carcinoma, tall cell papillary thyroid carcinoma (p < 0.05). Uncultured gammaproteobacteria bacterium and Anabaena sp. K119 were both less abundant in tumor samples. Red corresponds to upregulated pathways, and blue corresponds to downregulated pathways. (B) Heatmap of GSEA correlations associating microbe abundance to dysregulated pathways in both male and female patient cohorts (p < 0.05). The innermost ring corresponds to enriched pathways in male tumor samples, and the middle ring corresponds to enriched pathways in female tumor samples. The outermost ring specifies the microbe correlated to each dysregulated pathway. All microbes, except for Synechococcus sp. CC9311, were less abundant in tumor samples versus normal samples in both cohorts. In male samples, Synechococcus sp. CC9311 was less abundant in tumor tissue than in normal tissue, while in females, Synechococcus sp. CC9311 was more abundant in tumor tissue than in normal tissue. (C) Barplots visualizing most enriched pathways uniquely dysregulated (p < 0.05) in papillary thyroid cancer subtype cohorts. (D) Barplots visualizing most enriched pathways uniquely dysregulated (p < 0.05) in gender cohorts. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Microbe abundance correlation to copy number variations and mutations

To further investigate how microbes may regulate gene expression, we used the REVEALER algorithm to correlate microbe abundance to CNV and mutation rates (Fig. 8A). Microbe abundance in FVPTC was correlated with the greatest number of chromosomal alterations and mutations. Between the male and female patient cohorts, microbe abundance in males was more frequently and strongly correlated to a greater number of CNV. The Kruskal-Wallis test confirmed that microbes dysregulated in TCPTC and male patients were correlated with higher expression of the BRAF V600E mutation, which is highly implicated in PTC progression (Fig. 8B).

Fig. 8

(A) Associations between genomic alterations and microbe abundance of differentially abundant microbes in each subtype and gender cohort represented by mutation heatmaps. Red on the left side of the heatmap indicates that the microbe is overabundant, while red on the right side of the heatmap indicates that the microbe is underabundant. The following rows below the heatmap specify CNV type/genetic variants (alteration, deletion, mutation) across all samples. Each bar represents the presence or absence of a CNV in one sample. Correlations were measured using the Repeated Evaluation of Variables conditionAL Entropy and Redundancy (REVEALER) algorithm (CIC > 0.15, p < 0.05) (B) Boxplots correlating microbe abundance to BRAF V600E gene mutation, which is highly implicated in PTC. All boxplots were produced using the Kruskal-Wallis statistical test, and mutation expression data was obtained from the Broad Firehose database. Mutation expression was converted to binary and categorical variables of “HIGH” and “LOW” in order to conduct the Kruskal-Wallis statistical test. The median expression level for the BRAF_Mutation across all samples was calculated. Values lower than the median value were assigned an expression level of “LOW,” and values above the median value were assigned an expression level of “HIGH.” The x-axis refers to mutation type and level, and the y-axis refers to numerical abundance values of each microbe. The red plot indicates that differential Bradyrhizobium abundance in TCPTC tissue is correlated with higher expression of the BRAF mutation. The blue plots indicate that lower microbe abundance in male tumor tissue is correlated with higher expression of the BRAF mutation. (C) Boxplots correlating microbe abundance with extent of methylation at CpG sites. In this analysis, microbe abundance was converted to a binary and categorical variable of “HIGH” and “LOW” in order to conduct the Kruskal-Wallis test. Microbe abundance was converted to binary and categorical variables of “HIGH” and “LOW” in order to conduct the Kruskal-Wallis statistical test. The median abundance level for each microbe across all samples was calculated. Values lower than the median value were assigned an abundance level of “LOW,” and values above the median value were assigned an abundance level of “HIGH.” The x-axis refers to microbe abundance level, and the y-axis refers to numerical extent of methylate at a particular CpG site. Overall, lower microbe abundance, as observed in PTC tissue, was correlated with higher levels of methylation of cell-cycle related genes and tumor suppressors, including NEURL1B (cg07485775), POLE (cg27197524), and SRCIN1 (cg14122138). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Microbe abundance correlation to gene methylation

We found the most differentially methylated CpG sites between tumor and normal samples at promoter regions on chromosomes with the greatest number of copy number alterations that we identified from our REVEALER analysis (p < 0.05). The most significantly methylated sites occurred on chromosomes 1 and 17 and specifically at cell cycle related genes. For the majority of significant correlations, including 55 in CPTC, 17 in FVPTC, and 5 in TCPTC, lower microbe abundance was correlated with greater extent of methylation at known tumor suppressor genes including NEURL1B, POLE, and SRCIN1 [31], [32], [33] (Fig. 8C).

Validation with ENA dataset

All microbes identified from TCGA data and discussed in the manuscript were also present in samples from the ENA dataset (Fig. 9). Microbe abundance was comparable between the two groups of tumor samples from each dataset. The ENA dataset did not contain normal samples, and further in vitro and in vivo studies are necessary to validate our findings.

Fig. 9

Boxplots comparing microbe abundance between tumor samples from The Cancer Genome Atlas (TCGA) and a separate dataset with CPTC and FVPTC data obtained from the European Nucleotide Archive (ENA) (Project ID: PRJEB11591). “Normal-TCGA” refers to adjacent normal thyroid tissue matched to cancer tissue for corresponding patients.

Discussion

While recent studies have shown the influence of the intratumor microbiome on cancer development, the intratumor microbiome of papillary thyroid carcinoma has remained largely unexplored. Our results indicate that distinct differences in microbe abundance between tumor and normal samples can explain differences in prognosis between patients of varied PTC subtype and gender. PCA results revealed that Chroococcidiopsis sp. and other cyanobacteria species abundant in female tumor and normal samples produce great amounts of phycobiliproteins, which possess high antioxidant capacity that prevents cancer cell proliferation [34]. Overall, PCA showed major clustering, which suggests that intratumor microbial populations across all patient samples are fairly predictive. This highlights the potential of using microbes as a diagnostic or treatment measure for PTC. Micrococcus luteus, which is overabundant in TCPTC tissue, has been previously shown to cause severe infection in immunocompromised patients [35] as well as anaphylactoid reactions and early death in mice [36]. Presence of M. luteus in TCPTC tissue can explain the aggressiveness of TCPTC. On the other hand, Frankia sp., which has exhibited antimicrobial activity against phytopathogens, may act similarly against cancerous molecules in humans [37]. Among the species differentially abundant in only CPTC and FVPTC, it was recently found that presence of intratumor Trueperella pyogenes was significantly correlated with lower PSA levels in prostate cancer patients [38]. In addition, microbiota dysbiosis in CPTC was correlated with the greatest number of IA genes, potentially explaining why CPTC patients are low-risk, have superior prognosis, and respond better to immunotherapy compared to other PTC subtypes that are less immune-associated [39]. Bradyrhizobium sp. BTAi1, which we found to be uniquely dysregulated in TCPTC tissue and positively correlated with MACIS score, has been shown to be consistently overabundant in tumor tissue, especially in cervical cancer [40]. Significant difference of Synechococcus sp. CC9311 abundance based on gender could explain differences in PTC progression between males and females.

Overall, the trend of overabundance of microbe species in normal samples strongly suggests that intratumor microbes in CPTC and FVPTC are critical for stimulating tumor-suppressive activity. Our analyses support the anti-tumor role of microbes by showing that microbes may regulate CD4 and helper T-cells to defend against tumor growth in FVPTC. However, the ubiquitous overexpression of immune cells may contribute to the aggressiveness of PTC by causing systemic inflammation, thereby promoting tumor growth, as seen in the male patient cohort and in the case of M. luteus and Bradyrhizobium sp. BTAi1 abundance in TCPTC [41]. A recently published study found that microbiome-derived ligands and metabolites act directly on intestinal immune cells and travel to remote tissues via the systemic circulatory system to modulate immune cell expression. Some microbe species were found to prevent disease-promoting activity and limit the expansion of iNKT cells by producing sphingolipids. Other species were also found to promote CD4+ T cell differentiation and to balance Th1 and Th2 populations. On the other hand, dysbiosis of specific microbes and the overexpression of TGF-β resulted in increased regulatory T-cell expression and expansion of proinflammatory Th17 cells [42]. These results serve to validate our own findings that microbes can regulate immune cell expression differently to mitigate or aggravate disease.

In order to clarify microbes’ putative roles established by past literature, we conducted pathway, mutation, and methylation correlational analyses to elucidate mechanisms by which these microbes may regulate components of the immune landscape. We found microbes abundant in CPTC and FVPTC to be negatively correlated with oncogenic signatures, which is consistent with the two subtypes being less deadly than TCPTC and poorly differentiated thyroid cancer. We found microbe abundance in males to be negatively correlated with tumor suppressive pathways, which may explain worse prognosis in males, while microbe abundance in females was positively correlated with DNA damage. We found microbe abundance in males to correlate with a greater number of chromosomal alterations compared to female patient samples. Most significant methylation sites were found at cell-cycle and tumor suppressor related genes, and lower microbe abundance was correlated to higher levels of mutation, which suggests that higher microbe abundance, as observed in normal samples, may be necessary to properly regulate cell-cycle related genes.

To the best of our knowledge, we are the first to provide the largest comprehensive profiling of papillary thyroid carcinoma samples to date and show how the abundance of specific microbes may contribute to pathogenesis in different subtypes of PTC. We used numerous computational tools to rigorously collect and correlate microbial abundance data to cancer stage and severity, immune cell expression, cancer pathways, known cancer-related copy number variations and mutations, and methylation. Thorough integration and interpretation of our data suggests how specific microbes may work to mitigate tumor growth in CPTC and FVPTC or aggravate tumor growth in TCPTC. While we cannot definitively conclude if microbes directly cause differences in clinical outcome of thyroid cancer or if microbe abundance is a result of the hypoxic environment created by the tumor, previously published studies on the intratumor microbiome strongly suggest that microbes more likely contribute to pathogenesis than they are a mere result of it. One study found that differences in the gut microbiome composition directly related to p53’s tumor-suppressive effectiveness [43]. Another study found that microbiota contributed to colorectal cancer pathogenesis and genomic instability by secreting chemokines that recruit immunosuppressive myeloid-derived suppressor cells, tumor-associated macrophages, and tumor-associated neutrophils and by inhibiting the activity of cytotoxic activity of natural killer cells. The study also found some microbe species, including Clostridia, to play an anti-tumor role in colorectal cancer. These publications provide support to our own findings that specific microbe species can play an anti-tumor or pro-tumor role in thyroid cancer [44]. Nonetheless, since the results of our study are limited to correlative data, controlled in vitro and in vivo experiments are necessary to confirm our computational findings on the mechanisms by which specific microbes contribute to tumor growth or suppression in thyroid cancer subtypes.

Conclusion

In conclusion, our study significantly advances the understanding of the papillary thyroid carcinoma microbiome composition, its heterogeneity across samples of different PTC subtypes and patient gender, and its relationship with clinical and immunologic variables. Microbiota abundant in CPTC (Frankia sp. and Trueperella pyogenes) are highly correlated with immune-associated genes, suggesting that metabolite-gene interaction may result in better prognosis of CPTC compared to the more aggressive TCPTC phenotype. M. luteus and Bradyrhizobium sp. BTAi1, which are abundant in TCPTC, are positively correlated with uncontrolled cell growth pathways and oncogenic pathways (BIOCARTA_CTLA4_PATHWAY, p53 instability). Our results suggest that the overabundance of M. luteus and Bradyrhizobium sp. BTAi1 in TCPTC and male samples causes dysregulation of these critical pathways, leading to high levels of mutation and causing greater cancer severity (strong positive correlation to MACIS score). While validation analysis showed similar levels of abundance of the same microbes that we identified from TCGA data, deeper sequencing, spatial examination of tissue, and controlled experiments are necessary to confirm our findings on the effects of PTC microbiota and microbial heterogeneity between samples. In vitro and in vivo experiments must also be conducted to elucidate how varied physiological conditions of the tumor microenvironment could possibly determine and regulate microbe abundance and activity and how microbes can be utilized as therapeutic agents and prebiotic or probiotic treatment options for PTC.

Funding

UCSD Academic Senate, Grant, RG096651 to WMO.

CRediT authorship contribution statement

Aditi Gnanasekar: Methodology, Formal analysis, Investigation, Writing - original draft, Writing - review & editing, Visualization. Grant Castaneda: Formal analysis, Investigation, Writing - review & editing. Anjali Iyangar: Formal analysis, Investigation, Writing - review & editing. Shruti Magesh: Formal analysis, Investigation, Writing - review & editing. Daisy Perez: Formal analysis, Investigation, Writing - review & editing. Jaideep Chakladar: Methodology, Writing - review & editing. Wei Tse Li: Methodology, Writing - review & editing. Michael Bouvet: Writing - review & editing. Eric Y. Chang: Writing - review & editing. Weg M. Ongkeko: Conceptualization, Resources, Writing - review & editing, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.