Integrative analysis of gut microbiome and host transcriptome reveal novel molecular signatures in Hashimoto's thyroiditis

Abstract

Background

Hashimoto's thyroiditis (HT) is an autoimmune disorder with unclear molecular mechanisms. While current diagnosis is well-established, understanding of the gut-thyroid axis in HT remains limited. This study aimed to uncover novel molecular signatures in HT by integrating gut metagenome and host transcriptome data (miRNA/mRNA), potentially elucidating disease pathogenesis and identifying new therapeutic targets.

Methods

We recruited 31 early HT patients and 30 healthy controls in a two-stage study (discovery and validation). Blood and fecal samples underwent RNA and metagenomic sequencing, respectively. Integrative analysis included differential expression, weighted correlation network, correlation and random forest analyses. Regression models and ROC curve analysis were used to evaluate the significance of identified molecular signatures in HT.

Results

Integrative analysis revealed subtle changes in gut microbiota diversity and composition in early HT, increased abundance of Bacillota_A and Spirochaetota at the phylum level, and significant differences in 24 genera and 67 species. Ecological network analysis indicated an imbalance in the gut microbiota with reduced inhibitory interactions against pathogenic genera in HT. Functional analysis showed changes in infection- and immune-related pathways. Three characteristic species (Salaquimonas_sp002400845, Clostridium_AI_sp002297865, and Enterocloster_citroniae) were identified as most relevant to HT. Analysis of miRNA and mRNA expression profiles uncovered pathways related to immune response, inflammation, infection, metabolism, proliferation, and thyroid cancer in HT. Based on correlations with HT and interactions between them, six characteristic RNAs (hsa-miR-548aq-3p, hsa-miR-374a-5p, GADD45A, IRS2, SMAD6, WWTR1) were identified. Furthermore, our study uncovered significant gut microbiota-host transcriptome interactions in HT, revealing enrichment in metabolic, immune, and cancer-related pathways, particularly with strong associations among those 9 key molecular signatures. The validation stage confirmed improved HT classification accuracy by combining these signatures (AUC = 0.95, ACC = 0.85), suggesting their potential significance in understanding HT pathogenesis.

Conclusion

Our study reveals novel molecular signatures linking gut microbiome and host transcriptome in HT, providing new insights into the disease pathogenesis. These findings not only enhance our understanding of the gut-thyroid axis but also suggest potential new directions for therapeutic interventions in HT.

Background

Hashimoto's thyroiditis (HT) is a leading cause of hypothyroidism, characterized by an autoimmune attack on the thyroid gland [1]. It affects approximately 10–12% of the global population [2] and is influenced by factors like increased stress and dietary habits [3, 4]. The exact causes of HT are not fully understood but involve genetic, environmental, and epigenetic factors [5]. HT typically starts with the appearance of thyroid peroxidase antibodies (TPOAb) and/or thyroglobulin antibodies (TGAb), while thyroid function remains largely normal initially [6]. As the disease progresses, damage to thyroid follicular cells gradually occurs, leading to subclinical or clinical hypothyroidism [6]. The widespread use of autoimmune thyroid antibodies testing has increased the diagnosis of many early HT patients [6]. Unfortunately, effective treatments to prevent thyroid tissue damage and thyroid dysfunction in these patients are currently lacking [6].

Analyzing gut metagenome and transcriptome data is crucial for disease research. It helps unravel the complex interactions between the gut microbiota and host gene expression, providing insights into the mechanisms underlying disease onset and progression [7]. This approach has shown promise in identifying molecular signatures for various diseases [8, 9]. These advancements have revolutionized disease characterization and monitoring [10]. By applying these techniques to HT, we can uncover its pathogenesis and potential targets for treatment strategies. The "thyroid-gut axis" is crucial in understanding how the gut microbiota contributes to HT [11]. Evidence suggests that the gut microbiota plays a key role in triggering and advancing HT [12]. Advanced technologies like gut metagenome and transcriptome analysis have been used to study HT mechanisms and identify potential molecular signatures [13, 14]. While gut microbiota dysbiosis is recognized in HT [15, 16], there are still substantial knowledge gaps regarding its impact on host gene expression as measured by transcriptome analysis. This knowledge gap limits our understanding of the complex mechanisms of gut microbiota involvement in HT and hinders molecular signature identification. Additionally, comprehensive research is needed to integrate diverse data and gain a holistic understanding of HT's development, as well as to identify novel molecular signatures.

This study aims to uncover interactions between gut microbiota and host gene expression patterns in HT, identify novel molecular signatures, and enhance our understanding of its pathogenesis by integrating gut metagenome and host transcriptome data (miRNA/mRNA). We employed a two-stage case–control study design for robust results. Firstly, we analyzed the participants' gut microbiota composition using high-throughput sequencing techniques to explore its correlation with HT. Subsequently, we examined transcriptome data to identify differential gene expression regulation patterns in critical signaling pathways during the early stages of HT. By uncovering these distinct gene expression profiles, we gained valuable insights into the underlying mechanisms of HT progression. Finally, we integrated the gut metagenomics and host transcriptomics data to identify novel molecular signatures, which would serve as early characterization indicators and may suggest potential new directions for therapeutic interventions in HT.

Methods

Subjects

This study adopted a two-stage approach, utilizing two independent cohorts from different centers. From June 2022 to May 2023, we recruited 14 early HT patients and 13 healthy individuals from Qilu Hospital of Shandong University (discovery stage), and 17 early HT patients and 17 healthy individuals from the Second Affiliated Hospital of Xi'an Jiaotong University (validation stage). The natural course of HT starts with the presence of positive TPOAb or/and TGAb, often accompanied by morphological abnormalities in the thyroid, while thyroid function remains mostly normal. As the disease progresses, thyroid follicular cells gradually get damaged, leading to potential episodes of hyperthyroidism, but ultimately resulting in the development of hypothyroidism over time [6]. Thus, the diagnosis of early HT was based on the presence of thyroid enlargement, characteristic ultrasound images (such as hypoechogenicity or heterogeneous echo pattern), and elevated thyroid antibodies (TPOAb and/or TGAb), while maintaining normal thyroid function (normal TSH, FT3, and FT4 levels). Participants were included if they met these diagnostic criteria for early HT, were aged 18–65 years, and had no history of thyroid hormone replacement therapy. To further identify early HT, we recruited volunteers from the health examination centers of two hospitals for a six-month longitudinal assessment, with evaluations conducted every two months to observe changes in thyroid autoantibodies and ultrasound. During the initial assessment, thyroid autoantibodies had to be within normal levels. If both thyroid autoantibodies and ultrasound showed abnormalities during the second or third assessment, the participant's status at that assessment was included in this study, and peripheral blood and fecal samples were collected simultaneously. Participants were excluded if they experienced thyroid dysfunction during the evaluation period, had other thyroid diseases, a history of thyroid cancer or surgery, or any comorbid cardiovascular, autoimmune, infectious, musculoskeletal, or malignant diseases. Additionally, participants were excluded if they had a recent history of surgery or trauma, had used antibiotics, probiotics, or prebiotics within the past 3 months, or had significant changes in dietary habits in the past month. Healthy controls were age- and sex-matched individuals with normal thyroid function (normal TSH, FT3, and FT4 levels), negative thyroid autoantibodies (TPOAb and TGAb), and normal thyroid ultrasound findings.

Peripheral blood samples were collected after an overnight fast. Fecal samples were collected using sterile containers and immediately frozen at −80 °C until analysis. Thyroid function tests were performed using the Roche Cobas 6000 E601 module immunoassay analyzer (Roche, Basel, Switzerland). Demographic data were collected via a questionnaire. The study protocol was approved by the Medical Ethics Committees of both the Second Affiliated Hospital of Xi'an Jiaotong University (NO. 2020844) and Qilu Hospital of Shandong University (NO. 20200623SY), and conducted in accordance with the Helsinki Declaration. Written informed consent was obtained from each participant. Detailed demographic, clinical, and laboratory information is provided in Table 1.

Table 1 The demographic and clinical characteristics of the subjects

RNA sequencing and data analysis

Total RNA was extracted from peripheral blood samples of each participant using the PAXgene Blood RNA Kit (BD Biosciences, USA) following the manufacturer's instructions. Total RNA quality was assessed using agarose gel electrophoresis and quantified with a NanoDrop spectrophotometer (NanoDrop, USA). For mRNA-lncRNA-circRNA library construction, total RNA underwent ribosomal RNA depletion using the Epicenter Ribo-Zero kit. TruSeq RNA Sample Preparation kit processed 3 μg RNA/sample following Illumina's protocol. RT-PCR employed Phusion high-fidelity DNA polymerase, indexed (X) primers, and universal PCR primers. AMPure XP system purified the products, while library quality was evaluated on the Agilent Bioanalyzer 2100 system. The Illumina NovaSeq 6000 platform sequenced the RNA libraries, generating 150 bp paired-end reads. For miRNA library construction, 2 μg RNA/sample served as input material. TruSeq Small RNA Sample Preparation kit (Illumina, Inc., San Diego, CA, USA) generated miRNA sequencing libraries with added index codes for sample attribution. The library was PAGE gel-separated and the 140–160 bp region was selected. Library concentration was measured using the Qubit fluorometer, and quality assessed on the Agilent Bioanalyzer 2100 system. Indexed samples were clustered using TruSeq PE Cluster Kit v3-cBot-HS on the Illumina cBot clustering station. Illumina NovaSeq platform sequenced the libraries, generating 150 bp paired-end reads. Raw data quality control was performed using FastQC (version 0.11.8) [17]. Bowtie (version 1.1.2) was used to identify miRNAs [18], and the quantifier module of miRDeep2 generated miRNA expression profiles [19]. Duplicate samples' raw reads were processed using the R package edgeR with Trimmed Mean of M-values normalization [20]. Gene expression profiles were analyzed based on gene length, read count mapped to miRNAs, and fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) values were calculated. Clean reads were aligned to the reference genome using HISAT2 (version 2.0.4) [21], and StringTie (version 1.3.1) was used for transcript assembly [22]. Assembled transcripts were analyzed for coding potential with the RNA-seQc algorithm (http://www.broadinstitute.org/cancer/cga/rna-seqc). Cuffdiff (version 2.1.1) was employed to calculate FPKM for coding genes in each sample [23]. Gene FPKMs were determined by summing FPKMs for transcripts within each gene group.

Differential expression (DE) analysis of miRNA and mRNA from peripheral blood samples of HT patients and healthy subjects utilized the R package limma [24]. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis of DEmRNAs employed the R package clusterProfiler [25], while DEmiRNAs were analyzed using the web interface miEAA (https://ccb-compute2.cs.uni-saarland.de/mieaa2). [26] Weighted correlation network analysis (WGCNA) of DEmRNAs was performed using the R package WGCNA [27]. Heatmaps and Sankey diagrams illustrating clustering, pathways, and relationships were constructed using the R packages pheatmap and ggalluvial [28]. Data visualization was achieved using the R package ggplot2. [29]

Metagenomic sequencing and data analysis

DNA extraction from frozen feces samples using the QIAamp DNA Stool Mini Kit (QIAGEN, German) followed the manufacturer's protocol. The integrity, concentration, and quality of the DNA were determined using agar gel electrophoresis and NanoDrop spectrophotometry. PCR was used to amplify the sequencing libraries, which were then assessed for quality using a Qubit and Agilent Bioanalyzer. In the discovery stage, library construction failed for 3 HT case samples and 4 control samples. Therefore, most of the discovery stage samples (11 HT cases vs 9 controls) and all validation stage samples were used for metagenomic sequencing. High-throughput sequencing was performed on an Illumina NovaSeq platform, generating FastQ files with 150 bp paired-end reads. Raw data underwent quality control using FastQC (version 0.11.8). An average of 5,176,1899 clean reads per sample was obtained. Clean reads were assembled into Contigs using metaSPAdes (version 3.15.3) [30]. Bwa-mem (version 0.7.17) was used to compare clean reads to the Contigs [31]. MMseqs2 (version 13.45111) determined the sequence assembly efficiency [32], excluding contigs smaller than 500 bp. MetaGeneMark (version 3.38) predicted gene assembly efficiency [33], retaining genes with a length of ≥ 100 bp. R software calculated the expression transcripts per million (TPM) value of each Contig by comparing it to clean Reads using bwa-mem (version 0.7.17). TPM values were obtained for each gene by aligning nucleic acid sequences to the Genome Taxonomy Database for species annotations [34]. Functional annotations were obtained by aligning nucleic acid sequences to the Functional Protein Sequence Database using DIAMOND (version 2.0.14) [35].

α-diversity indices (Shannon, Ace, Simpson, Chao) and β-diversity analysis represented by non-metric multidimensional scaling (NMDS) were calculated using the R package vegan [36]. Stacked plots were generated using the R packages tidyverse and phyloseq to visualize the relative abundance of gut microbiota [37]. The Wilcoxon rank sum test (conducted by the R package doBy) was used to analyze the Firmicutes to Bacteroidetes (F/B) ratio and identify different gut microbiota [38]. Statistical analysis of KEGG Orthologs (KOs) level 3 pathways was performed using the software STAMP (version 2.1.3) [39]. Effects of environmental factors on the gut microbiota were assessed using redundancy analysis/canonical correspondence analysis (RDA/CCA) analysis and Monte Carlo permutations (999 permutations) test with the R package vegan [40]. Data visualization was done using the R package ggplot2 and GraphPad Prism (version 8.2.1). [29]

miRNA-mRNA network analysis

To establish a regulated network involving DEmiRNAs and DEmRNAs, we constructed a co-expression network of DEmiRNA-mRNA pairs. The mRNAs targeted by DEmiRNAs were obtained from databases such as "miRDB (http://mirdb.org)", "miRabel (http. //bioinfo.univ-rouen.fr/mirabel/index.php?page = mir)", and "rnainter (http://www.rnainter.org)". Further comparison was performed between the candidate target mRNAs and DEmRNAs to ensure their functional relevance in the network.

Correlation analysis

Spearman rank correlation analysis was employed to evaluate the relationships between gut microbiota and DEmiRNAs with clinical indicators, as well as between gut microbiota, DEmiRNAs, and DEmRNAs. This analysis was also utilized to construct ecological networks of the gut microbiota. Furthermore, Pearson rank correlation analysis was used to calculate correlations between WGCNA modules and genera. Both Spearman and Pearson rank correlation analyses were conducted using the R package Hmisc [41]. The visualization of results involved the use of R packages ggplot2 and graph, [29] as well as software tools Gephi (version 0.10.1) and Cytoscape (version 3.9.0). [42, 43]

Random forest analysis

Random Forest is a commonly utilized and highly effective machine-learning algorithm for screening molecular signatures. In this study, we utilized random forest analysis on the gut microbiota species to identify novel molecular signatures in HT. To ensure robustness, we performed internal validation using 5 iterations of tenfold cross-validation. The randomForest package in R was employed for conducting the random forest analysis [44]. By leveraging this approach, we assessed the predictive power of the gut microbiota species for HT characterization.

Assessment of molecular signatures

To assess the significance of identified molecular signatures in an independent validation cohort, we employed the R package autoReg to fit regression models. The capabilities of these molecular signatures in differentiating between early HT patients and healthy subjects were assessed using Receiver Operating Characteristic (ROC) curve analysis and accuracy (ACC) analysis. These analyses were conducted using the R packages pROC, rROC and caret. [45] The diagnostic potential of the molecular signatures were assessed by analyzing their discriminatory performance through multiple metrics, including receiver operating characteristic (ROC) curves, area under the curve (AUC), accuracy (ACC), sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV).

Results

Clinical characteristics of the subjects

No significant differences in age were found between groups (p > 0.05, Table 1) in both the discovery stage (14 early HT patients and 13 healthy subjects) and validation stage (17 early HT patients and 17 healthy subjects). In both cohorts, except for TPOAb and TGAb, which were remarkably elevated beyond normal levels in HT patients, other thyroid function indicators (FT3, FT4, and TSH) remained within the normal range (p values < 0.05, Table 1). This aligns with our definition of early HT, characterized by the presence of thyroid autoantibodies but maintaining normal thyroid function. The immune-mediated inflammation triggered by the elevated antibodies affected thyroid gland morphology, as evidenced by increased shear wave elastography (SWE) values and thyroid volume in early HT patients compared to healthy controls (Table 1). These findings are consistent with the characteristic ultrasound changes observed in early HT.

Dysbiosis of gut microbiota in individuals with early HT

Analysis of gut microbiota in both cohorts revealed no significant changes in α-diversity in early HT (p > 0.05, Fig. 1A), aligning with our focus on patients with elevated thyroid antibodies but normal thyroid function. Examination of β-diversity (Fig. 1B), stacked plots (Fig. 1C), and Firmicutes to Bacteroidetes ratio (Fig. 1D) showed suggestive compositional changes. These subtle alterations may reflect early immune responses and mild inflammation associated with early HT, affecting the gut environment without causing overt thyroid dysfunction. These results were consistent in both cohorts and the combined cohort (Fig. S1).

Fig. 1

Dysbiosis of gut microbiota in early HT. α-diversity of gut microbiota in healthy subjects and early Hashimoto's thyroiditis (HT) patients was assessed using Ace, Chao, Shannon, and Simpson indices (A). β-diversity analysis using non-metric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarity was applied. A stress value of < 0.1 indicated a reliable fit (B). Stacked plots displayed changes at the phylum level, presenting the top ten phyla based on relative abundance (C). The Firmicutes to Bacteroidetes (F/B) ratio was calculated using the Wilcoxon rank-sum test (D). The differential abundance of gut microbiota at the phylum (E), genus (F), and species (G) levels was determined using the Wilcoxon rank-sum test, with the top 24 genera and 30 species ranked by p-value magnitude (p < 0.05) (E–G). Spearman rank correlation analysis was used to construct ecological networks, depicting correlations between gut microbiota at the genus level (correlation > 0.6, p < 0.05) (H)

Despite subtle overall changes, significant differences were found at various taxonomic levels. Early HT patients showed increased abundance of Bacillota_A and Spirochaetota at the phylum level (p < 0.05) (Fig. 1E). At genus (Fig. 1F) and species (only top 30 species presented in Fig. 1G) levels, 24 genera and 67 species differed significantly (p < 0.05, Table S1). Notably, butyrate and acetate-producing genera (Catonella, Murimonas_intestini) and protective species (Barnesiella_intestinihominis) decreased [46, 47], while opportunistic pathogen Peptostreptococcus increased in early HT patients.

Ecological network analysis (Fig. 1H) revealed positive correlations among gut microbiota in both healthy and early HT individuals, with Bacillota_A playing a central regulating role. However, an imbalance was identified in early HT patients, with reduced inhibitory interactions against pathogenic genera (e.g., Klebsiella, Escherichia, Streptococcus), indicating disrupted gut microbiota. Enhanced inhibition was observed in genera like Alistipes (Bacteroidota) and Clostridium (Bacillota_A) (Fig. 1H).

Functional alternations in gut microbiota associated with early HT

Comparison of gut microbiota functional characterization using Welch's t-test revealed significant differences in five KO pathways between early HT patients and healthy subjects (p values < 0.05, Fig. 2A). Notably, Staphylococcus aureus infection [PATH:ko05150] and Coronavirus disease—COVID-19 [PATH:ko05171] pathways were enriched in early HT patients, both strongly associated with infection. Conversely, the natural killer cell mediated cytotoxicity pathway [PATH:ko04650], crucial in tumor immunity, was enriched in healthy subjects. Correlation analysis between species and pathways identified certain species from the Lachnospiraceae family as the closest relatives to these pathways (Fig. 2B). The genus Salaquimonas (family Rhizobiaceae) displayed the strongest correlation with natural killer cell mediated cytotoxicity [PATH:ko04650] (r = 0.74, p = 0.0002, Spearman rank correlation). Staphylococcus aureus infection [PATH:ko05150] was significantly associated with the genus CAG_95 (Lachnospiraceae) (r = 0.58, p = 0.0002), while Coronavirus disease—COVID-19 [PATH:ko05171] showed a strong correlation with the genus Duodenibacillus (family Burkholderiaceae_A) (r = 0.84, p < 0.0001) (Fig. 2C).

Fig. 2

Functional alterations in gut microbiota associated with early HT. Differential functional pathways of the gut microbiota between early HT patients and healthy subjects were assessed using Welch's t-test (p < 0.05) (A). Spearman rank correlation analysis identified correlations between 30 differential genera and functional pathways, represented by red (positive) and blue (negative) modules. Significance levels: *p < 0.05, **p < 0.01, ***p < 0.001. Those genera in the red box all belong to the Lachnospiraceae family (B). A scatter plot displayed genera with the highest correlation to pathways relevant to early HT (C)

Identification of molecular signatures for early HT by integrating gut microbiota and clinical indicators

Environmental factors analysis using RDA/CCA revealed that TPO-Ab, TG-Ab, SWE, TSH, age, FT3, FT4, and thyroid volume influenced the subtle differences in gut microbiota between early HT patients and healthy subjects. Among these, TPO-Ab, TG-Ab, TSH, and SWE showed significant correlations with gut microbiota structure (Fig. 3A). The positive correlation among these factors suggests their combined impact on gut microbiota composition. Early HT is characterized by elevated TPO-Ab and TG-Ab levels, while other thyroid-related parameters remain normal. This explains the limited impact on gut microbiota diversity and the non-significant changes observed. However, gut microbiota still holds promise as a molecular signature for early HT. The top 10 species significantly correlated with these clinical indicators were identified (p < 0.05, Spearman rank correlation) (Fig. 3B). Random forest analysis with tenfold cross-validation ranked Salaquimonas_sp002400845, Clostridium_AI_sp002297865, and Enterocloster_citroniae as the top three species with high importance metrics (Fig. 3C). The cross-validation curve was employed to examine the correlation between the model error and the number of species utilized in the fit. The results demonstrated that 2–5 significant species were sufficient to achieve optimal regression outcomes (Fig. 3D). These characteristic species may serve as molecular signatures in early HT, providing potential targets for early detection and intervention in HT development.

Fig. 3

Identification of molecular signatures for early HT using gut microbiota and clinical indicators. RDA/CCA analysis was utilized to examine the impact of environmental factors on the composition of gut microbiota in early HT patients and healthy subjects, with Monte Carlo permutations (999 permutations) determining the significance of this influence. CCA (Canonical Correspondence Analysis) and RDA (Redundancy Analysis) were employed, where sample groups were represented by different colors, and environmental factors were indicated by arrows originating from the origin. The choice between RDA and CCA was determined by the size of the first axis of Axis-lengths (a intermediate value during calculation, not shown here); a value ≥ 3.5 indicated CCA, while < 3.5 indicated RDA. Based on the intermediate value, we chose CCA for this analysis. The length of the arrows indicated the degree of influence of the environmental factors on the gut microbiota composition. The angle between the arrows reflected positive (acute angle) or negative correlation (obtuse angle), while the vertical distance from the sample point to the arrows represented the intensity of the influence on the samples. A shorter distance indicated a stronger influence. The right bar charts showed the significance of the correlation between each environmental factor (TPO-Ab, TG-Ab, SWE, TSH, age, FT3, FT4, and thyroid volume) and gut microbiota composition (A). Spearman rank correlation analysis identified the top 10 genera with the strongest correlation to clinical indicators, ranked by number of significant environmental factors and p-value magnitude. Positive and negative correlations were depicted by red and blue modules, respectively. Significance levels: *p < 0.05, **p < 0.01, ***p < 0.001 (B). Random forest analysis ranked species based on their importance in characterizing early HT (C). The cross-validation curve showed the correlation between model error and number of species used, indicating that 2–5 significant species are adequate for optimal regression outcomes (D)

Altered expression of the transcriptome in early HT

Comparison of miRNA and mRNA expression profiles between 14 early HT patients and 13 healthy subjects revealed 1975 down-regulated and 1821 up-regulated DEmRNAs. Additionally, 12 DEmiRNAs were up-regulated and 15 were down-regulated (|log2 FC|≥ 1, p < 0.05). The expression patterns of all DEmiRNAs are visualized in the heatmap (Fig. 4A, Table S2 and Table S3 for comprehensive lists of DEmiRNAs and DEmRNAs, respectively). GO and KEGG pathway analysis of DEmRNAs showed significant enrichment in pathways related to infection, metabolism, proliferation, and thyroid cancer (Fig. 4B). These findings provide insights into the biological and functional implications of transcriptome alterations in early HT, potentially contributing to the identification of molecular signatures associated with the disease's early stages.

Fig. 4

Altered expression of transcriptome in early HT. A volcano plot depicted the DEmiRNAs and DEmRNAs in early HT patients and healthy subjects, where blue dots represented down-regulated DEmiRNAs/DEmRNAs in early HT and red dots represented up-regulated ones (|log2FC|> 1, p < 0.05). A clustered heatmap displayed all up-regulated or down-regulated DEmiRNAs, with red modules indicating high enrichment and blue modules indicating low enrichment (A). The functional prediction of DEmRNAs was conducted using KEGG and GO analysis, presenting all differential KEGG pathways and the top 10 GO pathways (BP, CC, and MF), ranked by p-value magnitude (p < 0.05) (B)

Identification of molecular signatures for early HT by combining transcriptome analysis and clinical indicators

Significant correlations were observed between 27 DEmiRNAs and 7 clinical indicators in early HT patients, particularly with SWE, TSH, and antibodies (Fig. S2). The top 10 DEmiRNAs with highest correlations to clinical indicators were identified, with has-miR-548aq-3p showing the most correlations (Fig. 5A). GO pathway analysis of the top 10 DEmiRNAs revealed enrichment in multiple pathways, including T cell proliferation regulation and MAP kinase activity regulation, suggesting synergistic actions in early HT pathogenesis (Fig. 5B, p< 0.05). Has-miR-548aq-3p and has-miR-374a-5p were particularly enriched in immune- and inflammation-related pathways. Correlation analysis revealed significant associations among the 10 DEmiRNAs, except for hsa-miR-5094. This further supports the synergistic effect of these DEmiRNAs in their functional interplay (Fig. 5C).

Fig. 5

Identification of molecular signatures for early HT using transcriptome analysis and clinical indicators. Spearman rank correlation analysis identified the top 10 DEmiRNAs strongly correlated with clinical indicators, ranked by the number of significant environmental factors and the magnitude of p-values. Red and blue modules represented positive and negative correlations, respectively. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001 (A). Using miEAA, the enrichment and annotation of the 10 DEmiRNAs revealed GO pathways associated with immune and inflammation (B). Spearman rank correlation analysis was performed to calculate the correlations between the 10 DEmiRNAs. Positive correlations were represented by the color red, while negative correlations were represented by the color blue (C). A Venn diagram showed the intersection of predicted target mRNAs of hsa-miR-548aq-3p and hsa-miR-374a-5p with up-regulated DEmRNAs (D). The functional prediction of target mRNAs was conducted using KEGG and GO analysis, presenting all differential KEGG pathways and the top 10 GO pathways (BP, CC, and MF), ranked by p-value magnitude (p < 0.05). Four genes, namely ASTN2, MEX3B, YPEL2, and PTPN14, were not found to be associated with any significantly enriched pathways related to HT. Alternatively, their enriched pathways that were identified had only minimal or weak connections to HT (E). A Sankey diagram displayed the enrichment of 2 DEmiRNAs and their co-targeted 4 DEmRNAs in HT-related pathways (p < 0.05) (F)

A DEmiRNA-mRNA regulatory network was constructed for has-miR-548aq-3p and has-miR-374a-5p. We obtained target gene lists for each DEmiRNA using miRDB, miRabel, and rnainter, selecting the overlapped genes as the final target gene set (Table S4). Eight target mRNAs (ASTN2, MEX3B, YPEL2, PTPN14, WWTR1, IRS2, SMAD6, and GADD45A) overlapped between predicted targets and up-regulated DEmRNAs (Fig. 5D). These DEmRNAs were significantly enriched in pathways associated with early HT, including TGF-β, Hippo, and FoxO signaling pathways. Additionally, three pathways governing parathyroid, insulin, and GnRH exhibited significant enrichment (Fig. 5E). Four of the 8 target DEmRNAs (GADD45A, IRS2, SMAD6, and WWTR1) were enriched in HT-related pathways (Fig. 5F). Based on the interactions between these 4 DEmRNAs and 2 DEmiRNAs (has-miR-548aq-3p and has-miR-374a-5p), and their association with early HT, they could be considered transcriptomic molecular signatures in early HT.

Integration of interactions between gut microbiota and transcriptome in early HT

WGCNA analysis of DEmRNAs identified 12 significant modules (Fig. 6A, Table S5). Analysis of these modules and 24 significant genera revealed 13 sub-modules strongly associated with some microbiota (p < 0.05, Fig. 6A). KEGG pathway analysis of genes in these microbiota-associated modules showed enrichment in metabolism, immunity, and cancer-related pathways (Fig. 6B, Table S6). Significant correlations were observed between 27 DEmiRNAs and 24 differential genera. Salaquimonas showed strong correlations with most DEmiRNAs. Genera associated with up-regulated and down-regulated DEmiRNAs displayed differences in homogeneity. Pararoseburia and Harryflintia were predominantly associated with down-regulated DEmiRNAs, while Coprococcus was mainly related to up-regulated DEmiRNAs (Fig. 6C).

Fig. 6

Integration of interactions between gut microbiota and transcriptome in early HT. Pearson rank correlation analysis was conducted to examine the correlations between the 12 WGCNA modules of DEmRNAs and the 24 differential genera. The correlation coefficients and p-values were labeled within the modules. A total of 13 sub-modules displayed significant correlations (p < 0.05). Positive correlations were denoted by the color red, while negative correlations were represented by the color blue (A). Functional enrichment of DEmRNAs within 13 microbiota-associated modules was analyzed, revealing robust enrichment in pathways related to metabolism, immunity, and cancer (p < 0.05) (B). Spearman rank correlation analysis determined the correlation between 27 DEmiRNAs and 24 differential genera, with red and blue modules denoting positive and negative correlations, respectively (C). The association between the top 10 significant species (only 9 with interactions) and the top 10 DEmiRNAs with their target DEmRNAs was evaluated using Spearman rank correlation analysis, revealing miRNA-species correlations with p < 0.05, as well as miRNA-mRNA targeting relationships. However, it should be noted that the species Blautia_A_intestinalis did not show a significant correlation with the 10 DEmiRNAs and is therefore not shown in the analysis results (D)

These findings suggest that gut microbiota in early HT patients can influence various biological processes through cross-regulation with host genes. Notably, the immune pathway-enriched modules (yellow, Fig. 6A) and two characteristic DEmiRNAs (hsa-miR-548aq-3p, hsa-miR-374a-5p, Fig. 6C) showed correlations with Clostridium_AI and Salaquimonas genera, highlighting their potential regulatory roles in early HT. Visualization of interactions between the top 10 significant species (only 9 with interactions) and the top 10 DEmiRNAs with their target DEmRNAs (Table S7) revealed strong associations between 3 characteristic species (Salaquimonas_sp002400845, Clostridium_AI_sp002297865, and Enterocloster_citroniae) and 6 characteristic RNAs (hsa-miR-548aq-3p, hsa-miR-374a-5p, GADD45A, IRS2, SMAD6, and WWTR1). This suggests their potential as combined molecular signatures for characterizing early HT (Fig. 6D).

Validation of molecular signatures in an independent cohort

The performance of identified molecular signatures was evaluated in an independent validation cohort. When used independently, gut microbiota showed the best performance in characterizing early HT (Fig. 7A). Combining all molecular signatures (2 miRNAs, 4 mRNAs, and 3 bacterial species) resulted in the most significant improvement in early HT classification performance, achieving an AUC of 0.95 and an ACC of 0.85 (Fig. 7B). These results support the potential of integrating gut microbiota, transcriptome, and clinical indicators as molecular signatures and treatment targets for early HT.

Fig. 7

Validation of molecular signatures in an independent cohort. Discrimination capacity of each type of molecular signature in characterizing HT (A). Discriminative capacity of the combined molecular signatures in characterizing HT (B)

Discussion

In this study, we performed an integrative analysis of gut microbiome and host transcriptome in HT patients and healthy controls. We identified subtle alterations in gut microbiota composition and significant changes in host gene expression patterns in early HT. By investigating the correlations between gut microbiota, transcriptome (miRNA/mRNA), and clinical indicators, we uncovered a novel set of molecular signatures that effectively characterized early HT. Our findings revealed complex interactions between gut microbiota and host gene expression, providing new insights into the pathogenesis of HT. This integrative approach not only enhances our understanding of the gut-thyroid axis in HT but also suggests potential new directions for early molecular characterization and therapeutic interventions.

Our metagenomic analysis revealed subtle alterations in gut microbiota composition in early HT patients, despite normal thyroid function. While α-diversity remained unchanged, we observed specific compositional shifts, including increased abundance of certain phyla and a slight dysbiosis in the Firmicutes/Bacteroidetes ratio, consistent with previous studies. [48, 49] At genus and species levels, we noted a decrease in beneficial bacteria and an increase in opportunistic pathogens. Ecological network analysis revealed reduced inhibitory interactions against pathogenic genera in early HT patients. These findings suggest that even in early HT stages, the autoimmune process can influence gut microbiota composition, potentially contributing to disease progression. The relationship between gut microbiota and HT remains an active area of research, with larger studies needed to fully elucidate the impact of HT on gut microbiota diversity.

Our ecological network analysis revealed increased abundance of Spirochaetota and reduced inhibition of Pseudomonadota in early HT patients. Functional analysis showed enrichment of infection-associated pathways, suggesting a potential role for pathogenic bacteria in early HT. We found significant correlations between thyroid antibody levels and gut microbiota composition, consistent with previous studies. [50, 51] Elevated antibodies in early HT may contribute to gut microbiota dysbiosis, potentially accelerating disease progression. While gut microbiota has been proposed as risk signatures for HT, [52, 53] our study focused on species-level selection combined with transcriptome data for more accurate and specific signatures. Validation in an independent cohort achieved good classification based on three species (AUC = 0.88), further improved by incorporating transcriptome data (AUC = 0.95). Notably, Salaquimonas_sp002400845 showed potential as a therapeutic target due to its absence in HT patients and correlations with clinical indicators and HT-related pathways. Further investigation into supplementation of this species or its metabolites may offer new therapeutic strategies for early HT.

To better understand the gut microbiota's role in HT development, we investigated the interplay between gut microbiota and host gene expression patterns. Our analysis identified 27 DEmiRNAs in early HT patients, consistent with previous studies showing gut microbiota's influence on host gene expression through miRNA interactions. [54,55,56] We focused on 10 DEmiRNAs based on clinical indicators, finding they act synergistically [57] to regulate HT-associated target genes. These DEmiRNAs also showed significant correlations with gut microbiota. Two DEmiRNAs, hsa-miR-548aq-3p and hsa-miR-374a-5p, were particularly enriched in immune-inflammatory pathways, suggesting their potential involvement in HT. While these miRNAs have been identified as signatures in other diseases, [58,59,60] their role in HT was previously unexplored. Validation in an independent cohort showed modest classification performance for these DEmiRNAs alone, but combining them with shared target genes significantly improved accuracy. This emphasizes the value of using multiple molecular signatures and demonstrates the co-regulatory roles of gut microbiota, DEmiRNAs, and target DEmRNAs in early HT development.

Our study has limitations that future research should address. The small sample size may limit generalizability, necessitating replication with larger, more diverse cohorts for robust results. The inclusion of only female participants introduces a gender bias; future studies should include both sexes for a comprehensive understanding of the gut microbiota-HT relationship. Additionally, we only compared HT patients with healthy controls, without including other autoimmune diseases for comparison. While our multi-dimensional biomarker panel shows strong discriminatory power for HT, future studies including patients with other autoimmune conditions would be valuable to further validate the specificity of these markers and their potential utility in differential diagnosis. While we identified significant correlations, it's important to note that correlation does not imply causation. Further experiments, such as animal models or fecal transplantation studies, are needed to establish causal relationships and determine the direct impact of gut microbiota on host gene expression patterns and early HT pathogenesis.

In conclusion, our study revealed subtle gut microbiota alterations and significant transcriptome changes in early HT patients, highlighting their interaction and the potential influence of disrupted gut microbiota on HT development. We identified a novel combination of molecular signatures for early HT characterization, integrating gut microbiota species, miRNAs, and mRNAs. These findings provide new insights into HT etiology and suggest potential new directions for early molecular characterization and therapeutic interventions in HT.