This study recruited 50 TC patients (22 males, 28 females, ages 16-69 years, median age 41) from the Department of Nuclear Medicine, Shanghai Tenth People’s Hospital between May 2017 and February 2019. The control group consisted of 58 healthy volunteers (22 males, 36 females, ages 22-71 years, median age 44) who visited this hospital for their routine physical examination. The diagnosis of TC was based on the following criteria: (i) clinically diagnosed with differentiated thyroid cancer within six months (13), (ii) age between 18 and 75 years, female (nonpregnant) or male without history of other malignancies, and (iii) receiving surgical treatment. The healthy volunteers were age-, sex- and BMI-matched with the TC patients. The following patients were also excluded: (i) patients confirmed not to have thyroid cancer on postoperative pathology, (ii) subjects suffering from a gastrointestinal disease or any other severe physical or mental disease, and (iii) subjects with a history of probiotics, orally ingested antibiotics, or any other similar drug within two months before fecal sampling. Fecal samples of the TC patients were collected prior to surgery. The following indices were obtained from the eight-hour fasting blood samples: (i) thyroid function indices, including serum free triiodothyronine (fT3), serum free thyroxine (fT4) and serum thyroid stimulating hormone (TSH); (ii) thyroid antibodies, including serum anti-thyroid peroxidase antibody (TPOAb) and serum anti-thyroid autoantibodies (TgAb); and (iii) serum biochemical indicators, including alanine transaminase (ALT), aspartate transaminase (AST), triglycerides (TGs), serum total cholesterol (TCHO), high-density lipoprotein cholesterol (HDL), and low-density lipoprotein cholesterol (LDL). The diagnosis of TC was determined by two independent pathologists in Shanghai Tenth People’s Hospital. All research subjects were of Han nationality and resided in the eastern and southeastern provinces of China. This study was conducted in accordance with the principles of the World Medical Association and the Declaration of Helsinki and was approved by the Shanghai Tenth People’s Hospital Ethics Committee (No. SHSY-IEC-KY-4.0/16-18/01). Written informed consent was obtained from the subjects.

Fecal samples were collected the morning after an overnight fast (more than 8 hours). Fecal samples were snap-frozen with liquid nitrogen following collection and stored at -80°C (Haier, DW-86L626, China) (14). Microbial genomic DNA was extracted using the E.Z.N.A. soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA) in accordance with the manufacturer’s protocol. A Nanodrop 2000 (Thermo Scientific, USA) was used to measure the DNA concentration. Extracted DNA was stored at -80°C until experiments were conducted. The primers (338-F: 5’-ACTCCTACGGGAGGCAGCAG-3’) and (806-R: 5’-GGACTACHVGGGTWTCTAAT-3’) targeting the V3-V4 hypervariable region of the 16S rRNA were used to amplify the extracted DNA samples by the polymerase chain reaction system (ABI GeneAmp 9700, USA). The amplicons were subsequently purified by AxyPrep DNA Gel (Axygen, CA, USA) and quantified by QuantiFluor-ST (Promega, USA). Sequencing was performed on an Illumina MiSeq platform (Illumina, San Diego, CA) to combine the purified amplicons at equimolar concentrations and perform end pair sequencing analysis.

Fecal sample (50 mg) was weighted to an EP tube, and 400 µL extract solution (methanol: water = 4: 1) was added. Then they were homogenized by a Wonbio-96c high-throughput tissue crusher (Shanghai Wanbo Biotechnology Co., Ltd.) for 6 min at 50 Hz. Then, they were vortexed for half a minute and ultrasonicated for 30 min at 40 kHz (incubated in ice water). The mixture was incubated for 30 minutes at -20°C to precipitate proteins and then centrifuged at 13,000 g for 15 min at 4°C. The supernatant of the samples was transferred to fresh EP tubes and analyzed by ultra-performance liquid chromatography (UPLC) system coupled to a quadrupole-time-of-flight mass spectrometer (Triple TOFTM5600+, AB Sciex, USA) equipped with a electrospray ionization (ESI) source operating in both positive mode and negative ion mode. Additionally, 10 μL of supernatant from each sample was combined and mixed to serve as the quality control (QC) sample, which was injected randomly (every 8 samples) to screen the consistency of the analysis.

Nontargeted liquid chromatography-mass spectrometry (LC–MS/MS) analyses were performed for peak detection and comparison by importing the raw data into the Progenesis QI 2.3 software program (Nonlinear Dynamics, Waters, USA), and a data matrix composed of retention time (RT), peak intensity and mass-to-charge ratio (m/z) values was obtained. Additionally, qualitative and quantitative information for specific metabolites was acquired. Preprocessing was necessary before the multivariate statistical analysis, and the methods were as follows: (i) data were removed if missing values were more than 50% in the sample; (ii) based on one-half of the minimum metabolite value, the missing values were filled and the total peaks were normalized; and (iii) QC samples were eliminated if the relative standard deviation (RSD) was more than 30%. Finally, logarithmic transformation was performed for preprocessed datasets prior to further analysis. Principal component analysis (PCA) was used to analyze the cation and anion set data for QC analysis. Mass spectra of these metabolic features were identified by using the accurate mass, MS/MS fragments spectra and isotope ratio difference with searching in reliable biochemical databases as Human Metabolome Database (HMDB) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Specifically, the mass tolerance between the measured m/z values and the exact mass of the target components was ± 5 ppm. For MS/MS-confirmed metabolites, only metabolites with MS/MS fragments scores above 30 were considered reliably identified. Otherwise, the metabolites had only tentative assignments.

The analysis of 16S rRNA gene sequencing data was performed in the Majorbio Cloud Platform (https://cloud.majorbio.com/). Paired-end reads were first assembled using FLASH (v1.2.11) (15). Then, the sequences were optimized by filtering and quality control using Fastp (v0.19.6) (16). We then used UPARSE (v 7.0.1090) to cluster sequences by a similarity of 97% into operational taxonomic units (OTUs) (17).

The annotation of each taxonomy was performed by the RDP Classifier (v2.11) and the Silva database (v132). To decrease the effect of spurious sequences, sequences of OTUs <0.005% were removed (18).

The alpha diversity was analyzed to determine the richness of species at the OTU level based on the Ace index. The beta diversity was analyzed using the unweighted UniFrac algorithm to compare the composition of bacterial communities among samples. On the Majorbio Cloud Platform, the analysis of principal coordinates and similarities were performed to determine whether the intergroup differences were significantly greater than the intragroup differences or not. The Partitioning Around Medoids algorithm was used to cluster individuals into two enterotypes. The differences at the phylum or genus levels between the two groups were assessed by the Kruskal–Wallis H test. Linear discriminant analysis effect size (LEfSe) analysis was applied to determine the significantly different genera between the TC and HC groups (19).

Novel microbial and metabolic predictors were explored using logistic regression by forward stepwise (likelihood ratio) methods (SPSS Statistics v22), lasso regression and random forest model analysis (R environment v.4.1.0) (19, 20). The diagnostic ability of the selected substances was characterized by receiver operator characteristic (ROC) curve analysis with comparisons of the areas under the curve (AUCs). Differences between two AUCs were evaluated by the Z test. The final predictive models were validated in an additional cohort, consisting of 20 patients with thyroid cancer and 20 healthy controls via measuring microbiota and metabolites. Canonical correspondence analysis (CCA) was performed to conduct correlation analysis of clinical and metabolic factors among samples. The Spearman correlation between the two groups was determined for each factor.

Differential metabolites between the two groups were screened by orthogonal partial least squares discrimination analysis (OPLS-DA) and Student t-test (data conformed to a normal distribution and satisfied the homogeneity of variance) in the SIMICA software program (v 16.0.2) (21). The analysis of functional pathways was performed in the Majorbio cloud and MetaboAnalyst platform (v5.0) (22). The network between changed genera and lipid metabolism was visualized by Cytoscape (v3.6) software (23). A difference of statistical significance was considered only when the value of p was less than 0.05.

Our research enrolled 50 patients with TC and 58 healthy volunteers. There were no significant differences in age and sex between the two groups. The baseline characteristics of the TC patients and HCs are summarized in Table 1. Then we analyzed fecal samples from TC and HC groups using 16S rRNA gene sequencing and untargeted metabonomic analyses of their gut microbiota and metabolites.

After quality filtering, a total of 5,666,615 sequences from 108 samples were distributed into 1695 OTUs. According to the rank-abundance curves at the OTU level, the Good’s coverage indices of both TC and HC groups were higher than 99.5%, suggesting a satisfactory sequencing depth of the intestinal flora (Figure 1A). The microbial diversity characteristics and their relationship with clinical parameters are shown in Figure 1B. The top four phyla were Firmicutes, Bacteroidetes, Actinobacteria, and Proteobacteria (Figure 1B). At the genus level, Bacteroides was the most abundant genus, with an average relative abundance of 10.5%. Meanwhile, Prevotella_9 and Megamonas showed bias distribution in our samples. We found that the enrichment of Prevotella_9 (sequences in OTUs > 100) was lower in TC patients (n = 10, 20.0%) than in HCs (n = 26, 44.8%; χ² = 7.448, p < 0.05). Conversely, the enrichment of Megamonas was higher in TCs (n = 9, 18.0%) than in HC patients (n = 3, 5.2%; χ² = 4.474, p < 0.05) (Figure 1B). The α-diversity of gut microbes in TC patients was significantly decreased compared with that in HCs (p < 0.0001), which suggests that the richness of the gut microbiota is reduced in TCs (Figure 1C). The β-diversity was shown on the overall microbial composition of the microbiota based on the results of principal coordinate analysis (PCoA), which revealed obvious differences between the TC patients and HCs (R = 0.1270, p = 0.001, Figure 1D).

According to the PCoA analysis in the previous section, although the difference between the samples was greater than that within the group, the R value in the group analysis was very small. By partitioning around medoids algorithm, individuals can be robustly clustered into two enterotypes (type 1, n=64; type 2, n=44) (Supplementary Figure 1A), and their distribution of sample differences is shown in Figure 2A. Type 1 comprised 74% of TC patients, but the ratio of the two types were almost equal in HCs (Figure 2B). Moreover, at the phylum level, significant differences were only observed in the abundance of Proteobacteria (higher in type 1, p < 0.001, Figure 2C). At the genus level, Bacteroides, Ruminococcus_2, Ruminococcaceae_UCG-013 were enriched in type 1, and Prevotella_9 and Ruminococcaceae_UCG-002 had greater enrichment in type 2 (p < 0.05, Figures 2D, E).

Next, we performed preliminary discriminant analysis on bacterial community structure at different taxonomic levels to identify the differentially abundant taxa in the TC group. At the phylum level, the abundance of Proteobacteria was higher in the TC group than in HCs (p < 0.05, Figures 2F, G). Then, we calculated the ratio of Firmicutes/Bacteroidetes, and there were no differences in the ratio between the two groups (Supplementary Figure 1B). Meanwhile, we did not find significant correlations between the ratio and age or BMI in either group (Supplementary Figures 1C, D). At the genus level, LEfSe analysis demonstrated that a total of 25 changed genera could be regarded as biomarkers to distinguish TC patients from HCs (seven dominant genera in TC; 18 in HC, Figure 2H). In addition, the LEfSe results demonstrated that the g_Christensenellaceae_R-7_group (LDA score = 3.3, p < 0.001) and g_[Eubacterium]_coprostanoligenes_group (LDA score = 3.5, p < 0.001) levels were significantly decreased in the TC group (Figure 2H).

According to the 25 altered genera identified by LEfSe analysis, we subsequently explored novel microbial predictors for noninvasive screening and diagnosis of TC. Lasso regression validated that 23 genera could be used as combined markers. Three genera were indicated in forward stepwise logistic regression, and seven genera were indicated in the random forest model (Supplementary Table 1, Supplementary Figures 2A–C). As we preferred to include relatively few floras to set up the prediction model, we did not include the results of the lasso regression in the following research. Then, we compared the composition of genera between the other two analyses (Figure 3A). Based on the ROC curve, the AUCs (95% CI) of the two models were 0.828 (0.75-0.91, random forest), 0.813 (0.73-0.90, logistic regression)(Figure 3B). However, there were no significant differences among the pairwise comparisons of the AUCs of each model (Z test, p > 0.05; Supplementary Table 2).

To determine the influence of clinical factors on our sample distribution, CCA was performed. After filtering with variance inflation factor analysis, the CCA results shown in Figure 3C demonstrated that TGs, TCHO and LDL had the strongest contribution to the sample distribution (r² > 0.176, p < 0.01, Supplementary Table 3). Furthermore, Spearman correlation analysis was performed between 25 differential genera (by LEfSe analysis) and 17 clinical characteristics. It was obvious that the relationships between serum TG levels and differential genera were completely different from other clinical indicators at the genus level (Figure 3D). TG was significantly positively correlated with Prevotella_9 (r = 0.33, p < 0.001) and Ruminococcaceae_UCG-010 (r = 0.31, p < 0.01) and significantly negatively correlated with Eggerthella Igolzaea (r = -0.40, p < 0.001).

All intestinal bacteria constitute a complex ecosystem, where microorganisms are influenced not only by the physiological activities of the host but also by other genera in the microecological environment. Therefore, Spearman correlation was analyzed on the 25 genera that were altered between TCs and HCs. We found 10 clusters of flora with similar expression patterns in our samples by clustering. Interestingly, Christensenellaceae_R-7_group, which was found in one of the clusters, had a strong correlation with [Eubacterium]_coprostanoligenes_group (r > 0.60, p < 0.001) and five other genera within the same cluster (Figure 3E), which indicated a strong coexistence among them.

To evaluate whether the profiles of metabolites were related to TC, we performed metabolic profiling in all fecal samples. In our study, 4385 cationic ions (363 metabolites identified) and 4453 anionic ions (161 metabolites identified) were obtained (Supplementary Material). As shown in Figures 4A, B, the differences in classification between TC patients and HCs were not obvious; there were few differences in composition between the two groups. Since evidence from this study indicated that the intestinal flora might affect the metabolism of fatty acids, we subsequently focused on fatty acid metabolites.

HMDB annotation results identified 56 fatty acyl groups, 39 sterols and their derivatives (Figure 4C); KEGG annotation results identified three fatty acid compounds (Figure 4D). OPLS-DA analysis (Figures 4E, F, Supplementary Figures 3A, B) identified a total of 49 differential positive compounds involved in human metabolism between the TC and HC groups, and 21 negative compounds were identified (variable importance in projection (VIP) > 1.0, p < 0.05, Tables 2, 3). The dominant differential metabolites of positive and negative cations that were indicated by the HMDB database are shown in Figures 4G, H.

On the basis of the 70 differential metabolites between the two groups that were identified in the previous section, three methods were employed to construct metabolite models for TC prediction (Supplementary Table 4 and Supplementary Figures 3C–E). The metabolites included in the three models and the diagnostic effectiveness of the three models are shown in Figures 4K, I. These results indicate that the stepwise logistic regression classifier shows great predictive ability with relatively few metabolites (classifier II).

Next, we entered all 25 genera and 70 metabolites with apparent distinction between the TC and HC groups into forward stepwise logistic regression (classifier III, combination of genus and metabolites, Supplementary Table 5). The distinguishing ability of the models among genera (classifier I, the composition of seven differential genera generated by random forest model analysis), metabolites (classifier II) and their combination (classifier III, Figures 4L, M) was compared. Adlercreutzia was the only same genus between classifiers I and III, and the common metabolites between classifiers II and III were isopersin, erinapyrone B, 3-sulfodeoxycholic acid, capsaicin, salicyluric acid and uric acid. Moreover, as shown in Figure 4J and Supplementary Table 6, the metabolite model (AUC = 0.940, 95% CI: 0.90-0.98) was significantly better than the microbial model (AUC = 0.828, 95% CI: 0.75-0.91, Z = 2.504, p = 0.012). Additionally, the combined model (AUC = 0.967, 95% CI: 0.93-1.00) was markedly better than the microbial model (Z = 3.198, p = 0.001). However, there was no greater distinction between the combined model and the metabolite model (Z = 1.029, p = 0.304). Subsequently, the above results were validated by an additional validation cohort (TC vs HC, n = 20 and n = 20, Supplementary Figure 4), indicating the accuracy of our predictive models.

Next, we aimed to explore the functional influence of these metabolites and characterize which metabolic pathways they might be related to. There were seven enriched metabolic pathways in the MetaboAnalyst platform with no statistical significance (p > 0.05 of all, Supplementary Table 7 and Figure 5A), and 66 networks in the Major Cloud platform were determined to be significantly different by the KEGG pathway enrichment analysis (p < 0.05, Supplementary Table 8 and Figure 5B). The pathways related to steroid synthesis enriched by the metabolite cholesterol were highly significant in both platforms. For example, cholesterol and 7-dehydrocholesterol were enriched in the steroid biosynthesis [M00101] pathway in the MetaboAnalyst platform (p = 0.184, corrected p > 0.05), and cholesterol and DG (16:0/18:0/0:0) showed strong enrichment in the pathway of fat digestion and absorption by the Major Cloud platform (corrected p = 0.028). Notably, the abundance of metabolic compounds (cholesterol, 7-dehydrocholesterol, DG (16:0/18:0/0:0) enriched in the above pathways was decreased in the TC group. These results suggest that functional pathways, such as steroid biosynthesis and lipid digestion, might be inhibited, which would further lead to the dysfunction in the TC group.

Based on the above results, we subsequently analyzed the potential association between genera and metabolic compounds, especially lipid metabolites. There were a total of 10 lipid compounds among 70 differential metabolites that were identified using the KEGG pathway analysis of the Major Cloud platform (Supplementary Figures 5A–J). Aminocaproic acid, ricinoleic acid, linoelaidic acid, and DG (16:0/18:0/0:0) were more abundant in TC patients; however, 27-hydroxycholesterol (27HC), cholesterol, 7-dehydrocholesterol, deoxycholic acid glycine, conjugate, phytosphingosine, and solasodine were more abundant in HCs. After the filtering with variance inflation factor analysis, the CCA results (Figure 5C) demonstrated that more than 10 lipid compounds significantly impacted the distribution of flora in the samples, and 27HC had the most robust influence (r² = 0.3838, p = 0.001), followed by cholesterol (r² = 0.2637, p = 0.001) and 7-dehydrocholesterol (r² = 0.2497, p = 0.001, Supplementary Table 9). Significant relationships between the TC and HC groups were uncovered by Spearman correlation analysis of 28 differential bacteria and 10 differential lipid compounds (Figure 5D). 27HC positively correlated with g_Christensenellaceae_R-7_group (r = 0.529, p < 0.001), g_Eubacterium_coprostanoligenes_group (r = 0.504, p < 0.001), g_Ruminococcaceae_UCG-002 (r = 0.544, p < 0.001), and g_Ruminococcaceae_UCG-005 (r = 0.557, p < 0.001). Then, based on the above correlations, the intestinal fecal genera and lipid compounds were further analyzed by the Spearman correlation network (Figure 5E). The core metabolites in both groups were 27HC and cholesterol, and the key bacteria were g_Christensenellaceae_R-7_group, g_Eubacterium_coprostanoligenes_group, g_Rumenococcus UCG-002, g_Rumenococcus UCG-005, and g_Rumenococcus UCG-014.