Clinical Efficacy and Mechanism of Vitamin D2 in Treating Hashimoto’s Thyroiditis

Introduction

Hashimoto’s thyroiditis (HT) is one of the most prevalent autoimmune diseases, with the highest incidence rate among autoimmune thyroid disorders.¹ Studies have revealed the annually increasing incidence of HT. In China, the prevalence of HT ranges from 0.2% to 1.3%.^2,3 HT pathology mainly includes diffuse lymphocyte and plasma cell infiltration, interstitial fibrosis, glandular atrophy, and acidic degeneration of follicles.⁴ The onset of HT is indistinct; many patients are asymptomatic during the early stages, and the cardinal clinical feature is painless, diffuse goiter. Approximately half of the patients with HT develop hypothyroidism and require lifelong thyroid hormone replacement. The risk of papillary thyroid carcinoma in patients with HT is considerably higher than that in the general population. Furthermore, the presence of thyroid peroxidase antibody (TPOAb) and thyroglobulin antibody (TGAb) may lead to unexplained infertility, abortion, preterm delivery, and postpartum thyroid disorders.⁵ Therefore, the high incidence of HT confers immense burdens on patients, the families of patients, and society.

Vitamin D2 (ergocalciferol) from foods and vitamin D3 (cholecalciferol) from sunlight exposure collectively constitute vitamin D.⁶ Vitamin D further undergoes biotransformation into more biologically active metabolites. In the liver, vitamin D is hydroxylated by 25-hydroxylase into its major circulating form, 25-hydroxyvitamin D (25OHD). Thus, blood 25OHD levels reflect vitamin D status from both diet and sunlight exposure.⁷ Numerous studies have confirmed associations between HT and trace elements such as iodine and selenium;^8–11 however, the association between HT and vitamin D remains debatable. A recent review proposed that vitamin D deficiency may contribute to the pathophysiology, hypothyroidism, and autoimmunity of HT.¹² Moreover, a randomized controlled trial showed that vitamin D supplementation benefits HT remission; newly diagnosed patients with HT were randomly assigned to either an intervention group receiving weekly supplementation of 60,000 IU of vitamin D3 and 500 mg/day of calcium for 8 weeks or a control group receiving only 500 mg/day of calcium. The results of the trial showed that after 3 months, TPOAb decreased by 46.73% in the intervention group, which is significantly more than the 16.6% decrease in the control group (P = 0.028).¹³ Additionally, Muscogiuri et al found a higher prevalence of autoimmune thyroiditis in subjects with vitamin D <20 ng/mL versus those with ≥20 ng/mL vitamin D. TPOAb levels negatively correlated with vitamin D3.¹⁴ In a large cohort, Choi et al observed lower vitamin D3 levels in women positive for TPOAb versus those negative for TPOAb.¹⁵ Furthermore, vitamin D receptor polymorphisms, such as BsmI and TaqI, may affect autoimmune thyroiditis risk.¹⁶

Recently, some studies have reported that patients with HT presented lower vitamin D levels, suggesting that vitamin D2 may exert therapeutic effects on HT.^17–19 Reportedly, higher serum 25OHD levels have been associated with a decreased risk of HT. A previous study showed that for every 5 ng/mL increase in the serum 25OHD level, the risk of HT incidence decreases by 19%.²⁰ However, presently, the molecular basis of vitamin D2 therapy for HT remains unclear and debatable,^21–23 highlighting the requirement for comprehensive research on the same.

This study aimed to further elucidate the molecular mechanisms underlying vitamin D2 therapy for HT. We collected peripheral blood samples from patients with HT before and after vitamin D2 treatment and performed transcriptome sequencing analysis to compare the changes in gene expression profiles after the treatment. The results of this study may provide the key genes and signaling pathways associated with vitamin D2 therapy for HT, elucidate the underlying molecular mechanisms, and provide a theoretical basis for the pharmacological treatment of HT.

Methods

Data Collection

In total, 10 patients with HT were enrolled as study subjects from the Department of Endocrinology, General Hospital of Ningxia Medical University and randomly divided into the control group and the experimental group, with five cases in each group. Peripheral blood transcriptome sequencing was performed for Hashimoto’s basic treatment group (Before group) versus Hashimoto’s basic treatment combined with vitamin D2 treatment group (After group). Inclusion criteria: age, 18–70 years; patients with a confirmed diagnosis of HT; patients with hypothyroidism who need to correct thyroid function to the normal levels for 3 months; and patients who voluntarily signed informed consent and could follow medication instructions and regular follow-up visits. Exclusion criteria: patients with discomfort due to enlarged thyroid; women who are pregnant, planning for pregnancy, or lactating; patients with other thyroid diseases; and patients with other severe systemic diseases. The diagnostic criteria used for HT were as follows: diffusely enlarged thyroid; firm texture; enlarged isthmus and lobes, which were confirmed by a professional physician; and positive serum TPOAb and TGAb. In all cases, the assessments, including the determination of thyroid enlargement, were conducted by experienced healthcare professionals.

All study subjects were required to have an iodine-restricted diet, mood smoothening, and regular work and rest. Patients with hypothyroidism were allowed to continue using levothyroxine sodium tablets. The Control group was treated with the selenium yeast tablet intervention method; two tablets (100 μg) each were orally administered 30 minutes after breakfast and dinner every day for 12 weeks. The Experimental group was treated with the vitamin D2 intervention method, in which, the medicine was dissolved in approximately 200 mL of boiling water and taken orally 30 minutes after breakfast and dinner while it was still warm; patients were administered medicine for 3 weeks, stopped for 1 week, and followed up, and the process was repeated thrice during 12 weeks. This study was approved by the Ethics Committee of General Hospital of Ningxia Medical University. All patients signed informed consent. This study strictly adheres to the ethical principles outlined in the Declaration of Helsinki, ensuring the adequate protection of the rights and well-being of the participants.

RNA Sequencing (RNA-Seq)

After 12 weeks of treatment, 2 mL peripheral blood samples were collected from the participants in the morning after overnight fasting to extract RNA for transcriptome sequencing. Total RNA was extracted using TRIzol reagent (Invitrogen, Burlington, ON, Canada).²⁴ RNA integrity was assessed using the RNA Nano 6000 Assay Kit and Bioanalyzer 2100 system (Agilent Technologies, CA, USA).²⁵ The qualified samples were subjected to RNA-Seq analysis at Novogene Corporation (Beijing, China).²⁶

Differential Expression Analysis

Fragments per kilobase of transcript per million mapped reads (FPKMs) were used to estimate the gene expression.²⁷ DESeq2 was utilized for determining differential expression in the digital gene expression data, and those with adjusted P-values of <0.05 were considered differentially expressed genes (DEGs).²⁸ In addition to the FPKM hierarchical clustering analysis of DEGs, we also analyzed the subclusters based on log2 (ratios) values of their gene expression relative to that in the control group. The log2 (ratios) values of ≥1 or ≤−1 were used as a cutoff for subcluster analysis in the Before and After groups. The clustering algorithm divided DEGs with similar gene expression trends into several subclusters.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analyses of DEGs

GO and KEGG (http://www.genome.jp/kegg/) enrichment analyses of DEGs were performed by the clusterProfiler R package, in which, the gene length bias was corrected. Those with a P-value of <0.05 were considered significantly enriched by DEGs.²⁹

Protein–Protein Interaction (PPI) Network Analysis

A combined score of ≥0.4 was set as the cut-off criterion for DEGs to construct the PPI network. The PPI network was visualized using the Cytoscape software (version 3.8.0).³⁰

Single Sample Gene Set Enrichment Analysis (ssGSEA)

The ssGSEA method is an extension of the gene set enrichment analysis method that allows defining an enrichment score representing the absolute enrichment of a gene set within each sample in a given dataset.³¹ Through ssGSEA, we determined the correlation among the samples and immune cells in the Before and After groups, and the variability in immune cell scores between these two groups.

Quantitative Reverse Transcription Polymerase Chain Reaction (RT‑qPCR)

RT-qPCR was used to detect 10 hub gene mRNA levels. Total RNA was extracted from blood samples using the TRIzol method, under the following standard reaction conditions: pre-denaturation at 94 °C for 10 min, denaturation at 94 °C for 1 min, annealing at 47 °C for 1 min, extension at 72 °C for 1 min, with 30 cycles, and finally extension at 72 °C for 10 min. The amplified products were analyzed by 1.5% agarose gel electrophoresis. GAPDH served as the internal control for 10 hub gene expression. The following primer sequences were used: CXCL8, forward: 5′-AAGTTTTTGAAGAGGGCTGAGA-3′, reverse:5′-AACCAAGGCACAGTGGAACA-3′, CXCL1, forward: 5′-AGGCCCTGCCCTTATAGGAA-3′, reverse: 5′-AAGGTAGCCCTTGTTTCCCC-3′; MMP9, forward: 5′-ACGATGACGAGTTGTGGTCC-3′, reverse: 5′-TCGCTGGTACAGGTCGAGTA-3′; PTGS2, forward: 5′-ATCTACGGTTTGCTGTGGGG-3′, reverse: 5′-TTCTGTACTGCGGGTGGAAC-3′; FOS, forward: 5′-AGACTACGAGGCGTCATCCT-3′, reverse: 5′-AGGTTGGCAATCTCGGTCTG-3′; EGR1, forward: 5′-ACCCCTTGCTCCCTTCAATG-3′, reverse: 5′-GGTGAGCATGTCCCTCACAA-3′; FPR1, forward: 5′-GGTGAACAGTCCAGGAGCAG-3′, reverse:5′-GGCAACGGCCACATTTATCC-3′; FCGR3B, forward: 5′-CTGTCTAGTCGGCTTGGTCC-3′, reverse: 5′-GGGGTCATGTGTCTTGAGGG-3′; TREM1, forward: 5′-CCAAGCTCCACCCAAGTCAA-3′, reverse: 5′-GCTGGTGGAATGAAGGACCA-3′; CCR1, forward: 5′-GTAGTGGTGACTGTGGGCTC-3′, reverse: 5′-TCCAAGATGGCAGTCGGAAC-3′; GAPDH, forward:5′- CAGCCCCAGAGTGTGTATCC-3′, reverse: 5′-GAAGATGCGGTCACCTCACA-3′.

Western Blot

Protein samples (10 μg) were analyzed using sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) at a constant voltage of 120 V for 60 min. The samples were blocked with 5% non-fat milk powder and incubated overnight at 4 °C with primary antibodies against CXCL8 (1:1000, ab289967), CXCL1 (1:1000, ab307589), MMP9 (1:2000, ab76003), PTGS2 (1:1000, ab179800), FOS (1:1000, ab208942), EGR1 (1:10,000, ab194357), FPR1 (1:1000, ab113531), FCGR3B (1:1000, ab89207), TREM1 (1:100, ab90808), CCR1 (1:1000, ab233832), and GAPDH (1:1000, ab8245). After washing with TBST, the samples were treated with horseradish peroxidase-labeled secondary antibodies (1:5000, ab6741) and incubated at 25 °C for 1.5 h. ECL color development was visualized using an automatic chemiluminescence imaging analysis system (Tanon, Shanghai, China) and analyzed using ImageJ software (version 1.45 s, Wayne Rasband, National Institutes of Health, USA, public domain).

Receiver Operating Characteristic (ROC) Curve Analysis of Hub Genes

The R package pROC was used,³² along with the gene expression matrix and Hub gene list, to calculate the area under the curve (AUC) value of each hub gene. Additionally, we plotted the ROC curves to visualize the accuracy of the Hub genes in predicting the outcomes of the two treatment groups of HT.³³

Correlation Analysis Between Hub Genes and Immune Cell Infiltration

The CIBERSORT algorithm was used to evaluate the relative proportion of 22 types of immune cell infiltration for each sample (perm = 1000), and the results were filtered according to the P-value of <0.05.³⁴ We referenced the R packages limma, reshape2, ggpubr, and ggExtra to examine the expression values for each hub gene and immune cell infiltration matrix. Furthermore, we used the cort.test function to perform the Spearman correlation tests obtain results with statistically significant differences (P < 0.05), and visualize the results using correlation scatter plots.³⁵

Interactions Between Hub Genes and Target Drugs

The Drug–Gene Interaction Database (DGIdb, https://dgidb.genome.wustl.edu/) was used to access information regarding the association of genes with their known or potential drugs. The 10 hub genes were imported into the database, and the therapeutic drugs related to the core genes associated with those 10 hub genes were obtained and visualized using the Cytoscape software.³⁶

Construction of lncRNA-miRNA-mRNA (ceRNA) Regulatory Network Based on Hub Genes

To identify the interactions between lncRNAs and hub genes (mRNAs), the data of lncRNAs and mRNAs was combined with that of miRNAs to construct the ceRNA regulatory network. The Cytoscape software was used to visualize the interaction network.³⁶

Statistical Analysis

Data analysis was conducted by Statistical software SPSS version 26.0 (SPSS Inc., Chicago, IL), and the results were expressed as mean±standard deviation. Student’s t tests were applied to compare the two groups. p < 0.05 was deemed as the statistical significance.

Result

Differential Expression Analysis

We analyzed differently expressed mRNAs with the cutoff fold changes ≥ 1 or ≤ −1 along with p < 0.05 and false discovery rate (FDR) < 0.05 in HT basic treatment combined vitamin D2 group After group vs Before group. The analysis strategy and procedure of the current study are shown in Figure 1. A total of 102 DEGs (6 upregulated and 96 downregulated) were identified between the after and before group (Figure 1a). Heatmap and hierarchical clustering were used to determine the gene expression patterns in the group of the HT basic treatment combined vitamin D2 group (Figure 1b).

Figure 1 (a) Volcano Plot and (b) Heatmap analysis of differentially expressed mRNAs between samples from the Hashimoto’s basic treatment group (Before group) versus Hashimoto’s basic treatment combined with vitamin D2 treatment group (After group). They were differentially expressed with the cutoff fold changes ≥1 or ≤ −1 along with p < 0.05 and false discovery rate (FDR) < 0.05.

GO and KEGG Enrichment Analyses of DEGs

The numbers of GO terms significantly enriched by the upregulated and downregulated DEGs in the Before and After groups are shown in Figure 2a. The 10 richest GO terms were mainly enriched in “response to lipopolysaccharide”, “response to molecule of bacterial origin”, “cell chemotaxis”, “regulation of angiogenesis”, “positive regulation of angiogenesis”, “positive regulation of vasculature development”, “leukocyte chemotaxis”, “female pregnancy”, “response to chemokine”, and “cellular response to chemokine.” The eight richest KEGG pathways were mainly enriched in “Cytokine-cytokine receptor interaction”, “Kaposi sarcoma-associated herpesvirus infection”, “Lipid and atherosclerosis”, “IL-17 signaling pathway”, “Viral protein interaction with cytokine and cytokine receptor”, “Apoptosis”, “Staphylococcus aureus infection”, and “NF-kappa B signaling pathway” (Figure 2b).

Figure 2 (a) The bubble plot of the top 10 enrichment pathways with DEGs by GO analysis; (b) The bubble plot of the top 8 enrichment pathways with DEGs by KEGG analysis.

PPI Network Analysis

Based on the STRING database (version 11.0), we constructed a PPI network reflecting the functional association between DEGs and visualized the network analysis results shown in Figure 3a. To further explore, the Top 10 gene by Degree method in the network ranked were shown in Figure 3b, the most significant hub proteins were CXCL8, CXCL1, MMP9, PTGS2, FOS, EGR1, FPR1, FCGR3B, TREM1, CCR1.

Figure 3 (a) DEGs with a combined score ≥0.4 was set as the cut-off criterion to construct the PPI network; (b) Top 10 gene by degree method in the network ranked.

Correlations Between the HT Treatment Groups and Immune Cells

Through ssGSEA, we determined the correlations among the samples in the Before and After groups and immune cells, and the differences in immune cell scores between the two groups.

The differences in infiltration levels of 28 immune cells in the five samples each of the Before and After groups (quantified by ssGSEA enrichment scores) are shown in Figure 4a. Compared with the Before group, the immune enrichment scores of mast cell (P = 0.008), neutrophil (P = 0.016), and plasmacytoid dendritic cell (P = 0.016) were significantly decreased in the After group (P < 0.05) (Figure 4b).

Figure 4 (a) Heatmap and (b) Violin plot of the difference in immune cell enrichment scores between samples from the Hashimoto’s basic treatment group (Before group) versus Hashimoto’s basic treatment combined with vitamin D2 treatment group (After group). Red in the boxes represents the After group and blue represents the Before group.

Expression Validation of 10 Hub Genes

To further explore the difference in the expression of 10 hub genes treatment with Hashimoto’s basic combined vitamin D2 group, we examined the expression of 10 hub genes in the Before and After groups. The RT-qPCR results showed that CXCL8, CXCL1, MMP9, PTGS2, FOS, EGR1, FPR1, FCGR3B, TREM1, CCR1 were significantly downregulated in the After group compared with the Before group in Figure 5a–j (P < 0.05). In addition, we found that the 10 hub genes related proteins CXCL8, CXCL1, MMP9, PTGS2, FOS, EGR1, FPR1, FCGR3B, TREM1, CCR1 were significantly decreased in the After group compared with the Before group in Figure 6a–j (P < 0.05).

Figure 5 RT-qPCR was used to detect the mRNA expression level of the 10 hub genes CXCL1 (a), TREM1 (b), PTGS2 (c), FOS(d), EGR1 (e), CXCL8 (f), MMP9 (g), PCGR3B (h), FPR1 (i), CCR1 (j) in the Hashimoto’s basic treatment group (Before group) versus Hashimoto’s basic treatment combined with vitamin D2 treatment group (After group). **p < 0.01,***p < 0.001.

Figure 6 Western blot was used to detect the protein expression level of the 10 hub genes CXCL1 (a), TREM1 (b), PTGS2 (c), FOS(d), EGR1 (e), CXCL8 (f), MMP9 (g), PCGR3B (h), FPR1 (i), CCR1 (j) in the Hashimoto’s basic treatment group (Before group) versus Hashimoto’s basic treatment combined with vitamin D2 treatment group (After group). *p < 0.05, **p < 0.01, ***p < 0.001.

ROC Curves of 10 Hub Genes

Analyses of the ROC curves of the 10 hub genes revealed that based on these hub genes, it could be accurately predicted whether a given sample belonged to the Before or After group. The AUC values of these 10 hub genes were as follows: were as follows: CCR1 (0.920), CXCL1 (0.960), CXCL8 (0.960), EGR1 (0.960), FCGR3B (0.920), FOS (1.000), FPR1 (0.840), MMP9 (0.720), PTGS2 (0.960), and TREM1 (1.000) (Figure 7a–j). This suggests that the identified hub genes exhibited good diagnostic performance in differentiating different subgroups (AUC > 0.700).

Figure 7 ROC curves for CCR1 (a), CXCL1 (b), CXCL8 (c), EGR1 (d), FCGR3B (e), FOS (f), FPR1 (g), MMP9 (h), PTGS2 (i), and TREM1 (j).

Correlation Between Hub Genes and Infiltrating Immune Cells

The immune cell infiltration profiles in each sample of the Before and After groups are shown in Table 1. Figure 8a–o shows the correlations between the expression values of 10 hub genes and infiltration levels of various immune cells (P < 0.05). Among them, CCR1 and FCGR3B expression positively correlated with neutrophil infiltration (R = 0.83 and 0.83, P = 0.0027 and 0.0027, respectively); CXCL1 expression negatively correlated with M2 macrophage and cluster of differentiation 8⁺ T-cell infiltration (R = −0.68 and −0.7, P = 0.029 and 0.031, respectively); CXCL8, EGR1, FOS, and PTGS2 expression negatively correlated with M2 macrophage infiltration (R = −0.7, −0.68, −0.68, and −0.68, P = 0.024, 0.029, 0.029, and 0.029, respectively); FPR1 expression positively correlated with monocyte and neutrophil infiltration (R = 0.7 and 0.67, P = 0.031 and 0.036, respectively), whereas it negatively correlated with resting natural killer (NK) cell infiltration (R = −0.78, P = 0.012); MMP9 expression negatively correlated with naive B cell infiltration (R = −0.76, P = 0.016), whereas it positively correlated with neutrophil infiltration (R = 0.87, P < 0.001); and TREM1 expression positively correlated with neutrophil infiltration (R = 0.69, P = 0.027), whereas it negatively correlated with resting NK cell infiltration (R = −0.66, P = 0.044).

Table 1 Infiltration Profiles of Different Immune Cells in Each Sample in the Two Groups

Figure 8 (a-o) Scatter plots of correlation between (CCR1, CXCL1, CXCL8, EGR1, FCGR3B, FOS, FPR1, MMP9, PTGS2, and TREM1) expression and different immune cell contents, neutrophil, M2 macrophage, T cells CD8, monocyte, NK cell, naive B cell.

Regulatory Network of Hub Genes and Targeted Drugs

From DGIdb, we obtained the drugs interacting with the hub genes (Supplementary sTable 1) and constructed the regulatory network between the hub genes and targeted drugs using the Cytoscape software. The interactions are shown in Figure 9, green squares represent hub genes, blue ovals represent molecularly targeted drugs, and the arrows between them indicate their interaction types, including agonist, antagonist, inhibitor, antibody, and vaccine. The constructed drug regulatory network includes 8 hub genes, 108 molecular drugs, and 114 interaction relationship pairs.

Figure 9 Green nodes indicate Hub genes, blue nodes indicate target drugs and the arrows connecting them indicate interaction pairs.

CeRNA Regulatory Network

In total, 165 pairs of mRNA–miRNA relationships were obtained through the databases TargetScan, miRanda, and miRDB and a customized Perl script (Supplementary sTable 2). A total of 160 pairs of miRNA–lncRNA relationships were obtained from the database spongeScan and Perl scripts (Supplementary sTable 3). The ceRNA regulatory network was constructed using the Cytoscape software, and it includes 129 lncRNAs, 145 miRNAs, 8 mRNAs (hub genes), and 324 interaction relationship pairs (Figure 10).

Figure 10 Red oval nodes indicate mRNAs (hub genes), blue diamond nodes indicate lncRNAs, green triangle nodes indicate miRNAs, and the lines between them indicate interaction pairs.

Discussion

Both genetic and environmental factors combinedly contribute to the development of HT. However, the mechanism underlying the disease pathogenesis has not been elucidated. Presently, various studies have correlated vitamin D deficiency with HT development and showed that vitamin D2 supplementation can help alleviate HT.^15,37,38 This study employs transcriptome sequencing technology to delve into the molecular mechanisms of vitamin D2 treatment for HT. The results reveal significant molecular-level differences between patients treated with vitamin D2 and those without such treatment. The disparities primarily concentrate on 10 key genes, with analyses indicating a notable correlation between these genes and the level of immune cell infiltration. Furthermore, through an analysis of the interactions between hub genes and targeted drugs, it was observed that these differential genes play pivotal roles in the network. Interaction types encompass various forms, including agonists, antagonists, inhibitors, antibodies, and vaccines, providing diverse perspectives for a comprehensive understanding of the treatment mechanism. This not only offers novel insights into the molecular-level treatment of HT but also lays the groundwork for future research and potential therapeutic strategies.

Herein, many DEGs (102 DEGs; 6 upregulated and 96 downregulated) were identified among patients treated with vitamin D2 and those not treated with vitamin D2. These DEGs were mainly enriched in “Cytokine–-cytokine receptor interaction”, “IL-17 signaling pathway”, and “NF-kappa B signaling pathway.” The cytokine–cytokine receptor interactions mediate signaling communication between cells and are involved in processes such as immune regulation, inflammatory responses, and cell proliferation. The cytokine signal transduction is an important bridge connecting the immune system with other physiological systems and plays a key role in maintaining homeostasis; some drugs used to treat clinical conditions such as inflammation, cancer, and autoimmunity mimic or inhibit cytokine–receptor interactions as their mode of action.³⁹ The interleukin (IL)-17 signaling pathway is involved in the pathological injury process of autoimmune diseases.⁴⁰ IL-17 binds receptors on target cells, activating inflammatory pathways such as nuclear factor (NF)-κB and mitogen-activated protein kinase pathways, inducing the production of inflammatory factors such as IL-6 and tumor necrosis factor-α. Furthermore, IL-17 recruits neutrophils, aggravating tissue damage. In autoimmune diseases, aberrant immune activation causes IL-17 overproduction. IL-17 inhibitors can block inflammatory responses and tissue injury, providing potential treatment targets.⁴¹ The NF-κB signaling pathway is widely regulated in the immune system.^42–44

The results of ssGSEA revealed that the immune enrichment scores of mast cells, neutrophils, and plasmacytoid dendritic cells decreased in the vitamin D2 treatment group. Mast cells are important pro-inflammatory cells that rapidly release transmitters, such as histamine, to cause inflammatory responses.⁴⁵ Neutrophils are the main effector cells of non-specific inflammation.⁴⁶ Plasmacytoid dendritic cells promote T-cell proliferation and differentiation.⁴⁷ A previous study showed that tissue miR-142-3p shuttles from T lymphocytes to thyroid cells and is an important regulator of the T regulatory cell functions and thyroid cell activity in the pathogenesis of HT.²¹ We speculate that vitamin D2 may inhibit mast cell activation and release of granules, decreasing their accumulation in the inflamed tissues. Concurrently, it may inhibit chemotaxis and adhesion of neutrophils, thereby decreasing their infiltration to the inflammation sites. The decreased secretion of chemokines and antigen presentation may decrease neutrophil aggregation in the lymphoid tissues. Thus, vitamin D2 may inhibit the activation and function of these pro-inflammatory cells, either directly or indirectly, thereby decreasing their enrichment scores in the treated tissues. This indicates that vitamin D2 may inhibit the inflammatory response by downregulating the aggregation of these pro-inflammatory cells in the tissues.

Vitamin D2 treatment affected immune cell infiltration, suggesting its immunomodulatory role via immune cell function regulation. The ROC curve analysis showed that based on the 10 hub genes, the origin of the sample could be accurately predicted, suggesting that there was a notable differential expression pattern of these hub genes before and after vitamin D2 treatment. Herein, we investigated the differences in the expression of these 10 hub genes in the After and Before groups by RT-qPCR and Western blot experiments, respectively. We hypothesized that this differential expression was related to the immunomodulatory effects of vitamin D2. This hypothesis was tested by the correlation analysis between the hub genes and infiltrating immune cells, which revealed that the expression of these hub genes was significantly positively or negatively correlated with the infiltration of immune cells such as lymphocytes and macrophages. These findings suggested that the hub genes may be involved in the process of vitamin D2 regulation of immune cell infiltration. The expression of hub genes such as CXCL1, CXCL8, EGR1, FOS, and PTGS2 decreased after vitamin D2 treatment, and the changes in their expression were negatively correlated with the M2-type macrophage infiltration, suggesting their involvement in the vitamin D2-mediated inhibition of M2-type macrophage infiltration. M2-type macrophages mainly exhibit anti-inflammatory and tissue repair functions.^48–50 Vitamin D2 seems to inhibit the aggregation of M2-type macrophages at the inflammatory sites by downregulating the expression of the aforementioned hub genes. We speculate that CXCL8 is an important chemokine that recruits macrophages, EGR1 may play a role by activating key transcription factors required for M2-type macrophage differentiation, and FOS and PTGS2 may be involved in generating chemical signals that direct M2-type macrophages. By inhibiting the abovementioned hub genes, vitamin D2 reduces the signals that direct M2-type macrophages to reach the inflammatory tissues, thereby inhibiting their infiltration. These findings help further elucidate the molecular mechanism by which vitamin D2 regulates the macrophage function. CCR1 is a chemokine receptor on the neutrophil surface, which is involved in chemokine-mediated neutrophil movement and activation,⁵¹ and MMP9 is an important matrix metalloproteinase secreted by neutrophils, which degrades the basement membrane and extracellular matrix, aiding neutrophil chemotaxis and tissue infiltration.⁵² Vitamin D2 treatment may enhance neutrophil responsiveness to chemokines by upregulating CCR1 expression and promoting their aggregation to the inflammatory sites. Furthermore, vitamin D2 may also increase the expression and secretion of MMP9, enhancing the ability of neutrophils to degrade extracellular matrix and make them more accessible to the tissues and target sites. Therefore, the upregulation of CCR1 and MMP9 can promote neutrophil infiltration and aggregation, which is positively correlated with the increase in neutrophil infiltration. Altogether, vitamin D2 appeared to improve neutrophil chemotaxis and tissue penetration by synergistically upregulating the expression of CCR1 and MMP9, thus promoting neutrophil infiltration into the inflammatory site. FPR1 is a chemokine receptor that mediates monocyte chemotaxis and activation and can promote monocyte migration to the inflammatory sites.⁵³ FPR1 expression was positively correlated with monocyte infiltration, suggesting that vitamin D2 may have improved the infiltration and aggregation of monocytes by upregulating FPR1 expression, which may favor the function of monocytes in antigen presentation and clearance at the inflammation site. Further studies with functional validation experiments are warranted to determine the role of these hub genes in the aforementioned processes.

The hub gene-targeted molecular drug network includes 8 hub genes, which are significantly changed after vitamin D2 treatment and may be related to the immunomodulatory effects of vitamin D2. They are also related to 108 molecular drugs (small molecule drugs or biologics) acting on the targets encoded by the hub genes. Furthermore, the network presented 114 pairs of interactions between the hub genes and molecular drugs. These interactions included effectors that can activate hub gene-encoded proteins, antagonists that can inhibit the function of hub gene-encoded proteins, inhibitors that can inhibit hub gene expression or the function of the encoded proteins, monoclonal antibodies that can specifically bind to the hub gene-encoded proteins to inhibit their function, and vaccines that are developed against the targets encoded by the hub genes.^54,55 The construction of this network allowed for the search of molecular drugs that could act on the hub genes that can be potentially regulated by vitamin D2, providing targets and candidate compounds for the development of drugs that mimic the immunoregulatory function of vitamin D2. Additionally, the lncRNA–miRNA–mRNA (ceRNA) regulatory network consisted of 129 lncRNAs, 145 miRNAs, 8 mRNAs (hub genes), and 324 pairs of interactions. The ceRNA network is based on a regulatory mechanism in which lncRNAs, by attracting miRNA binding, exert a regulatory effect on the expression of hub genes.^56–58 The construction of this network allowed for the search of vitamin D2-regulated non-coding RNAs and the prediction of their regulatory relationship with the hub genes, thus, exploring the non-coding RNA level mechanism of vitamin D2 immunoregulation and providing a theoretical basis for the subsequent regulation of hub genes by non-coding RNAs.

This study systematically investigated the mechanisms associated with the use of vitamin D2 for treating HT at the transcriptomic level. Through differential expression and enrichment analyses, we comprehensively revealed the changes in immune characteristics after vitamin D2 treatment, identified immune-related hub genes with significant differential expression, and proposed a hypothesis that vitamin D2 can regulate different immune cells based on correlation analyses. This study has some limitations. The small sample size led to limited statistical accuracy and the conclusions should be verified by performing follow-up experiments. This study proposes potential mechanisms underlying vitamin D2 regulation of immunity through transcriptomic and bioinformatic analyses and provides important insights and theoretical basis for subsequent research. To summarize, this study is of immense significance in exploring the regulatory mechanisms of vitamin D2 on immunity.

Inevitably, our study also has some limitations. The primary constraint of this research is the lack of an investigation into the treatment efficacy at different stages of HT. Through examinations and tests conducted at various treatment stages, further exploration of the molecular regulatory mechanisms of vitamin D in the treatment of HT could be undertaken. This could involve identifying differentially expressed genes, exploring pathways associated with differential gene enrichment, aiming to elucidate its predictive value in practical clinical applications. Additionally, clinical outcomes are typically multidimensional, encompassing aspects such as survival, disease progression, and treatment response. In the future, we anticipate conducting an in-depth investigation into the molecular regulatory mechanisms of vitamin D in the treatment of HT from multiple dimensions. This approach aims to comprehensively affirm the practical clinical value of vitamin D in the treatment of HT.

Conclusion

This study identified significant differences in molecular level between patients who received vitamin D2 treatment and those who did not. These key differentially expressed genes exhibit a pronounced correlation with the level of immune cell infiltration. From a molecular perspective, the study elucidates the regulatory role of supplementing vitamin D in alleviating and treating HT. This not only provides novel insights into the molecular-level treatment of HT but also establishes a foundation for future research and potential therapeutic strategies.