Gut microbiota changes and its potential relations with thyroid carcinoma.

Introduction: Emerging evidence suggests that the essence of life is the ecological balance of the neural, endocrine, metabolic, microbial, and immune systems. Gut microbiota have been implicated as an important factor affecting thyroid homeostasis.
Objectives: This study aims to explore the relationship between gut microbiota and the development of thyroid carcinoma.
Methods: Stool samples were collected from 90 thyroid carcinoma patients (TCs) and 90 healthy controls (HCs). Microbiota were analyzed using 16S ribosomal RNA gene sequencing. A cross-sectional study of an exploratory cohort of 60 TCs and 60 HCs was conducted. The gut microbiota signature of TCs was established by LEfSe, stepwise logistic regression, lasso regression, and random forest model analysis. An independent cohort of 30 TCs and 30 HCs was used to validate the findings. Functional prediction was achieved using Tax4Fun and PICRUSt2. TC patients were subsequently divided into subgroups to analyze the relationship between microbiota and metastatic lymphadenopathy.
Results: In the exploratory cohorts, TCs had reduced richness and diversity of gut microbiota compared to HCs. No significant difference was found between TCs and HCs on the phylum level, though 70% of TCs had increased levels of Proteobacteria-types based on dominant microbiota typing. A prediction model of 10 genera generated with LEfSe analysis and lasso regression distinguished TCs from HCs with areas under the curves of 0.809 and 0.746 in the exploration and validation cohorts respectively. Functional prediction suggested that the microbial changes observed in TCs resulted in a decline in aminoacyl-tRNA biosynthesis, homologous recombination, mismatch repair, DNA replication, and nucleotide excision repair. A four-genus microbial signature was able to distinguish TC patients with metastatic lymphadenopathy from those without metastatic lymphadenopathy.
Conclusion: Our study shows that thyroid carcinoma patients demonstrate significant changes in gut microbiota, which will help delineate the relationship between gut microbiota and TC pathogenesis.

Introduction

Thyroid carcinoma is the most common head and neck endocrine malignancy. In China, thyroid carcinoma ranks 7th overall amongst malignant tumors threatening the health of residents [1]. Microecology theory states that the essence of life is the ecological balance of the neural, endocrine, metabolic, microbial, and immune systems [2]. Gut microbiota have been implicated in various thyroid disorders, including thyroid cancer [3], [4], [5], [6]. Obesity, a risk factor for thyroid cancer, influences the balance of intestinal bacteria and nutrients [7]. Moreover, recent studies have reported that gut microbiota may act as a regulator of immunity, and therefore may indirectly play a role in tumorigenesis and certain related immune therapies [8], [9]. It has been suggested that gut microbiota may be a promising biomarker of immune-related adverse events [10]. Therefore, it is worthwhile to explore the relationship between gut microbiota and thyroid cancer.

Recent evidence has linked thyroid cancer to gut microbiota. In 2017, a study using gas chromatography-time-of-flight mass spectrometry analyzed the serum of thyroid cancer patients with and without distant metastatic disease (n = 37 and n = 40). Thyroid cancer patients with distant metastatic disease had elevated levels of gamma-aminobutyric acid and aminooxy acetic acid, perhaps secondary to diet and/or gut microbiota [11]. Functional prediction analysis in these patients suggested a potential interaction between their intestinal flora and their thyroid glands. The following year, Feng et al examined the gut microbiota and metabolites in stool samples from thyroid cancer patients and normal controls (n = 30 and n = 35) using 16S rRNA sequencing and liquid chromatography-mass spectrometry [12]. A gut signature that included six genera was proposed as a microbial predictor to distinguish thyroid cancer patients from normal controls. The authors proposed that gut microbial dysbiosis affects lipid metabolism in thyroid cancer patients, thereby contributing to the development of cancer. In the same year, Zhang et al also observed a significant difference in gut microbiota between 20 patients with thyroid cancer and 36 normal controls, also using 16S rRNA gene sequencing [13]. Nevertheless, because these previous studies employed small sample sizes, they only offer limited information with regards to the relationship between gut microbiota and thyroid cancer. Moreover, these prior studies neglected to examine the relationship between gut microbiota and additional traits of thyroid cancer patients, such as metastatic lymphadenopathy.

Thus, we conducted a cross-sectional study examining the stool samples of thyroid cancer patients with a relatively large sample size. We employed 16S rRNA gene sequencing and multiple bioinformatics methods to comprehensively investigate the relationship between gut microbiota and thyroid cancer.

Materials and methods

Ethics statement

The protocols involved with regards to the human patients in this study were approved by the Ethics Committee of Shanghai Tenth People's Hospital (Approval no. SHSY-IEC-KY-4.0/16-18/01). All subjects were voluntarily recruited and informed of the nature of the study before sample collection. Written informed consent was obtained from all study subjects.

Study design and sample collection

This study enrolled thyroid carcinoma patients from Shanghai Tenth People's Hospital from May 2017 to March 2019 and recruited healthy people as controls from the resident community in northern Shanghai, China. Selection criteria for the thyroid carcinoma patients included the following: (i) clinically diagnosed with thyroid cancer within six months, and (ii) planning to undergo total thyroidectomy. The healthy controls were age, gender, and body mass index (BMI) matched with the thyroid cancer patients. All of the following were excluded: (i) patients not pathologically confirmed to have thyroid cancer following surgery, (ii) subjects with gastrointestinal diseases or any other severe mental or physical diseases, (iii) subjects with a history of orally ingested antibiotics, prebiotics, probiotics, or any other similar drug within the two months prior to stool sampling, and (iv) subjects with abnormal thyroid function, as indicated by the serum free triiodothyronine (fT3), serum free thyroxine (fT4), serum thyroid-stimulating hormone (TSH), serum anti-thyroid peroxidase antibody (TPOAb), and serum anti-thyroid autoantibodies (TgAb) tests.

The stool samples of thyroid patients were collected prior to surgery. Demographic information including age, gender, and BMI was obtained during subject recruitment. The following indicators were obtained from eight-hour fasting blood samples: (i) thyroid function indicators including fT3, fT4, and TSH, and (ii) thyroid antibodies including TPOAb and TgAb. The diagnosis of thyroid cancer was confirmed by two independent pathologists from Shanghai Tenth People's Hospital. Pathological TNM staging was performed for each thyroid cancer patient based on the AJCC 8th edition [14].

In total, 90 thyroid carcinoma patients (TC group) and 90 healthy controls (HC group) were included in the study. All study subjects were of Han nationality and lived in the eastern coastal provinces of China, where the typical diet includes rice, meat, vegetables, beans, and seafood, etc. The diets of all subjects were varied, and no subjects were vegetarians. The demographic information and clinical characteristics of the enrolled subjects are listed in Table 1. The TC group included 88 cases of papillary thyroid cancer (PTC) and two cases of follicular thyroid cancer (FTC). In terms of staging, the TC group included 56 cases of TCs with local lymph node metastasis (N1 group) and 34 cases of TCs without lymph node metastasis (N0 group). There were no cases with distant metastatic disease.

Table 1

Clinical characteristics of TCs and HCs in this study.

Characteristics	Exploratory cohort			Validation cohort			p-value#
	TC (n = 60)	HC (n = 60)	p-value*	TC (n = 30)	HC (n = 30)	p-value*
Demographic
Gender, male/female (female %)	20/40 (66.67%)	26/34 (56.67%)	0.260	13/17 (56.67%)	12/18 (60.00%)	0.793	0.666
Age, years, median (min–max)	41 (16–72)	44.5 (22–71)	0.964	42 (22–69)	38.5 (23–65)	0.487	0.070
BMI, kg/m2, median (min–max)	22.3 (18.8–25.2)	22.4 (19.2–24.6)	0.686	22.3 (20.2–25.2)	21.8 (18.8–23.9)	0.277	0.078
Characteristics of thyroid cancer
Pathological diagnosis, PTC/FTC	58/2	/	/	30/0	/	/	0.551a
Pathological T stage, T1-2/T3-4	57/3	/	/	29/1	/	/	0.717
Pathological N stage, N0/N1	24/36	/	/	10/20	/	/	0.539
Pathological M stage, M0/M1	60/0	/	/	30/0	/	/	/
Extrathyroidal invasion, positive/not reported	7/53	/	/	1/29	/	/	0.190
BRAF, V600E mutation/wild type/not available	16/10/34	/	/	5/4/21	/	/	0.453
Thyroid function, [range], median (min–max)
fT3, [2.8–6.3 pmol/L]	4.74 (3.29–5.99)	4.72(3.51–5.97)	0.850	4.955 (2.68–5.92)	4.585 (2.85–5.72)	0.482	0.064
fT4, [10.5–24.4 pmol/L]	12.79 (10.18–19.85)	13.09 (10.16–18.4)	0.322	13.08 (10.18–19.79)	12.81 (10.43–19.36)	0.935	0.564
TSH, [0.38–4.34 IU/L]	2.16 (1.28–3.17)	2.065 (1.12–3.26)	0.261	2.155 (1.28–3.17)	2.15 (1.5–2.84)	0.122	0.114
Thyroid antibody, [range], median (min–max)
TgAb, [<110 IU/mL]	10 (10–94)	10 (10–92)	0.726	10 (10–94)	10 (10–24)	0.246	0.073
TPOAb, [<40 IU/mL]	2 (2–35.82)	2 (2–38.76)	0.871	2 (2–25.13)	2 (2–25.18)	0.704	0.450

The statistical significance of other characteristics was tested by the Wilcoxon rank-sum test or chi-square test.

Comparisons between TCs and HCs.

comparisons between all subjects in the exploratory cohort and those in the validation cohort.

The statistical significance of diagnosis between two cohorts was tested by Fisher's exact test.

DNA extraction and 16S rRNA gene sequencing

Each stool sample was snap-frozen with liquid nitrogen following collection and was stored in a sterile container (Sarstedt, 80.734.311, Germany) within a −80 °C refrigerator (Haier, DW-86L626, China) until use [15]. Stool samples were transported to the Shanghai Majorbio Bio-pharm Technology Co., Ltd. laboratory on dry ice within three months following freezing. DNA extraction was performed using the E.Z.N.A.® Soil Kit (Omega, Bio-tek, Norcross GA, USA). The concentration of bacterial DNA was measured using the Nanodrop 2000 (Thermo Scientific, USA). DNA samples were stored at −80 °C (Haier, DW-86L626, China) until all samples were ready for sequencing. The V3-V4 region of the bacteria’s 16S rRNA was amplified by polymerase chain reaction with barcode-indexed primers (338-F: 5′-ACTCCTACGGGAGGCAGCAG-3′ and 806-R: 5′-GGACTACHVGGGTWTCTAAT-3′) using TransStart FastPfu DNA Polymerase (TransGen BioTech, AP221-02, Beijing, China) on an ABI GeneAmp®9700 polymerase chain reaction system (USA). Amplicons were then purified by gel extraction (AxyPrep DNA GelExtraction Kit, Axygen Biosciences, Union City, California, USA) and were quantified using QuantiFluor-ST (Promega, USA). The adapter sequences were added on the purified amplicons using a TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, USA), and paired-end sequencing was performed using an Illumina MiSeq instrument (Illumina, San Diego, USA).

Analysis of 16S rRNA gene sequencing data

The Majorbio Cloud Platform (https://cloud.majorbio.com/) was used to analyze the 16S rRNA gene sequencing data. Paired-end reads obtained through Miseq sequencing were first spliced according to the overlapping relationships (Flash, v1.2.11). Quality control was performed simultaneously (Fastp v0.19.6). In brief, the bases with a mass value<20 at the tail of the reads were filtered. Reads under 50 bp were also filtered, thereby removing the reads containing N-base. According to the overlapping relationships, the paired-end reads were merged into a sequence with a minimum overlap length of 10 bp and a maximum mismatch ratio of 0.2. The samples were identified based on the barcode and primers with a maximum mismatching primer of 2 permitted and no mismatching barcode accepted. Quantitative Insights Into Microbial Ecology (QIIME v1.9.1) was used to filter the sequencing reads and to construct a table featuring the 16S rRNA sequencing data. Sequences with a similarity of 97% were clustered into operational taxonomic units (OTUs) (Uparse v7.0.1090). Taxonomy was assigned using the RDP Classifier (v2.11) and Silva database (v132).

After generating an OTUs table for all subjects, a cross-sectional study was performed using a randomly generated cohort that included 60 TCs and 60 HCs. An independent cohort of 30 TCs and 30 HCs was used to validate the initial findings. The assignment of subjects was done prior to further data analysis.

Statistical analysis

A p-value of <0.05 was considered statistically significant in this study. The Kolmogorov-Smirnov test was used to test the normality (SPSS® Statistics v22). The Wilcoxon rank-sum test, independent-samples t-test, Chi-square test, and Fisher's exact test were used to compare the differences in richness, diversity, and clinical characteristics between the two groups respectively (SPSS® Statistics v22). Principal coordinates analysis (PCoA) and analysis of similarities (ANOSIM) with a permutation of 999 were used to determine whether or not inter-group differences were significantly greater than intra-group differences (Majorbio Cloud Platform). Using typing analysis, the samples were clustered into three types with the highest Calinski-Harabasz (CH) index determined by the Partitioning Around Medoids algorithm (Majorbio Cloud Platform). The Kruskal-Wallis H test was used to compare differences in taxa between the three groups (Majorbio Cloud Platform). Linear discriminant analysis (LDA) effect size (LEfSe) was used to identify the microbial taxa that significantly differed between the TCs (n = 60) and HCs (n = 60), with a threshold LDA score of > 3.5 used. LEfSe was also used to identify the microbial taxa that significantly differed between the TCs with and without metastatic lymphadenopathy, with a threshold LDA score of > 3.0 used (Majorbio Cloud Platform) [16].

Logistic regression was performed using both the forward and backward stepwise selection methods (likelihood ratio) (SPSS® Statistics v22). Lasso regression analysis was performed with the R environment (v3.6.0) using glmnet package [17]. Random forest model analysis was performed on the Majorbio Cloud Platform with default settings. Functional pathway analysis was performed using Tax4Fun and PICRUSt2, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, with the 16S rRNA gene sequencing data obtained from the TCs (n = 90) and HCs used to infer the metagenomes (n = 90) (Majorbio Cloud Platform). Differences in pathways between the TCs and HCs were identified using LEfSe with a threshold LDA score of > 2.0 used (Galaxy Platform, http://huttenhower.sph.harvard.edu/galaxy) [18].

Constructing the co-occurrence network

Based on the OTU arcsine square root transformation relative abundance data, a co-existing network was generated through weighted gene co-expression network analysis (WGCNA) with the R environment (v3.6.0) and using the WGCNA package (v1.69). A matrix table that included the top half of the OTUs (894 of 1788) with the greatest variance across all samples was utilized. Ten samples were excluded as outliers according to sample clustering. An unsigned network with nearly scale-free topology based on a soft threshold power of 6 (R2 = 0.84, slope = -1.13) was generated. The minimum number of OTUs in each module was 10. The co-occurrence network was constructed based on the inter-modular connectivity weighted values with a threshold of>0.01 utilized (Cytoscape v3.6.0). Among the modules, a grey module containing 706 OTUs was excluded, as its co-existing relationship had no significance. The heatmap correlation analysis was generated using the Spearman methods. A phylogenetic tree was also generated with IQ-TREE (v1.6.8) using the setting of maximum likelihood method.

Results

Overview of the 16S rRNA gene sequencing data

In performing 16S rRNA gene sequencing of all stool samples (n = 180), a total of 8,508,520 sequences were sorted into 1788 OTUs according to phylogenetic taxonomic levels. The good’s coverage indices for the observed OTUs in the HC and TC groups were 99.72% ± 0.078% and 99.79% ± 0.062% (Mean ± SD) respectively, confirming the thoroughness of the sampling. A broad overview of the clinical traits and taxonomic data of all subjects is provided in (Fig. 1A–C). Firmicutes, Bacteroidetes, Actinobacteria, and Proteobacteria were the four dominant phyla, with total abundances accounting for 98.95% ± 2.29% and 98.86% ± 3.78% (Mean ± SD) of the microbiomes of the HCs and TCs, respectively (Fig. 1B). The genera Prevotella_9 and Megamonas were detected in high numbers in a relatively high number of subjects, 48.33% (87/180) and 23.33% (42/180) respectively, though their abundances (with observed OTUs < 5) were quite low in the remainder of subjects (Fig. 1C), which resulted in an uneven distribution of these organisms amongst our study population.

Fig. 1

Taxonomic features of the gut microbiota of the healthy controls (n = 90) and the thyroid cancer patients (n = 90). Clinical phenotype annotations (A) for all subjects, with an overview of the 10 most abundant phyla (B) and 30 most abundant genera (C) in the stool samples identified by 16S rRNA gene sequencing. The taxa levels were determined by summing the OTUs and are presented by relative abundance.

Decreased richness and diversity of gut microbiota in TCs

To study the differences in gut microbiota between TCs and HCs, an exploratory cohort (n = 180) was randomly-generated from the two groups. This included a TC group (n = 90) and an HC group (n = 90). TC samples demonstrated reduced gut microbiota richness and diversity, as measured by the Ace index and Shannon index (p = 3.74e-05, and p = 0.0025, Fig. 2A,B) respectively. PCoA analysis on the unweighted UniFrac distance was performed to assess the overall diversity of the microbiota. ANOSIM demonstrated significant differences between the TC and HC groups (R = 0.1667, p = 0.001, Fig. 2C). These findings suggest gut microbial dysbiosis in TC patients.

Fig. 2

Bioinformatic analysis of 16S rRNA gene sequencing data in the exploratory cohort. (A-B) Box and whisker plots of the alpha diversity indices for richness (Ace index) and diversity (Shannon index) of the bacterial communities on the OTU level in thyroid carcinoma (TC) patients (n = 60) and healthy controls (HC, n = 60), respectively. (C) Principal coordinates analysis based on the unweighted UniFrac distance showed that the overall microbial diversity differed significantly between TC patients and controls (ANOSIM, R = 0.1667, p = 0.001). (D–F) Variations of stool microbiota composition in TCs. (D) The relative proportions of bacterial phyla in TC patients (n = 60) and healthy controls (n = 60). Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria are presented in different colors with other phyla grouped as “Others.” (E) The Firmicutes to Bacteroidetes (F/B) ratio did not significantly differ between the TC and HC groups. (F) The samples were clustered into three types with the highest Calinski-Harabasz (CH) indices using the Partitioning Around Medoids algorithm according to the abundance of the dominant microbial taxa. (G) The distribution of samples in the three cluster types is presented based on the unweighted Unifrac distance. The bar-plots show (H) the proportions of cases of different types in the TC and HC groups and (I) the relative abundances of the dominant phyla in the three types. (J) The microbial differences on the phylum level across the three types show that the abundance of Proteobacteria is lower in type 1 than in the other types (Kruskal-Wallis H test. ***, corrected p-value < 0.001). (K) LEfSe analysis identified the microbes whose abundances significantly differed between the TC and HC groups. The findings with regards to phylum, family, and genus are shown in the plot (LDA score > 3.5, p < 0.05). Model candidates for disease discrimination were established using multivariable binary logistic regression, random forest model analysis, and lasso regression analysis. The performance of the model candidates was assessed by area under the ROC curve for (L) the exploration data (TC: 60, HC: 60) and (M) the validation data (TC: 30, HC: 30), respectively.

Differences in microbial abundances between TC and HC subjects

On the phylum level, no significant difference in relative abundance of the top four dominant phyla was observed between the two groups (Fig. 2D). Additionally, there was no significant difference in the Firmicutes/Bacteroidetes ratio between the two groups (Fig. 2E). Thus, we further analyzed the microbiota according to the distribution of dominant organisms on unweighted UniFrac distance, and the samples were clustered into three types by the Partitioning Around Medoids algorithm (Fig. 2F–G). There were a greater number of type 1 cases in the HC group (n = 43) compared to the TC group (n = 24, Fig. 2H). Further, significant differences in Proteobacteria abundances were observed among the three types (corrected p-value = 9.26e-04, Fig. 2I and J). These findings suggest that TC patients have a relatively higher abundance of Proteobacteria in the gut.

In the exploratory cohort, significant differences in microbial abundances were observed between the two groups. Three genera were over-represented in TCs, and seven genera were over-represented in HCs on the genus level (Fig. 2K). In addition, the LEfSe results validated that Proteobacteria levels were significantly higher in the TC group (LDA score = 4.2, p = 0.009).

A gut microbiota-based signature can predict thyroid cancer status

Based on the initial LEfSe findings, we next assessed the value of gut microbiota as biomarkers using forward stepwise logistic regression, random forest model analysis, and lasso regression. The performance of each of the candidates was assessed using a receiver operating characteristic (ROC) curve, firstly on the exploratory cohort data. Using the stepwise logistic regression model and the random forest model as signature-predictors yielded areas under the curves (AUC) of 0.797 and 0.804, respectively. Lasso regression suggested that all 10 genera could be employed as combined predictors, which yielded an AUC of 0.809 (Fig. 2L).

The above findings were then validated using the independent validation cohort (TC vs. HC, n = 30 and n = 30, Fig. 2M). The lasso model performed the best, with an AUC of 0.746. Therefore, this signature combination of 10 genera has the highest potential in terms of TC diagnostic value. This signature includes increased levels of g__Bacteroides, g__Lachnoclostridium, and g__norank_f__Lachnospiraceae, and reduced levels of g__Prevotella_9, g__Collinsella, g__Faecalibacterium, g__Dorea, g__Ruminococcaceae_UCG-014, g__Ruminococcaceae_UCG-002, and g__Subdoligranulum.

Microbial functional dysbiosis in TCs

To study the functional and metabolic-related changes of the gut microbiota in TC patients, we conducted functional prediction analysis. A total of 273 KEGG categories were obtained by Tax4Fun, and 337 KEGG categories were obtained by PICRUSt2 (Fig. 3A). LEfSe analysis identified 30 Tax4Fun KEGG categories that significantly differed between the two groups and 52 PICRUSt2 KEGG categories that significantly differed between the two groups (Fig. 3B). The fewer KEGG categories obtained using Tax4Fun and their respective p-values suggest that Tax4Fun's estimates are more conservative. Twenty one common KEGG categories that significantly differed between the two groups were identified by both Tax4Fun and PICRUSt2 (Fig. 3C and D). Among these 21 categories, the majority were related to metabolism (13/21, 61.9%). Of note, five categories related to the processing of genetic information, including “Aminoacyl-tRNA biosynthesis,” “Homologous recombination,” “Mismatch repair,” “DNA replication,” and “Nucleotide excision repair” were significantly increased in the HCs, which suggests deficient genetic information processing in TCs.

Fig. 3

Functional analysis prediction. (A) The Venn plot shows the predicted KEGG function counts for the TCs (n = 90) and HCs (n = 90) based on Tax4Fun and PICRUSt2. (B) 21 KEGG terms significantly differed between the TCs and HCs (LEfSe analysis, LDA score > 2, p < 0.05). The bubble plot shows the 21 KEGG categories at the intersection of the (C) Tax4Fun and (D) PICRUSt2 findings, with significant differences between TCs and HCs. The categories increased in the TC group are marked with a triangle, and the categories increased in the HC group are marked with a circle.

Gut microbiota and local lymph node metastasis in TCs

Lymph node metastasis is common in thyroid cancer. In this study, 62.22% (56/90) of TC patients had metastatic lymphadenopathy at presentation. Therefore, we subdivided the TC patients into N0 (n = 34) and N1 (n = 56) subgroups. No significant difference in gut microbiota richness or diversity was observed between the two subgroups (Ace index: p = 0.111; Shannon index: p = 0.632), and no significant difference in overall distribution on the unweighted Unifrac-based PCoA was observed either (ANOSIM, R = −0.0019, p = 0.507).

We next analyzed the relative abundances of microbiota in the N1 and N0 subgroups. The relative abundances of the 10 TC-related genera showed no significant difference across the two groups (p-value > 0.05). LEfSe analysis identified five genera with increased abundances in the N1 group, and five genera with increased abundances in the N0 group (Fig. 4A). We applied back-forward stepwise logistic regression, random forest model analysis, and lasso regression to establish a microbial model that could predict the presence of local metastatic lymphadenopathy. An eight-genera model established using lasso regression performed the best with an AUC of 0.797. Meanwhile, a four-genera model (g__Hungatella, g__Alistipes, g__Fusobacterium, and g__Phascolarctobacterium) generated using back-forward stepwise logistic regression, performed nearly as well, with an AUC of 0.778 (Fig. 4B). Because of the increased simplicity of the second model, which requires fewer genera, the four-genera model may be most suitable to characterize the gut microbiota changes that occur in thyroid cancer patients with metastatic lymphadenopathy.

Fig. 4

The relationship between gut microbiota and lymph node metastasis in TCs. (A) LEfSe analysis identified the genera whose abundances significantly differed between the N0 and N1 subgroups (LDA score > 3, p < 0.05). (B) Three model candidates were established using multivariable binary stepwise logistic regression (backward), random forest model analysis, and lasso regression analysis, respectively. The performance of the above three model candidates was assessed by area under the ROC curve. (C) The Venn plot shows the predicted KEGG functions of the microbiota in the N1 and N0 subgroups respectively based on Tax4Fun and PICRUSt2. The bubble plots show the KEGG categories that significantly differed between the N1 and N0 groups according to (D) Tax4Fun and (E) PICRUSt2 (LEfSe analysis, LDA score > 2, p < 0.05). The categories increased in the N1 subgroup are marked with a triangle, and the categories increased in the N0 subgroup are marked with a circle.

Using functional prediction analysis, a total of 272 KEGG categories were obtained by Tax4Fun, and 313 KEGG categories were obtained by PICRUSt2 (Fig. 4C). LEfSe analysis identified only 18 KEGG categories that significantly differed between the N0 and N1 subgroups. (Fig. 4D–E).

The OTU co-existing network identified a potential regulator

Gut bacteria constitute a complex ecosystem where certain species are not only impacted by the host but also by other bacteria within the community. A total of 188 co-expressed OTUs were clustered by WGCNA into eight modules, each of which was labeled by color and constituted a distinct network (Fig. 5A). After mapping the OTUs to corresponding phyla, the phylogenetically related OTUs clustered into the same modules preferentially, although each module also contained OTUs from different taxa. Firmicutes was the dominant phyla within the co-existing network (Fig. 5B).

Fig. 5

OTU co-existing network and module-trait associations. (A) OTU co-existing network where OTUs (nodes) are colored according to the WGCNA clusters. (B) OTU co-occurrence network where OTUs (nodes) are colored according to the phyla to which they are classified. The thickness of the line (edges) represents the weighted value of node connectivity. (C) Module–trait associations are presented by heatmap. Each cell in the matrix contains the Spearman correlation coefficient between one OTU module and a clinical trait as well as the corresponding p-value. (D) Phylogenetic trees based on the maximum likelihood method showing the phylogenetic relationships between the OTUs in the pink module. The confidence of the branch generated by the bootstrap test is illustrated. (E) The membership score of each OTU in the pink module is shown in the dot-plot.

The associations between modules and clinical characteristics are presented by heatmap (Fig. 5C). A negative correlation between the pink module and disease status (r = -0.18, p = 0.02) was observed. In the pink module, OTU457 and OTU1399 display a relatively high numbers of reads (Fig. 5D). Five OTUs (OTU1022, OTU1245, OTU511, OTU457, and OTU490) potentially play a dominant role, with membership scores > 0.75 (Fig. 5E). Of note, because of its high levels and dominant role, the OTU457 (g__Ruminococcaceae_UCG-002) may act as a regulator within the pink-module community. In addition, these findings also serve as validation of the significant difference in g__Ruminococcaceae_UCG-002 levels between TCs and HCs.

Discussion

Recent evidence suggests that the gut microbiome, a sophisticated and metabolically active community, likely influences thyroid homeostasis in the host [19]. Herein, we characterized the microbiome of TC patients by examining the stool samples of a relatively large TC cohort by 16S rRNA gene sequencing. Our data demonstrate that the microbiomes of TC patients are reduced in richness and diversity, and there are significant changes in the relative abundances of 10 genera. Based on inferred functional analysis of these genera, our results suggest that the microbiomes of TC patients lack certain genetic information processing capabilities, such as base excision repair and nucleotide excision repair.

The microbiomes of patients with TC exhibited decreased microbial richness and diversity, compared to healthy controls. Reduced microbial diversity has been documented in a variety of diseases and is considered to be one of the major characteristics of gut microbial dysbiosis. However, upon reviewing the two previous studies linking gut microbiota and thyroid cancer [12], [13], some of our findings are somewhat contradictory. Zhang et al found that the gut microbiomes of patients with thyroid cancer are significantly increased in richness compared to controls, as measured by the Ace index [13]. Feng et al also found an increased richness and diversity of gut microbiota in TCs, as indicated by the Chao and Shannon indices. Our somewhat contradictory results may stem from various factors including demographics of the controls, tumor TNM status, diet, etc. In the previous studies, only a few of the TC patients presented with metastatic lymphadenopathy, however, 62.22% of our TC population presented with metastatic lymph nodes. Though we found no significant relationship between gut microbiota richness and diversity in thyroid cancer patients with and without metastatic lymphadenopathy, our findings do suggest a gut dysbiosis exists in TC patients. In addition, it has been proposed that the microbiome can be readily reshaped by dietary exposures [20]. The noted prior studies were conducted in northeastern China, while our study was conducted in Southeastern China. Therefore, the diet habits of the subjects in our study likely differ significantly from those of subjects in the previous studies. We speculate that regional and environmental variables may influence the gut microbiome and therefore research in this area.

In addition to changes in richness and diversity, we found elevated levels of Proteobacteria in the majority of TCs through typing analysis. Proteobacteria is considered a signature of microbial dysbiosis [21]. The overgrowth of these bacteria may be related to obesity and a high-fat diet. Although the TCs and HCs in our study showed no significant difference in BMI index, it is possible that changes in intestinal microbiota related to diet may occur prior to increases in BMI. The elevation in Proteobacteria levels may signify gut dysbiosis related to the underlying pathogenesis of the thyroid cancer, especially because obesity is one of the risk factors for thyroid cancer.

In our study, the intestinal microbial changes that occurred in thyroid carcinoma patients were best characterized by 10 genera. The decreased relative abundance of Prevotella in TCs is in accordance with prior studies [12], [13]. It has been proposed that Prevotella is one of the dominant microorganisms in the human gut, especially in hosts who follow a plant-based diet [22]. The low levels of Prevotella in TC patients may therefore suggest a relationship between the ingestion of animal-derived foods and thyroid carcinoma. However, some components of the model we derived to predict TC differed from previous studies, which showed greater heterogeneity of the intestinal flora of TC patients. Research on the role of intestinal flora in thyroid cancer is still rare. A meta-analysis based on all of the prior research in this area may be helpful in explaining the heterogeneity of gut microbiota in TC patients.

TC patients with metastatic lymphadenopathy demonstrated a microbial signature that consisted of increased levels of g__Fusobacterium and g__Alistipes, as well as decreased levels of g__Hungatella and g__Phascolarctobacterium. None of these four genera were a component of the microbial signature that we used to distinguish TCs from HCs, but perhaps they play a role in promoting metastasis. For example, Fusobacterium has been reported to be persistently associated with distant metastases in human colorectal cancer [23], which suggests its potential as a noninvasive biomarker of metastatic disease. Moreover, the genus Alistipes has also been documented as a microbial regulator of pathogenicity. [24]. Some studies indicate Alistipes is pathogenic in colorectal cancer and is related to mental depression [25], [26], [27]. Further, Fusobacterium and Hungatella have been implicated as methylation-regulating bacteria in colorectal cancer, which may lead to upregulation of DNA methyltransferase based on mechanistic validation studies using cell lines and animal models [28]. Aberrant DNA methylation is one of the features of thyroid carcinoma and is related to its development [29]. Therefore, we infer that gut microbiota changes may influence circulating DNA methyltransferase levels, thereby promoting the development and spread of thyroid carcinoma.

We deduce that the loss of short-chain fatty acid (SCFA)-producing bacteria may promote the development of thyroid cancer. These bacteria include Faecalibacterium, Ruminococcaceae_UCG-002, and Phascolarctobacterium. Faecalibacterium, one of the dominant genera in our population, is known as a commensal biological indicator of human health, as well as an SCFA producing genus [30], [31]. SCFA-producing genera are typically found in abundance within the human gut, benefiting the host by participating in the fermentation of carbohydrates into SCFAs, including acetate and propionate. It has been reported that SCFAs may induce apoptosis or cell injury in colon cancer cells [32]. Thus, we believe a low abundance of Faecalibacterium may result in gut microbial dysbiosis either prior to or following the development of thyroid cancer. Further, accumulating evidence suggests that the abundance of Phascolarctobacterium drops in many diseases. For example, Sung et al. suggested that reduced Phascolarctobacterium levels can predict clinical recurrence of hepatic encephalopathy in patients with cirrhosis of the liver [33]. TCs in our study had relatively low levels of Ruminococcaceae_UCG-002, which was identified as a potential regulator (OTU457). This may be a novel finding regarding the relationship between thyroid cancer and gut microbiota. Future studies with a larger sample size would be useful to validate this hypothesis.

Our results suggest that the microbiomes of TCs are deficient in their genetic information processing capabilities compared to HCs. A previous study proposed that the inhibition of aminoacyl-tRNA biosynthesis is an effective antimicrobial strategy [34], which may explain the decreased richness and diversity we saw in the gut microbiomes of TCs. A loss of the homologous recombination function has been documented as one of the key factors of DNA repair failure [35]. Moreover, the biological functions of mismatch repair, DNA replication, and nucleotide excision repair play important roles in genomic stability and tumorigenesis [36], [37]. Genetic events are common in the development of thyroid carcinoma including BRAF mutations, RAS mutations, etc [38]. Therefore, based on the predicted functions of the bacteria, changes in the microbiome may be involved in thyroid cancer pathogenesis.

Though the concept of studying the relationship between gut microbiota and thyroid cancer was not very novel, our study has the following strengths. Our sample size was relatively larger than previous studies that examined the relationship between gut microbiota and thyroid cancer [12], [13]. As such, an independent cohort was able to be reserved for validation in our study. Using the independent cohort to confirm the performance of the established microbial models increases the credibility of the results and boosts our confidence in the models' abilities to distinguish the respective groups. Further, our study was conducted in a relatively comprehensive manner. In addition to using multiple methods to generate microbial models, we conducted a subgroup analysis of thyroid cancer patients with and without metastatic lymphadenopathy and used the WGCNA method to explore the potential co-existing relationship among intestinal flora. Given the results of this study, there are potential clinical applications. Stool is a promising biological sample for non-invasive diagnosis. Further, the regulation of intestinal microecology with probiotics and prebiotics may become a valuable approach to the management of thyroid cancer.

However, there are a few limitations to this single-center cross-sectional study. Because of the study design, we were unable to observe changes in microbiota that occur over time or after the development of cancer. A longitudinal study would be useful to demonstrate that changes in an individual's microbiota occur following the development of disease. In addition, some of our findings were contradictory to the results of two previous studies [12], [13]. A multi-center study with a larger number of subjects would be useful to determine what trends occur across a more varied patient population. Additional factors could also be considered when designing future research such as diet, the Bristol stool scale, etc. Further, when using 16S rRNA gene sequencing to study gut microbiota, neither the microbial information identified at the species level nor the predictions of microbial function are comprehensive enough. Accordingly, adoption of the shotgun metagenomics method could help to overcome these limitations. Finally, our findings should be verified using animal models. The relationship between certain microbes and thyroid cancer could be further explored through fecal microbiota transplantation in thyroid cancer-bearing mice.

Conclusion

Our results offer greater insight into the host-microbiota relationship and may serve as an important resource for researchers who plan to perform studies in this field in the future. Our findings help to better characterize the intestinal microbial dysbiosis that occurs in TC patients. Functional prediction based on the microbial dysbiosis suggests that faulty processing of genetic information may link certain gut microbiota to the pathogenesis of TC.

Declaration

Availability of data and materials: The datasets generated and analyzed in this study are not publicly available, as they are also being used in an ongoing study, but are available from the corresponding author on reasonable request.

Funding

This study was sponsored by National Natural Science Foundation of China [grant number: 82071964], Shanghai Shenkang Three-year Action Project [grant number: SHDC2020CR2054B], Program of Shanghai Academic Research Leader [grant number: 18XD1403000], Shanghai Leading Talent Program sponsored by Shanghai Human Resources and Social Security Bureau [grant number: 03.05.19005], and Shanghai Science and Technology Commission Support Plan [grant number: 18441903500].

Authors' contributions

X.Y., W.J., Z.L., and D.L. conceived and designed the project. Z.L. obtained funding. Y.S., Q.L., S.L., and S.Q. performed clinical diagnosis. X.Y., W.J., Y.S., T.Q., J.T., S.L., C.D., and Q.S. collected samples. X.Y., Y.S., W.J., R.O.K., T.Q., and J.T. contributed to data collection. X.Y., W.J., R.O.K., Y.S., D.L., and Z.L. analyzed and interpreted data and performed result visualization. X.Y. and W.J. wrote the manuscript. R.O.K., D.L., and Z.L. reviewed and edited the manuscript. All authors approved the final version of the manuscript. X.Y., W.J., and R.O.K. contributed equally to this work and share co-first authorship.

Ethics statement

The protocols involved in human patients that used in this study were approved by the Ethics Committee of Shanghai Tenth People’s Hospital (Approval no. SHSY-IEC-KY-4.0/16–18/01). All subjects were voluntarily recruited and informed of the nature of the study before sample collection. Written informed consents were provided by all subjects.

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.