Shared Biomarkers and Potential Mechanisms Between Hashimoto’s Thyroiditis and Recurrent Miscarriage Revealed by Transcriptomics Analysis

Introduction

Hashimoto’s thyroiditis (HT) is one of the most common autoimmune thyroid diseases.¹ Studies have shown that the incidence of HT is approximately 0.3 to 1.5 cases per 1,000 population per year and is gradually increasing worldwide.² The primary features of the disease include chronic inflammation of the thyroid tissue and the presence of thyroid autoantibodies (anti-TPO and anti-TG). This autoimmune response leads to hypothyroidism, which eventually progresses to euthyroidism in approximately 20% to 30% of patients.³ Although many patients may be asymptomatic, prolonged uncontrolled HT can lead not only to disruption of thyroid hormone levels but also to a variety of clinical complications, including pericardial effusion.³ The pathogenesis of HT is complex and influenced by genetic, environmental, immunologic, and epigenetic factors; however, the exact mechanisms involved remain unclear.⁴ In addition, HT is significantly associated with the development of various malignancies, especially papillary thyroid cancer, and chronic inflammation is thought to be a possible important trigger for thyroid cell transformation.⁵ As HT is an autoimmune disease with complex and diverse etiological mechanisms, there is no effective treatment strategy.⁶ Therefore, there is an urgent need to explore its pathogenesis in detail, and its association with other diseases is significant.

Recurrent miscarriage (RM) refers to two or more spontaneous miscarriages occurring before 20–24 weeks of gestation, affecting approximately 1% to 5% of pregnant women.⁷ The occurrence of RM is associated with multiple factors, including anatomical defects, chromosomal abnormalities, endocrine disorders, the immune microenvironment, and psychological factors.⁸ Although maternal immunologic factors, thrombotic tendencies, uterine anatomic abnormalities, and endocrine abnormalities are the major known etiologic factors, there is still no clear etiology for approximately 50% of RM cases.⁹ Studies have shown that women over 30 years old are at a higher risk of RM. Furthermore, as the number of miscarriages increases, the probability of a subsequent successful pregnancy in RM patients may further decrease.¹⁰ In addition, a high percentage of RM patients also have polycystic ovary syndrome, suggesting that RM may interact with multiple pathological mechanisms.¹¹ Therefore, an in-depth study of the etiology of RM and its mechanisms is essential to improve pregnancy outcomes and enhance female fertility.

In recent years, the association between HT and RM has gradually attracted the attention of researchers.However, the direct pathophysiological links and shared molecular mechanisms remain elusive. Although autoimmune dysfunction is recognized as an etiological factor in RM, and the presence of thyroid autoantibodies (eg, TPOAb) significantly increases the risk of miscarriage,¹² the mere presence of antibodies cannot fully account for the recurrent nature of pregnancy loss. Notably, not all patients with HT experience RM, nor are all RM cases associated with thyroid autoimmunity,¹³ suggesting the involvement of multifactorial interactions more complex than a simple antibody effect. While both diseases share a background of immune dysregulation and chronic inflammation, the specific pathways connecting thyroid autoimmunity to endometrial receptibility or maternal-fetal immune tolerance are yet to be defined.¹⁴ As a systemic autoimmune disorder, HT may induce a pro-inflammatory state characterized by a Th1/Th17 immune bias and cytokine imbalance (eg, elevated IFN-γ, IL-17),¹⁵ which could disrupt immune tolerance at the maternal-fetal interface, potentially serving as a key mechanism in RM pathogenesis.¹⁶ Furthermore, HT leads to hypothyroidism via chronic thyroid inflammation and impaired hormone synthesis/secretion.¹⁷ Given the role of thyroid hormones in trophoblast function and early pregnancy maintenance,¹⁸ hypothyroidism, including subclinical forms, is associated with an increased risk of RM. Current clinical management, however, often focuses solely on correcting thyroid hormone levels, overlooking these underlying interactive mechanisms. This gap between theory and clinical practice underscores the necessity to investigate the common genomic features that concurrently drive both conditions. Elucidating these shared etiological factors will not only provide deeper insights into the HT-RM synergy but may also reveal novel molecular targets for early risk prediction and targeted therapeutic interventions, ultimately improving pregnancy outcomes in this high-risk population.

In summary, to comprehensively understand the association between HT and RM and the possible mechanisms, this study was conducted to identify candidate genes co-expressed in HT and RM by transcriptomic analysis, combined with protein-protein interaction (PPI) network and expression level validation, to obtain biomarkers co-expressed in HT and RM, and to validate the biomarkers to explore the influence of biomarkers on the pathogenesis of HT and RM and the possible mechanisms, which will provide a new theoretical basis for the prevention and treatment of these two diseases.

Materials and Methods

Data Source

From the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds), transcriptome datasets about Hashimoto’s thyroiditis (HT) and recurrent miscarriage (RM) were obtained. The GSE138198 dataset (GPL6244 platform) included 13 HT (HT group) and 3 normal (control group) thyroid tissue samples,¹⁹ which served as a training set 1. Meanwhile, the GSE165004 dataset (GPL16699 platform), consisting of 24 RM (RM group) and 24 normal (control group) endometrial tissue samples was used as a training set 2. The GSE26787 dataset (GPL570 platform), consisting of 5 RM (RM group) and 5 normal (control group) endometrial tissue samples,²⁰ It was used for verification to assist in the screening of biomarkers. The GSE286332 dataset (GPL34281) consists of thyroid tissue samples from 9 RM patients (RM group) and 9 normal (control group) for external validation. The clinical information of each dataset is presented in Tables S1–S4.

Differential Expression Analysis

In the GSE138198 dataset, to identify differentially expressed genes in HT (HT-DEGs), the limma package (v 3.54.0)²¹ was used to perform differential expression analysis between HT and control groups (|log₂Fold Change (FC)| > 0.5, adjusted p-value (adj.p) < 0.05). Subsequently, the ggplot2 package (v3.4.4)²² and the complexheatmap package (v 2.18.0)²³ were applied to draw a volcano plot and a heatmap, respectively. The top 5 up- and down-regulated HT-DEGs, ranked by log2FC from high to low, were labeled in the volcano plot, and their expression trends were visualized in the heatmap. Similarly, in the GSE165004 dataset, RM-DEGs were identified between the RM and control groups.

Functional Pathways Analyses of Candidate Genes

Moreover, up-regulated HT-DEGs and up-regulated RM-DEGs were intersected using VennDiagram (v 0.1.10)²⁴ to obtain up-DEGs. Similarly, the down-regulated HT-DEGs were intersected with the down-regulated RM-DEGs to obtain the down-DEGs. Subsequently, the up-DEGs and down-DEGs were combined to generate a set of candidate genes.

To explore the biological roles of candidate genes, GO term annotation and KEGG pathway enrichment were performed carried out using the clusterProfiler package (v 4.10.0)²⁵ and org. Hs. eg.db package (v3.16.0)²⁶ (p < 0.05). The 10 GO pathways that were most significantly enriched were visualized.

Construction of Protein‒Protein Interaction (PPI) Network

To explore the interactions of the candidate genes at the protein level, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) website (https://cn.string-db.org/) was employed to establish a PPI network (interaction score > 0.15). Cytoscape software (v 3.10.2)²⁷ was employed to visualize the protein interactions. Subsequently, connectivity scores of the candidate genes were calculated using 2 algorithms (edge percolated component (EPC) and clustering coefficient) of Cytoscape software (v 3.10.2). The top 10 scored genes for each algorithm were intersected to obtain core genes.

Identification of Biomarkers

Expression variations of essential genes across sample groups were analyzed through Wilcoxon testing in three datasets: GSE138198, GSE165004, and GSE26787 (p < 0.05). The genes with remarkable differential expression across groups and consistent trends in 3 datasets were chosen as biomarkers. In addition, the biomarkers were also validated for their expression levels and subjected to ROC analysis using the external validation set GSE286332 dataset.

Construction and Evaluation of a Nomogram

In the GSE165004 dataset, a nomogram comprising biomarkers was constructed using the rms package (v 6.7.1)²⁸ to predict the likelihood of developing RM. By computing the total score, the likelihood of having RM was evaluated. In addition, a calibration curve was plotted using the ResourceSelection package (v 0.35).²⁹ Hosmer-Lemeshow (HL) p > 0.05 indicated that the predicted and actual values did not differ significantly, and the model fit was good. Besides, ROC curves were created by pROC (v 1.18.0)³⁰ to quantify the diagnostic effect of the nomogram (AUC > 0.7).

Localization and Gene-Gene Interaction (GGI) Network Analyses of Biomarkers

The location of biomarkers on human chromosomes was analyzed by the RCircos package (v 2.4.4).³¹ Afterwards, the corresponding proteins of the biomarkers were analyzed for subcellular localisation using mRNALocater (http://bio-bigdata.cn/mRNALocater/). Meanwhile, genes functionally related to biomarkers were projected in the GeneMANIA database (https://genemania.org/), and a GGI network was constructed.

Gene Set Enrichment Analysis (GSEA)

To delineate the biological mechanisms associated with which the biomarkers might be involved during disease progression, GSEA was carried out. The “C2: KEGG gene sets” from the MSigDB were regarded as the reference gene set. We conducted Spearman correlation analyses between each biomarker and other genes using the psych package (v 2.2.9)³² in samples from the GSE138198 and GSE165004 datasets. The correlation coefficients were ranked (from high to low). Then the GSEA for biomarkers was examined by the clusterProfiler package (v 4.10.0) (False Discovery Rate (FDR) < 0.25, p < 0.05, |Normalized Enrichment Score (NES)| > 1). The top 10 significant pathways of each biomarker were depicted individually using the enrichplot package (v1.18.4).³³

Immune Infiltration Analysis

The abundance of 28 immune cells³⁴ for all samples in the GSE138198 dataset was assessed using the ssGSEA method to investigate how immune cells impact the course of HT. Additionally, the Wilcoxon test was performed to assess the infiltration differences of these immune cells between the HT and control groups (p < 0.05). Finally, the psych package (v 2.2.9) was used to calculate the correlation between differential immune cells and the correlation between biomarkers and immune cells (|correlation coefficient (cor)| > 0.3, p < 0.05). Meanwhile, to explore immune cells affecting RM progression, the same immune infiltration analysis was performed on all samples in the GSE165004 dataset.

Regulatory Network Analysis

To understand the molecular control systems governing biomarkers, the elmmo database in the multiMiR package (v 1.20.0)³⁵ was used to predict microRNAs (miRNAs) targeting biomarkers. Moreover, the Starbase database (http://starbase.sysu.edu.cn/) was applied to predict long non-coding RNAs (lncRNAs). Finally, the lncRNA-miRNA-mRNA network was generated utilizing the Cytoscape package (v 3.10.2).³⁶

The transcription factors (TFs) were computationally predicted using the NetworkAnalyst database (https://www.networkanalyst.ca), followed by the construction of a TF-mRNA network using the Cytoscape package (v 3.1).

Drug Prediction and Molecular Docking

The DrugBANK Database (https://go.drugbank.com/) was employed to identify potential drugs that interact with biomarkers. The Uniprot (https://www.uniprot.org/) was then used to acquire the biomarkers’ protein structures, and the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) was used to obtain the drugs’ three-dimensional structures. The proteins and drugs were uploaded into the cb-dock website (https://cadd.labshare.cn/cb-dock/php/blinddock.php) to perform molecular docking and calculate binding energy (binding energy below −5 kcal/mol indicated high binding affinity).

Statistical Analysis

The R software (v4.2.2) was used for all analyses. Additionally, the Wilcoxon test was used to assess differences between groups, with a significance threshold of p < 0.05.

Animal Experiment

Animal

Twenty 4-month-old SPF-grade female BALB/c mice were procured from Henan Skbeis Biotechnology Co., Ltd. (Animal License No: SCXK (Henan) 2020–0005). They were raised at Guizhou Medical University and acclimatized for 7 days in an environment with a constant temperature of 22°C ± 2°C, a constant humidity of 45% - 55%, and 12 hours of light and darkness each. This study was approved by the Ethics Committee of Guizhou Medical University (Approval No.: 2403682, Date:June 2,2024), and the approval opinion clearly stated that the experiment could be conducted in accordance with the Guidelines for the Welfare and Ethical Review of Laboratory Animals (GB/T 35892–2018) (Chinese National Standard) and the Guidance for the Care and Use of Laboratory Animals (issued by the Ministry of Science and Technology of China in 2006).

Establishment and Grouping of HT-RM Mouse Model

Twenty mice were divided into a control group (n=10) and a model group (n=20). A HT mouse model was established by high iodine water drinking combined with subcutaneous injection of porcine thyroglobulin. High iodine water with a concentration of 0.6 g/L (sodium iodide as solute and distilled water as solvent) and antigen solution with a concentration of 1 g/L (thyroglobulin as solute and phosphate buffer as solvent) were prepared. The primary immunization emulsion (Freund’s complete adjuvant: antigen solution = 1:1, concentration of 0.5 g/L) and the booster immunization emulsion (Freund’s incomplete adjuvant: antigen solution = 1:1, concentration of 0.5 g/L) were prepared. The model mice drank high iodine water during the modeling period. After one week of adaptive feeding, in the second week, 200 μL of the primary immunization emulsion was subcutaneously injected at multiple points on the back of the model mice once a week, whereas the control group received an equivalent volume of normal saline at corresponding sites. In the third to fifth weeks, 200 μL of the booster immunization emulsion was subcutaneously injected at multiple points on the back of the model mice once a week, while control mice received equivalent amounts of normal saline at corresponding locations. In the ninth week, venous blood was obtained from mice in each group to detect the levels of TGAb and TPOAb. The success of the modeling was confirmed when the serum levels of TGAb and TPOAb in the model mice reached or exceeded the mean values of the control group.

ELISA Assay

Four weeks after the model was established, blood samples were collected from the orbital cavities of 20 mice (10 mice in each group) during the estrus period after inhalation anesthesia with isoflurane. Peripheral blood specimens were incubated at ambient temperature (22–25°C) for 120 minutes before refrigerated centrifugation (3000 × g, 20 min, 4°C) to obtain cell-free supernatant for subsequent analyses. According to the instructions, the levels of serum thyroid peroxidase antibody (TPOAb), thyroid globulin antibody (TGAb), thyroid stimulating hormone (TSH), and free thyroid hormones (FT3/FT4) were detected using an ELISA kit.

HE Staining of Thyroid Tissue

After the animals died peacefully, thyroid tissues were collected. Some of them were fixed with 4% paraformaldehyde for 24 hours, then dehydrated and embedded in paraffin for pathological analysis. The remaining thyroid tissues were stored at −80°C in a freezer. The paraffin blocks were sectioned continuously with a thickness of 5 mm. The sections were stained with hematoxylin and eosin (H&E), and then observed under an optical microscope for the morphology and structure of the thyroid.

qPCR Experiment

Take the thyroid and uterine tissues from each group, add 1–2 steel balls, then use the high-throughput tissue grinder (New Zhong Scientists-48) to grind them, and extract the tissue RNA according to the instructions of the TRIzol RNA extraction kit. Reverse transcribe the RNA into cDNA. Prepare the reaction system, conduct qPCR reactions, and measure the Cq value of gene expression levels to detect the expression of candidate biomarkers and cytokine genes related to miscarriage.

Immunohistochemistry Experiment

IHC detection of TGF-β was performed on uterine sections from both groups. After citrate-based antigen retrieval and peroxidase blocking (3% H₂O₂, 37°C, 20 min), sections were probed with primary antibody (4°C, overnight) and HRP-polymer secondary (37°C, 20 min). DAB visualization was followed by standard hematoxylin counterstaining and histological processing for microscopic evaluation.

Immunofluorescence Experiment

Immunofluorescence staining was performed to detect the expression of TGF-β in uterine tissues of the control group and the HT-RM group. Antigen retrieval was conducted using 0.01 mol/L citrate buffer solution. Each sample was treated with 3% hydrogen peroxide and incubated at 37°C for 20 minutes. A 5% bovine serum albumin (BSA) solution was generously applied to cover the ovarian tissue sections and incubated at room temperature for 2 hours. Subsequently, Primary antibody incubation was performed overnight at 4°C. Following primary antibody incubation, slides were retrieved and rinsed three times with PBS (10 minutes per wash). For secondary antibody incubation, Excess liquid around tissues was carefully removed using absorbent paper. A cocktail solution containing species-matched fluorescent secondary antibody (1:200 in PBS) and DAPI (1:1000 in PBS) was applied to the sections. Slides were then incubated in a humidified chamber at 37°C for 1 hour, protected from light. After incubation, the fluorescent antibody solution was discarded, and slides were washed three times with PBS (10 minutes per wash). For mounting (performed under light-protected conditions), An appropriate amount of antifade mounting medium was applied onto the tissue sections. Cover slips were carefully placed, ensuring no air bubbles remained within or around tissues. Confocal microscopy was performed using a laser scanning microscope system.

Statistical Analysis

Data analysis was conducted with GraphPad software. Two-group comparisons were evaluated using the unpaired t-test. Each experiment was performed independently in triplicate (n = 3). Statistical significance was defined as P < 0.05.

Results

Exploring the Biological Functions of 47 Candidate Genes

In the GSE138198 dataset, the differential expression analysis revealed the presence of 3770 HT-DEGs, including 1,775 up- and 1,995 down-regulated HT-DEGs in the HT group (Figure 1a–b and Table S5). Besides, there were 595 RM-DEGs between the RM and control groups in the GSE165004 dataset, including 220 up- and 375 down-regulated RM-DEGs in the RM group (Figure 1c–d and Table S6). Afterwards, 20 up-DEGs and 27 down-DEGs were identified (Figure 1e–f). Ultimately, the integration of up-DEGs and down-DEGs yielded a total of 47 candidate genes.

Figure 1 Heatmap of differentially expressed genes among samples and Venn diagram of candidate genes, as well as GO enrichment pathways. (a) Differentially expressed genes between HT patient tissue samples and control samples. (b) Heatmap of differentially expressed genes between HT sample groups. (c) Heatmap of differentially expressed genes among RM sample groups. (d) Volcano plot of differentially expressed genes between RM sample groups. (e and f) Differential genes with the co-expression of HT and RM.G The GO enrichment pathways of candidate genes.

Notably, the candidate genes were significantly enriched in 162 GO pathways, like “protein stabilization”, “negative regulation of chemokine production”, “ATP−dependent protein folding chaperone” (Figure 1g and Table S7). Moreover, candidate genes were significantly enriched in 1 KEGG pathway, namely “cytoskeleton in muscle cells” (Table S8). These enrichment results suggested that the candidate genes might play crucial roles in both HT and RM by modulating key biological processes and pathways.

CFL1 and TRAPPC1 as Biomarkers

The PPI network included 29 proteins with 42 distinct protein-protein interactions (Figure 2a). After that, a total of 8 core genes were obtained, namely SUGT1, CAPZA2, CFL1, TCP1, TRIM23, SYNM, TRAPPC1, and FBN1 (Figure 2b–d). Unexpectedly, CFL1 and TRAPPC1 expressions were significantly higher in disease groups compared to controls (p < 0.05) and consistent expression trends in 3 datasets (Figure 2e–g). Therefore, CFL1 and TRAPPC1 were considered as biomarkers in this study. The validation of the biomarkers’ expression levels in the external validation set was consistent with that in the training set (Figure S1a), and the ROC analysis results showed that the AUC values of both biomarkers were greater than 0.8 (Figure S1b–c). Both genes showed consistent and significant upregulation in thyroid and endometrial tissues, indicating that they may play similar pathogenic roles in different organs. This provides an unprecedented candidate target for the development of universal biomarkers that cross tissue barriers.

Figure 2 Screening of candidate biomarkers and validation of expression levels. (a) 47 candidate genes underwent PPI (protein-protein interaction) network analysis. Each circle represents a gene. The larger the dot, the redder the color. The thicker the line, the stronger the connectivity. (b) ClusteringCoefficient algorithm, and the top 10 gene interaction network. (c) EPC algorithm, TOP10 gene interaction network. (d) The top 10 intersecting genes of the two PPI algorithms. (e) Validation of Expression Levels of Candidate Genes - HT Training Set GSE138198. * indicates p< 0.05, ** indicates p< 0.01. (f) Validation of Expression Levels of Candidate Biomarkers-RM Training Set GSE165004. ***Indicates p< 0.001, ****Indicates p< 0.0001. (g) Validation of expression levels of candidate biomarkers - RM validation set GSE26787. ns indicates not significant, *Indicates p< 0.05, **Indicates p< 0.01.

Nomogram with Excellent Performance for Diagnosis of HT

In GSE165004 dataset, a diagnostic nomogram for RM was constructed based on CFL1 and TRAPPC1 (Figure 3a), and the calibration curve demonstrated the nomograms’ substantial predictive accuracy for RM (p= 0.473) (Figure 3b). The ROC curve emphasized the excellent diagnostic value of the nomogram, with an AUC of 0.927 (Figure 3c). It’s excellent predictive performance indicates that the biomarkers have strong clinical transformation potential, providing potential molecular targets for the early non-invasive diagnosis of RM.

Figure 3 Construction and evaluation of the flowchart. (a) Biological marker nomogram. **Indicates p< 0.01. (b) Calibration curve of the nomogram. The x-axis represents the probability predicted by the nomogram for developing RM disease, while the y-axis represents the actual probability of developing RM disease. (c) The ROC curve has the false positive rate (1 - specificity) on the x-axis and the true positive rate (sensitivity) on the y-axis. The AUC represents the area under the ROC curve.

Chromosomal Localisation, Subcellular Localisation and GGI Network Analyses of Biomarkers

The CFL1 was located on chromosome 11, and TRAPPC1 on chromosome 17 (Figure 4a). In addition, both CFL1 and TRAPPC1 are predominantly localized in the cytoplasm (Figure 4b). Moreover, a GGI network was constructed based on the 20 genes whose biomarkers were related to their functions. There was a certain correlation between them, such as CFL1-ERBB4 and TRAPPC1-LIMK1. The GGI network demonstrated that biomarkers were involved in some functions, such as “actin polymerization or depolymerization” and “regulation of protein modification by small protein conjugation or removal” (Figure 4c).

Figure 4 Location analysis of biomarkers. (a) Biological marker chromosome localization analysis. (b) Subcellular-Biomarker Localization Map. (c) GeneMANIA Network Analysis.

Functional Pathways of Biomarkers

In the GSE138198 dataset, the GSEA results showed that biomarkers CFL1 and TRAPPC1 were enriched in 82 and 71 pathways, respectively (Table S9). It was worth mentioning that these biomarkers were co-enriched in some pathways like “graft-versus-host disease”, “autoimmune thyroid disease”, and “antigen processing and presentation” (Figure 5a–b). This co-enrichment suggested that CFL1 and TRAPPC1 might play significant roles in the pathogenesis of HT by modulating critical immune-related processes. For instance, the “autoimmune thyroid disease” pathway highlighted the involvement of these biomarkers in the chronic inflammation and immune dysregulation characteristic of HT. Additionally, their involvement in “antigen processing and presentation” underscored the potential impact on immune recognition and response, which were central to the autoimmune nature of HT.

Figure 5 GSEA enrichment analysis. (a and b) GSEA enrichment analysis of the CFL1 gene and the TRAPPC1 gene in GSE138198 (HT). (c and d) GSEA enrichment analysis of the CFL1 gene and the TRAPPC1 gene in GSE165004(RM).

In the GSE165004 dataset, the GSEA results showed that biomarkers CFL1 and TRAPPC1 were enriched in 50 and 52 pathways, respectively (Table S10). It was worth mentioning that these biomarkers were co-enriched in several pathways like “ribosome”, “Huntington’s disease”, and “cell adhesion molecules (cams)” (Figure 5c–d). These co-enriched pathways highlighted the potential roles of CFL1 and TRAPPC1 in the pathogenesis of RM.

Association of Biomarkers with the Immune Microenvironment

In the GSE138198 dataset, the 28 immune cells’ abundance in the HT and control groups was explored (Figure 6a), and the abundance of 13 immune cells was found to differ significantly (p < 0.05) (Figure 6b). These 13 differential immune cells were significantly more infiltrated in the HT than the control group. Moreover, significant positive correlations were found between most of the 13 differential immune cells (Figure 6c). The biomarkers CFL1 (cor = 0.738, p < 0.01) and TRAPPC1 (cor = 0.526, p < 0.05) were all significantly positively correlated with monocytes (Figure 6d). This finding underscored the importance of immune cell dysregulation in HT pathogenesis, as immune cells infiltration were a hallmark of the disease, contributing to thyroid tissue damage and inflammation.

Figure 6 Immune infiltration analysis. (a) Heatmap showing the infiltration of immune cells between the RM group and the normal group. (b) Differences in immune cells between the RM group and the normal group. *Indicates p< 0.05, **Indicates p< 0.01, ***Indicates p< 0.001. (c) Analysis of the correlation of immune cells between the RM group and the normal group. *Indicates p< 0.05, ***Indicates p< 0.001. (d) Analysis of the correlation between immune cells and biomarkers between the RM group and the normal group. *Indicates p< 0.05, **Indicates p< 0.01, ***Indicates p< 0.001. (e) Heatmap showing the infiltration of immune cells between the HT group and the normal group. (f) Differences in immune cells between the HT group and the normal group. *Indicates p< 0.05, **Indicates p< 0.01. (g) Analysis of the correlation of immune cells between the HT group and the normal group. **Indicates p< 0.01, ***Indicates p< 0.001. (H) Analysis of the correlation between immune cells and biomarkers between the HT group and the normal group. ns indicates not significant, *Indicates p< 0.05, **Indicates p< 0.01.

In the GSE165004 dataset, the 28 immune cells’ abundance in the RM and control groups was explored (Figure 6e). 6 immune cells displayed significant infiltration difference (p < 0.05) (Figure 6f). For example, monocytes exhibited markedly higher infiltration levels in the RM group. In comparison, eosinophils exhibited markedly lower infiltration levels in the RM group. Most of these 6 differential immune cells exhibited significant positive correlations with each other, indicating a complex interplay of immune responses in RM (Figure 6g). In addition, CFL1 (cor = −0.529, p < 0.001) and TRAPPC1 (cor = −0.497, p < 0.001) were all significantly negatively correlated with eosinophils (Figure 6h). While TRAPPC1 (cor = −0.497, p < 0.001) showed significant positive correlations with monocytes only (cor = 0.463, p < 0.001). These results demonstrated the possible role of immune cell dysregulation in RM pathogenesis, as immune cell infiltration was known to influence maternal-fetal immune tolerance and pregnancy maintenance.

Complex Regulatory Network

The constructed lncRNA-miRNA-mRNA network contained 18 miRNAs targeting CLF1, 16 miRNAs targeting TRAPPC1, and 60 lncRNAs targeting the above miRNAs (Table S11−12) Complex interaction pairs such as CFL1-hsa-miR-107-NEAT1 and TRAPPC1-hsa-let-7a-5p-XIST were included in the lncRNA-miRNA-mRNA network (Figure 7a). In addition, a total of 122 biomarker-related TFs were predicted. The TF-mRNA network revealed that SPI and PPARG, two transcription factors, co-regulated CFL1 and TRAPPC1 (Figure 7b). These intricate regulatory networks highlighted the multifaceted molecular mechanisms underlying HT and RM. Specifically, the involvement of lncRNAs and miRNAs in these networks underscored their potential roles in modulating gene expression and immune responses, which were critical in the pathogenesis of both HT and RM. This not only reveals the upstream regulatory mechanism, but more importantly, identifies that the two key transcription factors, SP1 and PPARG, can jointly regulate the two biomarkers. This provides a novel and efficient potential therapeutic target for the development of intervention strategies targeting transcription factors.

Figure 7 Molecular regulatory network. (a) The lncRNA-miRNA-mRNA molecular regulatory network, pink circular nodes represent biomarkers, Orange squares represent interacting miRNAs, and blue squares represent lncRNAs that interact with miRNAs. (b) TF-mRNA (biomarker) regulatory network. The pink circles represent four biomarkers (mRNAs), and the blue squares represent TFs.

The HT-RM Model Mice Were Successfully Constructed

Measurement of thyroid-related hormones by ELISA revealed that compared with the control group, the model group exhibited significantly elevated levels of TGAb and TPOAb antibodies (Figure 8a–b), increased TSH (Figure 8c), and decreased FT3 and FT4 levels (Figure 8d–e). To further investigate structural changes in the thyroid gland, hematoxylin-eosin (HE) staining of thyroid tissues was performed. The control group displayed normal histological architecture, whereas the model group showed disorganized thyroid structure with follicular atrophy and lymphocyte infiltration between follicles (Figure 8f). Subsequently, qPCR validation was conducted for eight characteristic genes (SUGT1, CAPZA2, CFL1, TCP1, TRIM23, SYNM, TRAPPC1, FBN1). In thyroid tissues of the HT-RM group, SUGT1, CAPZA2, CFL1, TRIM23, and TRAPPC1 exhibited expression trends consistent with our previous validation results (Figure 8g). In uterine tissues of the HT-RM group, CAPZA2, FBN1, CFL1, TRIM23, and TRAPPC1 showed concordant expression patterns with prior findings (Figure 8h). Notably, CFL1 and TRAPPC1 were consistently upregulated in the disease group. Given reports of dysregulated cytokines in HT patients potentially linked to miscarriage, we analyzed mRNA levels of IL-2, IL-4, IL-10, IL-17A, IL-23, IFN-γ, and TGF-β in uterine tissues of HT-RM mice. Compared with the control group, the model group demonstrated increased expression of IL-2, IL-17A, IL-23, and IFN-γ, alongside decreased expression of IL-4, IL-10, and TGF-β (Figure 8i). Additionally, we concurrently detected two cytokines, TGF-β and IL-17A, using both immunohistochemistry (IHC) and immunofluorescence (IF) techniques. The results demonstrated that compared with the control group, TGF-β expression was significantly decreased in the HT-RM group (Figure 8j–k), while IL-17A expression was markedly increased (Figure 8L). These differences were statistically significant (Figure 8m–n). The predictions of bioinformatics were verified through in vivo experiments, which improved the reliability of the research conclusions.

Figure 8 Validation in animal models. (a–e) Serum levels of TPOAb, TGAb, TSH, FT3, and FT4 in HT-RM mice. (f) HE staining of thyroid tissues demonstrated follicular atrophy and inflammatory infiltration in HT-RM mice. (g–h) mRNA validation of biomarkers in HT-RM tissues. (i) Verification of the expression levels of immune-related cytokines at the mRNA level. (j) Immunohistochemical IHC image of TGF-β. (k and l) Immunofluorescence (IF) images of TGF-β and IL-17A. (m and n) Fluorescence intensity statistical plot of IL-17A and TGF-β.

Biomarker-Targeted Drug Prediction

In addition, 2 and 1 drugs targeting CFL1 and TRAPPC1 were identified, respectively. Molecular docking revealed that CFL1 exhibited good binding affinity with artenimol (−6.2 kcal/mol). In comparison, TRAPPC1 exhibited strong binding affinity with S-palmitoyl-L-cysteine (−5.4 kcal/mol) (Table 1 and Figure 9). These findings suggested that CFL1 and TRAPPC1 might serve as potential therapeutic targets for HT and RM. The identified drugs, artenimol and S-palmitoyl-L-cysteine, could potentially modulate these pathways and alleviate the pathogenic mechanisms underlying HT and RM. Future studies should further explore the therapeutic potential of these compounds in preclinical models to validate their efficacy and safety.

Table 1 Statistics on Biomarker-Drug Binding Energy

Figure 9 Biological marker - drug molecule docking diagram.

Discussion

Integrated bioinformatics analysis coupled with experimental validation identified CFL1 and TRAPPC1 as shared biomarkers between Hashimoto’s thyroiditis (HT) and recurrent miscarriage (RM). Multiple lines of evidence support this central finding: both CFL1 and TRAPPC1 were significantly co-upregulated in tissues from HT and RM patients, occupied central positions in the protein-protein interaction (PPI) network, and were validated in an animal model. Furthermore, these genes were closely associated with disease-specific immune microenvironment alterations—specifically, monocyte infiltration in HT and eosinophil dynamics in RM. Functional enrichment analysis revealed their involvement in immune-inflammatory pathways and pathways critical for embryonic development, such as cell adhesion and ribosomal function, respectively. These results systematically demonstrate that CFL1 and TRAPPC1 may serve as key molecular bridges linking thyroid autoimmunity to dysfunction at the maternal-fetal interface, providing crucial theoretical and experimental insights into the comorbidity of HT and RM.

The roles of the CFL1 and TRAPPC1 genes in HT and RM have not yet been fully clarified. CFL1 (Cofilin-1) is a highly conserved actin-binding protein that is widely present in animal cells and mainly participates in the dynamic regulation of the cell skeleton.³⁷ CFL1 can bind to and regulate the polymerization and depolymerization of actin filaments, thereby influencing key processes such as cell morphology, migration, adhesion, division, and intracellular substance transportation. Moreover, CFL1 plays an essential role in various physiological and pathological processes,³⁸ such as in tumor cells,³⁹ where high expression of CFL1 is associated with tumor proliferation, migration, and invasion ability.⁴⁰ In the nervous system, CFL1 is involved in the pathological process of neurodegenerative diseases (such as Alzheimer’s disease). Dysfunction of CFL1 is also related to various diseases, including epilepsy, muscle atrophy, and vitamin B12 deficiency. Studies have shown that CFL1 can participate in actin cytoskeleton regulation through interaction with the subcortical maternal complex (SCMC), thereby influencing early embryonic development. However, no significant differential expression of the CFL1 gene was observed in the endometrium, chorionic villi, or decidual tissues of RM patients in the literature.⁴¹ Our findings demonstrate that in the HT-RM group, CFL1 was significantly upregulated in both thyroid and uterine tissues, showing a consistent trend with elevated levels of inflammatory factors IL-17A and IFN-γ. This phenomenon suggests that when RM is accompanied by HT, the function of CFL1 may shift from “embryonic development regulation” to “immune-inflammatory mediation.” Notably, current research on CFL1 in HT and RM remains limited, and its specific mechanisms require further validation. Therefore, CFL1 holds significant importance in cell biology research and has potential application value in the biomedical field, particularly in disease diagnosis, treatment, and drug development.

TRAPPC1, a core component of the TRAPP complex, is responsible for intracellular vesicular trafficking and protein transport, playing a critical role in autophagy and protein folding.⁴² TRAPPC1 is essential for various cellular functions,⁴³ including autophagy, secretion processes, and the proper folding of intracellular proteins. Variations in the TRAPPC1 gene are linked to a variety of diseases, including severe neurodevelopmental disorders and myopathies.⁴⁴ Moreover, the role of TRAPPC1 in maintaining the stability of the immune system suggests that it may be involved in the development of autoimmune diseases.⁴⁵ Although no current literature has reported on the role of the TRAPPC1 gene in HT and RM, the sustained upregulation of TRAPPC1 observed in both thyroid and uterine tissues of the HT-RM model may reflect a compensatory increase in intracellular protein synthesis and trafficking demands under autoimmune inflammatory conditions. Combined with the altered inflammatory cytokine profile observed in this study—characterized by elevated pro-inflammatory factors (IL-2, IL-17A, IL-23, IFN-γ) and reduced anti-inflammatory factors (IL-4, IL-10, TGF-β)—these findings suggest that TRAPPC1 may not only regulate intracellular trafficking but also modulate the local immune microenvironment by influencing the expression and secretion of immunoregulatory molecules. This finding offers novel understanding of the molecular mechanisms involved in RM among HT patients.

In HT, two biomarkers are co-enriched in “graft-versus-host disease”⁴⁶, “autoimmune thyroid disease”⁴⁷ and “antigen processing and presentation”, These enriched pathways are all closely associated with the activation and maintenance of autoimmune responses,⁴⁸ involving processes such as immune cell activation,⁴⁹ antigen recognition, and immune-mediated attack on self-tissues.⁵⁰ For example, the autoimmune thyroid disease pathway: HT is an autoimmune disease characterized by lymphocyte infiltration in thyroid tissue and the production of autoantibodies (such as thyroid peroxidase antibody TPOAb and thyroid globulin antibody TGAb). The enrichment of CFL1 and TRAPPC1 in this pathway indicates that they may participate in regulating the infiltration of immune cells and the production of autoantibodies,⁵¹ thereby affecting the inflammatory response of thyroid tissue; In RM, the two biomarkers are co-enriched in “ribosome”⁵², “Huntington’s disease” and “cell adhesion molecules”: Among them, the “ribosome” pathway is crucial for protein synthesis and intracellular homeostasis, and its dysregulation may affect embryonic development and implantation. In addition, the “cell adhesion molecules” pathway is involved in the interaction between the embryo and maternal tissues, which is crucial for successful implantation and the implantation of the placenta.⁵³ These findings suggest that CFL1 and TRAPPC1 may promote RM by regulating key biological processes related to embryonic development and maternal-fetal interaction.

Immunoinfiltration analysis in this study revealed significant differences in the infiltration abundance of 13 immune cell types in HT and 6 immune cell types in RM. Among these, monocytes play a pivotal role in the immune system, serving as key effector cells in maintaining homeostasis and defending against pathogens.⁵⁴ Previous studies have demonstrated significantly elevated IL-βlevels in peripheral blood mononuclear cells (including monocytes) of HT patients, indicating their critical involvement in disease pathogenesis.⁵⁵ Furthermore, monocytes in HT exhibit a pro-inflammatory signature due to aberrant gene expression, mediating immune tolerance breakdown and tissue injury. This mechanism shares commonality with other autoimmune diseases and operates independently of thyroid hormone regulation.⁵⁶ Notably, dysregulated monocyte chemotaxis, adhesion, and pro-inflammatory cytokine release constitute key pathogenic events in HT.

During normal pregnancy, monocytes infiltrate the decidua and differentiate into macrophages or dendritic cells. Decidual macrophages, in particular, contribute to fetal antigen tolerance and participate in critical processes such as placental development.⁵⁷ In contrast, RM patients exhibit an altered monocyte-associated immune microenvironment characterized by pro-inflammatory polarization, subset imbalance, and dysregulated crosstalk with NK/T cells, collectively contributing to pregnancy failure.⁵⁸ Intriguingly, our study observed a positive correlation between TRAPPC1 and monocytes in both HT and RM groups. We thus hypothesize that TRAPPC1 may participate in shared pathological mechanisms of HT and RM by modulating monocyte activation states, cytokine secretion, or intracellular protein trafficking functions.

Although this study successfully identified CFL1 and TRAPPC1 as shared biomarkers for HT and RM through bioinformatics analysis and preliminarily verified their expression trends and partial immune characteristics via animal experiments, several limitations remain. First, the sample size of the transcriptomic datasets (eg, GSE138198) was limited, which may affect statistical power and the generalizability of the conclusions. Second, the findings are primarily based on mRNA expression levels and bioinformatic predictions. The protein expression, functional activity, and causal roles of CFL1 and TRAPPC1 in the diseases have not been directly validated through genetic manipulation experiments. Furthermore, the predicted regulatory networks and potential targeted drugs (eg, Artenimol) require further assessment of their efficacy and safety through in vitro and in vivo studies. Additionally, although lymphocyte infiltration in the thyroid and hormonal abnormalities were observed in the animal model, this complex model cannot fully recapitulate the complexity of human disease. Future research should expand clinical sample cohorts for multi-center validation and utilize technologies such as single-cell RNA sequencing and spatial transcriptomics to delineate the heterogeneity of the immune microenvironment. By constructing conditional knockout/overexpression models and conducting preclinical drug trials, future work should functionally elucidate the mechanisms of CFL1/TRAPPC1 in the pathogenesis of HT and RM and explore their translational potential for targeted intervention, ultimately providing a solid foundation for developing cross-disease diagnostic strategies and therapeutic approaches.

Conclusions

In summary, this study identified the common biomarkers CFL1 and TRAPPC1 in HT and RM diseases, which may affect cell function and immune regulation, further influencing the occurrence of HT and RM. The results provide a basis for understanding the pathogenesis of HT and RM, as well as for improving clinical diagnosis and treatment. However, due to the insufficient sample size when predicting the incidence of HT, the identified biomarkers still need to be verified through a large number of clinical experiments, and the predicted drugs also need to be further verified through a series of in vivo and in vitro experiments. In the future, the sample size can be further expanded, the experimental period can be prolonged, and other related signaling pathways and targets can be continuously explored to provide more ideas and methods for researching the mechanism and treatment of this disease.