Introduction

Hashimoto thyroiditis (HT) has emerged as the most prevalent autoimmune disease.1 This condition is a lymphocyte-related chronic inflammation of the thyroid gland.1,2 It is reported that more than 20%–30% of the population suffers from this disease, with a continually increasing incidence.3,4 The most prevalent feature of Hashimoto’s thyroiditis (HT) is the elevated presence of two thyroid autoimmune antibodies: anti-thyroperoxidase antibody (TPOAb) and anti-thyroglobulin antibody (TGAb), which are detected in 95% of patients with HT and are rarely found in healthy individuals.5 Serum markers and ultrasonography remain the main diagnostic methods for HT.6–8 Although some patients are asymptomatic, prolonged and uncontrolled HT not only causes hypothyroidism or subclinical hypothyroidism in most cases but can also result in hyperthyroidism in some cases.9

HT has several other additional consequences. It has been reported that HT can, independently or synergistically with other diseases, contribute to severe cardiac effusion or cardiac tamponade.10 Furthermore, excessively high titers of TPOAb and TGAb may cross the blood-fetal barrier and lead to antenatal and neonatal disorders in pregnant patients.11 In addition, HT has been identified as a risk factor for thyroid cancer.12–14 Therefore, addressing HT at the source is crucial for reducing the prevalence and progression of thyroid cancer.

Currently, the main treatment for HT with permanent hypothyroidism or subclinical hypothyroidism is the daily, lifelong administration of oral levothyroxine. This treatment, however, addresses the symptoms rather than the underlying pathogenesis of HT.15,16 Since the mechanisms of autoimmune diseases are complex and multifaceted, there are no available strategies for curing HT. Therefore, it is essential for clinicians to seek effective therapies to treat HT at the level of pathogenesis and to prevent its progression before the onset of carcinogenic effects and other complications.

Bioinformatics, which integrates medicine and computational biology, contributes significantly to medical research by facilitating the exploration of disease pathogenesis and the identification of potential diagnostic as well as prognostic markers related to various tumors. However, there is a scarcity of published bioinformatics analyses related to HT. In the present study, we found that CD27 contributed to predicting clinical outcomes in patients with HT, including disease condition and response to immunotherapy. In addition, CD27 could potentially be used as a reliable biomarker for assessing the remodeling of the HT microenvironment. This could shed light on the underlying mechanisms of HT pathogenesis and progression, leading to new therapeutic strategies and improved treatment options.

Materials and MethodsData Selection

The training group datasets (GSE138198 and GSE54958) and test group datasets (GSE29315) were downloaded from the Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/). The GPL6244[hugene-1_0-st] Affymetrix human gene 1.0 ST array[transcript (gene version)] as well as the GPL8300[HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array were used for the training and test expression matrices, respectively.

The training group consisted of the GSE138198 dataset, which included 3 normal tissues and 13 hT tissues, and the GSE54958 dataset, which comprised 7 normal tissues. These were merged using the sva R package17 to correct for batch effects and non-biotechnological biases. The GSE29315 dataset, containing 8 samples of thyroid physiological hyperplasia and 6 hT samples, was selected as the test group.

Differential Analysis of Training Groups

DEGs between normal tissues and HT tissues within the training group were identified using the limma R package18 with the criteria ∣log2 FC∣ > 1.5 and an adjusted P value of < 0:05. Pheatmap and ggplot2 R packages were used for generating heatmaps and volcano plots of DEGs.

Enrichment Analysis of DEGs

Functional enrichment analysis of DEGs was conducted using the org.Hs.eg.db, enrichplot, clusterProfiler, DOSE, and GSEABase packages in R. GO functional enrichment analysis, KEGG pathway enrichment analysis, and Disease Ontology (DO) enrichment analysis were performed. The top 10 most significantly enriched functions in biological processes (BP), cell components (CC), and molecular functions (MF) were visualized by drawing bubble plots. The top 20 pathways or diseases with the most significant enrichment were selected for KEGG and DO enrichment analysis and visualized using bubble diagrams.

Weighted Gene Co-Expression Network Analysis (WGCNA)

A co-expression network was constructed using the weighted gene co-expression network analysis (WGCNA) package from R.19 First, the samples were clustered to identify significant outliers. Next, the automatic network constructor was utilized to develop the co-expression network. The soft thresholding power (β) was calculated using the pickSoft threshold function in R, and the common expression similarity was improved to determine the adjacency degree. Subsequently, modules were detected through hierarchical clustering and by using the dynamic tree-cutting function. Finally, the gene significance and module membership of the modules correlated with the HT groups were calculated. The gene information of the corresponding module was retained for subsequent analysis.

Screening of Hub Genes in HT

The intersection genes between DEGs and WGCNA module genes, particularly those highly correlated with HT groups, were identified. Then, a PPI network encompassing these intersection genes was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://cn.string-db.org/). The top 10 hub genes were then selected based on their degree of centrality using the cytoHubba function in Cytoscape 3.6.0 software.

Assessing the Clinical Diagnostic Value of Hub Genes

The differential expression analysis of the top 5 hub genes was conducted using the limma R package. Violin plots of the top 5 hub genes were generated using the ggpubr R package. Then, diagnostic ROC curves of the top 5 hub genes were constructed with the pROC R package. Nomograms were subsequently developed to predict the diagnostic value and risk assessment of these hub genes in HT using the rms R package.

Exploring the Potential Role of CD27 in HT Based on Functional Enrichment Analysis

The HT samples in the training groups were divided into the high- and low-CD27 expression groups based on the median value of CD27 levels. DEGs were identified using the limma R package under the criteria of ∣log2 FC∣ > 1.5 and an adjusted P value < 0.05. Subsequently, DEGs associated with CD27 were also analyzed. Utilizing the org.Hs.eg.db, enrichplot, clusterProfiler, DOSE, and GSEABase packages, GO functional enrichment analysis, KEGG pathway enrichment analysis, and GSEA enrichment analysis were performed. The top 10 most significantly enriched functions in BP, CC, and MF were visualized by drawing bubble plots. The top 13 pathways with the most significant enrichment were selected for KEGG enrichment analysis and visualized in a bubble diagram. Additionally, GSEA enrichment analysis was carried out on high- and low-CD27 expression groups, and the top five functions and pathways that were most significantly enriched were visualized.

Comprehensive Analysis of the Relationship Between CD27 and Immune Cell Infiltration in HT

To analyze the relationship between CD27 and immune cell infiltration in HT, the relative proportions of immune cells in each sample were quantified. The CIBERSORT algorithm, along with the e1071 and preprocessCore packages, were used to process the gene expression data and determine the proportion of different immune cells in each sample, which were subsequently visualized through bar charts.

The corrplot package was used to analyze the correlation of immune cells, and a correlation heat map was generated. Using the vioplot package, differences in immune cell infiltration between normal and HT tissues were analyzed, and the results were illustrated with violin plots.

Using the CIBERSORT and single-sample GSEA (ssGSEA) algorithms, the relationship between CD27 and immune-infiltrating cells was evaluated in order to explore the effect of CD27 on the immune microenvironment in HT. In addition, the relationship between CD27 and the functions of immune-infiltrating cells was assessed using ssGSEA algorithms.

Construction of the CD27 mRNA-Related Competing Endogenous RNAs (ceRNA) Network

For constructing the CD27 mRNA-related competing endogenous RNAs (ceRNA) network, the miRanda, miRDB, and TargetScanHuman online databases were utilized to calculate the CD27 mRNA-miRNA interaction pairs. Only the CD27 mRNA-miRNA pairings predicted by all three databases were retained for further analysis. The spongeScan online database was used to predict lncRNA-miRNA pairs based on the identified CD27 mRNA-miRNA pairings. The resulting network of interactions was visualized using Cytoscape 3.6.0 software.

Identification of CD27 Expression in Thyroid Carcinoma (THCA)

mRNA expression data for thyroid carcinoma (THCA) in fragments per kilobase of transcript per million mapped reads (FPKM) format were downloaded from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). The expression levels of CD27 in THCA were extracted from this dataset, and differential expression analysis between normal and tumor tissues was performed using the limma R package. Additionally, immunohistochemical staining figures of normal thyroid tissues and tumor tissues were retrieved from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) to visually compare CD27 protein expression in these tissues.

Clinical Sample Collection

Five clinical samples were obtained by ultrasound-guided needle biopsies from patients with HT, and three normal samples were obtained from patients with thyroid cancer who had undergone a total thyroidectomy. Prior to sample collection, all patients gave informed consent for the diagnostic procedures and the proposed treatment. All protocols were reviewed and approved by the Ethical Committee of The Second Affiliated Hospital of the Air Force Medical University.

Real-Time Quantitative PCR (Real-Time qPCR)

Total RNA was extracted from the collected samples using a Total RNA Extraction Kit (Beijing Solarbio Science & Technology Co., Ltd., Beijing, China). The concentration and purity of RNA were measured at A260 and A260/280, respectively. Following the manufacturer’s instructions, the RNA was then reversely transcribed into cDNA using HiScriptII Q Select RT SuperMix for qPCR (Vazyme Biotech Co. Ltd).

The primer sequences for CD27 and GAPDH were as follows: human CD27-Forward: 5′-GGACAAGCAGTGACCATCAAG-3′; human CD27 -Reverse: 5′-CCCAGAATTACCAAGTGAGTCCT-3′; human GAPDH-Forward: 5′- TGACTTCAACAGCGACACCCA −3′; human GAPDH -Reverse: 5′- CTACATGGCAACTGTGAGGAG −3′.

The expression levels of each gene were quantified using the 2−ΔΔCt equation, with GAPDH serving as the internal control for CD27 expression.

Measurement of CD27 Protein Levels in Clinical Samples

Total proteins were extracted from the clinical samples using radioimmunoprecipitation assay (RIPA) buffer (Beyotime, Shanghai, China), which was supplemented with a protease inhibitor cocktail (Sigma-Aldrich). CD27 was then measured by using a CD27 human Instant ELISA™ Kit (Invitrogen, BMS286INST) as per the manufacturer’s instructions. A volume of 100 μL of distilled water and 50 μL of the protein sample were added into each well of the micro-ELISA plate and incubated at room temperature for 3 hours. The microwell strips were washed six times with approximately 400 µL of Wash Buffer per well. Subsequently, 100 ul of TMB solution was added to the plates and incubated at room temperature for 15 minutes. Finally, 100 ul of the stop solution was added, and within two minutes of adding this, the absorbance of each well was measured at a wavelength of λ = 450 nm.

Statistical Analysis

All statistical analyses were performed using the R language (version 4.1). The data represented results from at least three independent experiments and were expressed as the mean ± SD. Additional statistical analyses were performed using SPSS 19.0 and mapping with GraphPad Prism 8.0. Group means were compared using the Student’s t-test for independent data. All P values were two-tailed, and a P value < 0.05 was considered to reflect statistical significance for the difference.

ResultsThe Flow Chart of This Study

As shown in Figure 1, the GSE138198 and GSE54958 datasets were utilized as training groups. Differentially expressed gene (DEG) analysis and weighted gene co-expression network analysis (WGCNA) were used to identify DEGs most associated with HT. Subsequently, hub genes were obtained by constructing protein-protein interaction (PPI) networks and further validated with the help of diagnostic receiver operating characteristic (ROC) curves and nomograms. CD27 was ultimately selected as a key gene in HT, and its relevance was confirmed using clinical samples. Then, groups with high and low expression of CD27 were subjected to Gene Ontology analysis (GO), Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis (KEGG), as well as Gene Set Enrichment Analysis (GSEA), and the components and functions of HT-infiltrating immune cells were examined in these groups. Furthermore, Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) analysis revealed that the proportion of HT-infiltrating immune cells (HTICs) was strongly associated with CD27 expression, suggesting that CD27 has the potential to be a biomarker for the immune status of HT. Consequently, CD27 contributed to predicting clinical outcomes in patients with HT, including disease condition and response to immunotherapy.

Figure 1 Study workflow.

Identification of DEGs Positively Associated with HT Identified Through WGCNA

The expression data from the GSE138198 and GSE54958 datasets were merged after being standardized. Data heterogeneity was assessed using principal component analysis (PCA) (Figure 2A and B). A total of 270 DEGs were identified in patients with HT compared to healthy controls, with a total of 183 upregulated and 87 downregulated genes. A hierarchical clustering heat map and a volcano plot of these DEGs were generated, as shown in Figure 2C and D.

Figure 2 Screening the DEGs highly correlated with HT. (A) Before standardization of training groups. (B) After standardization of training groups. (C and D) Heat map and volcano plot of DEGs. (E-G) GO, KEGG, and DO enrichment analysis of DEGs.

In addition, enrichment analyses of DEGs, including GO, KEGG, and DO enrichments, were conducted. Go analysis, consisting of BPs, CCs, and MFs, was used to identify the top 10 enriched terms for each category (Figure 2E). We observed that DEGs were mainly involved in lymphocyte differentiation, mononuclear cell differentiation, muscle organ development, activation of immune response, muscle cell development, contractile fiber, myofibril, sarcomere, the external side of the plasma membrane, actin binding, carbohydrate binding, cytokine binding, and immune receptor activity.

The KEGG pathways (Figure 2F) included hematopoietic cell lineage, NOD-like receptor signaling pathway, motor proteins, coronavirus disease COVID-19, measles, Epstein-Barr virus infection, cardiac muscle contraction, human T-cell leukemia virus 1 infection, and T cell receptor signaling pathways.

The results of the DO analysis (Figure 2G) included Graves’ disease, hepatitis, nephritis, chronic lymphocytic leukemia, primary immunodeficiency disease, systemic lupus erythematosus, and lupus erythematosus.

Subsequently, weighted gene co-expression network analysis (WGCNA) was used to identify the genes with the highest correlation to HT. All samples were within acceptable quality standards, and hence, all the samples were included in the analysis (Figure 3A). The soft thresholding power (β) was set to 6, with the scale independence reaching 0.8 (Figure 3B).

Figure 3 Screening the module genes highly correlated with HT using WGCNA. (A) Clustering dendrogram of the samples in training groups based on their Euclidean distance. (B) Analysis of network topology for various soft thresholding powers. The x-axis reflects the soft thresholding power. The y-axis reflects the scale-free topology model fit index. The soft thresholding power (β) was set at 6, with the scale Independence reaching 0.8. (C) The interaction analysis of gene modules in normal and HT groups. (D) Clustering dendrogram of genes with dissimilarity, based on the topological overlap, together with assigned module colors. (E) 107 intersection genes between WGCNA and DEGs. A total of 1247 WGCNA genes genes in module_turquoise (1051) and module_salmon (196) were selected to screen the intersection genes with DEGs(269).

The interactions between gene modules in both the normal and HT groups were analyzed (Figure 3C). A cluster dendrogram of genes was constructed based on the dissimilarity of the topological overlaps, with specific colors to distinguish the modules (Figure 3D). A total of 1247 WGCNA genes from module_turquoise (1051) and module_salmon (196) were identified to screen for the intersection genes with DEGs(269). This intersection yielded 107 genes (Figure 3E) that were selected for further exploration.

CD27 Was Found to Be Highly Expressed in HT and Demonstrated Diagnostic Significance

To further investigate the 107 intersection genes identified in the previous analysis and examine the underlying protein levels, a PPI network was constructed using the STRING database (Figure 4A). Next, the network data was analyzed using Cytoscape software to identify and construct cluster gene modules based on the degree of interaction, which reflected the importance of each gene. The top 10 hub genes in the network were identified (Figure 4B), and the top 5 hub genes from among these were used for further analysis. The expression levels of these genes were compared across the training and test groups. Line charts were used to depict the patterns of expression levels of these genes in both groups (Figure 4C).

Figure 4 CD27 is highly expressed in HT and is of great value in the diagnosis of HT. (A) A PPI network of 107 intersection genes constructed using the STRING database. (B) Network of top 10 hub genes constructed using the Cytoscape software. (C) Expression of the top 5 hub genes in the training and test groups. (D) Differential expression analysis of the top 5 hub genes between the normal group and the HT group, and diagnostic ROC curves of the top 5 genes in HT in the training group.**P < 0.01,***P < 0.001. (E) Differential expression analysis of the top 5 hub genes between the normal group and the HT group, and diagnostic ROC curves of the top 5 genes in HT in the test group. **P < 0.01,***P < 0.001. (F) Predicting the accuracy of CD27 in diagnosing HT in training and test groups using nomograms.

In addition, differential expression analysis was conducted on the top 5 hub genes to compare their expression levels between HT patients and healthy controls in both the training and test groups, and diagnostic ROC curves were generated for each of these genes (Figure 4D and E). Excitingly, in the diagnostic ROC curve of CD27, we found that the area under the curve (AUC) of CD27 reached 1 in both the training and validation groups. Based on these findings, CD27 emerged as a key hub gene for further analysis, with significant differential expression and high diagnostic accuracy for HT. Nomograms were used to evaluate the accuracy and predictive value of CD27 in diagnosing HT in the training and test groups (Figure 4F). All these results indicated high expression of CD27 in HT, with significant clinical utility in the diagnosis of HT.

Exploring the Potential Roles of CD27 in the Pathogenesis of HT

We further explored the mechanisms underlying the role of CD27 in the pathogenesis of HT. The comparison of CD27-high and CD-27-low groups resulted in the identification of 154 DEGs. A total of 126 genes were upregulated and 28 genes were downregulated in the CD27-high group compared to the CD27-low group. A hierarchical clustering heat map and a volcano plot of the DEGs were generated (Figure 5A and B). Notably, the DEGs that were highly associated with CD27 included GVINP1, SLAMF6, AIM2, and SMAP2, among others (Figure 5C).

Figure 5 Exploring the potential roles of CD27 in HT. (A and B) Heat map and Volcano plot of DEGs in CD27-high and CD27-low expression groups based on the median expression value of CD27 in the training group. (C) Correlation heat maps of DEGs with high correlation to CD27 expression. (D and E) GO and KEGG enrichment analysis of DEGs in CD27-high and CD27-low expression groups based on the median expression value of CD27 in the training group. (F and G) GSEA enrichment analysis of groups with high CD27 and low CD27 expression.

Further functional enrichment analyses of DEGs were performed using GO, KEGG, and GSEA. GO analysis of the three main categories (BP, CC, and MF) revealed the top 10 enriched terms for each category (Figure 5D). We observed that DEGs related to CD27 were mainly involved in lymphocyte differentiation, mononuclear cell differentiation, B cell differentiation, leukocyte cell-cell adhesion, lymphocyte proliferation, B cell activation, mononuclear cell proliferation, and positive regulation of leukocyte cell-cell adhesion.

The KEGG pathway analysis identified several key pathways in which CD27-associated DEGs were involved (Figure 5E). These included hematopoietic cell lineage, malaria, the T cell receptor signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, primary immunodeficiency, Th17 cell differentiation, cytokine-cytokine receptor interaction, and the NOD-like receptor signaling pathway.

The GSEA results (Figure 5F and G) showed that CD27 expression was positively related to cell adhesion molecule scams, cytokine-cytokine receptor interaction, hematopoietic cell lineage, natural killer cell-mediated cytotoxicity, and the T cell receptor signaling pathway, and negatively correlated with drug metabolism (cytochrome p450), metabolism of xenobiotics (cytochrome p450), glutathione metabolism, the TGF-beta signaling pathway, and the Wnt signaling pathway.

High Expression of CD27 Was Positively Correlated with Immune Cell Infiltration in Patients with HT

To explore the relationship between CD27 and the infiltration of immune cells in the immune microenvironment of HT, the proportions of 22 types of immune cells in samples of both HT and normal tissues were assessed (Figure 6A and B). Significant differences were observed in the proportions of several immune cell types between HT tissues and normal tissues, including B cells naïve, B cells memory, T cells CD4 memory resting, T cells regulatory (Tregs), monocytes, macrophages M0, and macrophages M2 (Figure 6C).

Figure 6 Relationship between CD27 and immune cell infiltration in HT tissues. (A and B) The distribution of various immune cells in normal and HT tissues visualized through histograms and heat maps. (C) Differential analysis of various immune cells between normal and HT tissues. (D) Correlation analysis of immune cells in HT tissues in the training group. (E and F) Correlation analysis between CD27 expression and immune cells. (G) Analysis of the difference of immune cell functions between CD27-high and -low expression groups in the training group.*P < 0.05,**P < 0.01.

Correlation analysis revealed that the positive correlation between T cells CD4 memory resting and NK cells activated was the strongest (r = 0.9), while the strongest negative correlation was between T cells CD8 and CD4 memory resting (r = −0.79) in HT tissues (Figure 6D). Notably, CD27 expression level was strongly positively correlated with T cells, CD4 naïve, and dendritic cells activated, and negatively correlated with monocytes (Figure 6E and F). CD27 was found to positively regulate immune cell functions, such as APC co-stimulation, B cells, cytolytic activity, immature dendritic cells (iDCs), inflammation-promoting, neutrophils, T cell co-inhibition, T cell co-stimulation, T helper cells, T follicular helper cells (Tfh), and tumor-infiltrating lymphocytes (TIL) (Figure 6G).

Construction of a ceRNA Network

Additionally, in this study, a lncRNA-miRNA-CD27 ceRNA network was established to identify potential regulatory mechanisms involving CD27. The ceRNA regulatory network included 13 long non-coding RNA (lncRNA) nodes and 3 microRNA (miRNA) nodes (Figure 7A). CD27 was found to be regulated by has-miR-509-3-5p, hsa-miR-214-3p, and hsa-miR-193a-5p. These miRNAs were, in turn, regulated by various lncRNAs.

Figure 7 CD27 ceRNA network and CD27 expression in clinical samples and thyroid cancer. (A) CD27 ceRNA network based on the differentially expressed lncRNAs, a ceRNA network was constructed: red represents CD27, green represents miRNAs, and blue represents lncRNAs. (B).CD27 mRNA expression levels in clinical normal and HT samples Data are shown as the mean ± SD, ***P < 0.001. (C) CD27 protein levels in clinical normal and HT samples. Data are shown as the mean ± SD, ***P < 0.001. (D and E) CD27 mRNA expression in thyroid cancer and normal tissues of TCGA data. (F and G) Human Protein Atlas analysis of CD27 protein expression levels in thyroid normal and cancer tissues.

Validation of CD27 Expression Using Clinical Samples

Clinical samples, including three normal thyroid samples and five HT samples, were analyzed to further verify the expression levels of CD27. Both mRNA and protein levels of CD27 were higher in HT samples compared to normal thyroid samples (Figure 7B and C).

Expression of CD27 in Thyroid Cancer Using TCGA Data

We used TCGA and HPA databases to explore the expression level of CD27 in thyroid cancer. The mRNA level of CD27 in tumor tissues was lower than that in normal thyroid tissues (Figure 7D and E). The lower CD27 protein expression levels in tumor tissues compared to normal thyroid tissues were further validated by immunohistochemical staining (Figure 7F and G).

Discussion

CD27 has been found to be involved in various types of tumors, such as breast cancer,20–23 lung cancer,24–26 gastric cancer,27–29 colon cancer,30,31 liver cancer,32 and esophageal cancer.32 It is closely associated with the tumor immune microenvironment. Studies have shown that a high expression of CD27 promotes immune cell infiltration and enhances the efficacy of immune therapy.33–35 However, there is limited research on the role of CD27 in non-tumor diseases such as HT. Although the exact etiology of HT has not been fully elucidated, immunoreaction plays a key role in the development of the disease.15 Therefore, further research into the immune effect of HT can improve the understanding of the pathogenesis and treatment of this disease.

In this study, DEGs between the normal and HT groups were compared. GO enrichment analysis showed that DEGs in HT groups were mainly involved in lymphocyte differentiation, mononuclear cell differentiation, and activation of the immune response. KEGG pathway analysis suggested that DEGs in HT groups were strongly associated with the NOD-like receptor signaling pathway, the T cell receptor signaling pathway, and the primary immunodeficiency pathway. In addition, DO enrichment analysis revealed that DEGs in HT groups were related to various immunological diseases, such as primary immunodeficiency disease, human immunodeficiency virus infectious disease, and Graves’ disease. These results suggest a significant association between the pathogenesis of HT and immune dysregulation. Immunotherapy may offer new therapeutic opportunities for the clinical management of HT.36

WGCNA (Weighted Gene Coexpression Network Analysis) is a tool for analyzing gene expression patterns and is widely used in disease and trait related studies. WGCNA uses the strategy of heavy co-expression. This clustering method is more in line with biological phenomena and can reflect the interaction between genes more accurately. Through WGCNA analysis, it is possible to identify hub genes at the center of regulatory networks that may play a key role in biological processes. The WGCNA method is suitable for large sample sizes, and the larger the sample size, the more reliable the results. This gives WGCNA an advantage when processing large-scale gene expression data. WGCNA can use module eigenvalues and hub genes to conduct association analysis with specific traits and phenotypes, so as to analyze biological problems more accurately. However, when the sample size is too low, the calculation of the correlation coefficient may be unreliable, resulting in little value of the obtained regulatory network. Therefore, the WGCNA method requires a sufficient sample size to ensure the reliability of the results. In general, when the number of independent samples is ≥8, the WGCNA co-expression network method based on Pearson correlation coefficient can be considered. When the sample size is ≥15, the WGCNA method has better results. In this study, using WGCNA to identify gene modules that were positively associated with HT, we found that module_turquoise and module_salmon were strongly positively associated with HT, with correlation coefficients of 0.69 and 0.7, respectively. A total of 107 intersection genes between module genes and DEGs were identified for further research. The top 10 hub genes (namely, PTPRC, CD4, CD2, CD27, CD48, CD69, CD3E, CD5, CD3D, and CD53) were screened by constructing a protein-protein interaction network. Comparative analysis between the HT and control groups revealed a markedly higher expression of these genes, particularly CD27, in both the training and test groups of HT samples.

Diagnostic ROC curve analysis revealed that these genes had good diagnostic value for HT in training groups and test groups. The cut-off values of CD27 were 7.171 and 7.012 in the training and test groups, respectively. Notably, the area under the ROC curves (AUC) for CD27 reached 1 in both the training and test groups. Our nomogram analysis indicated that the high expression of CD27 played a considerable role in the diagnosis of HT. The above results show that CD27 is upregulated in HT, indicating that it plays a crucial role in the diagnosis of HT.

To further investigate the role of CD27 in HT and its potential mechanisms, the HT samples in the training group were divided into CD27-high and CD27-low groups based on the median value of CD27. Subsequently, differential expression analysis was conducted to identify DEGs in these two groups. GO results suggested that the DEGs were mainly involved in the proliferation and differentiation of immune cells, such as lymphocyte differentiation, mononuclear cell differentiation, B cell differentiation, lymphocyte proliferation, B cell activation, and mononuclear cell proliferation. The KEGG analysis revealed that the pathways mainly involved included the T cell receptor signaling pathway, the PD-L1 expression and PD-1 checkpoint pathway in cancer, primary immunodeficiency, Th17 cell differentiation, cytokine-cytokine receptor interaction, and the NOD-like receptor signaling pathway. Additionally, GSEA results showed that CD27 expression was positively related to cytokine-cytokine receptor interaction, natural killer cell-mediated cytotoxicity, and the T cell receptor signaling pathway. These results suggest that CD27 plays an important role in the development of HT, mainly by regulating the function of immune cells, especially B cells and T cells.

The production of thyroid peroxidase antibodies in the tissues or serum of patients with HT is causally related to the infiltration of immune cells. In this context, we analyzed the relationship between immune cell infiltration levels in HT groups and CD27 levels. Our results revealed a significant association between CD27 and a variety of immune cells, particularly T cells, with a correlation coefficient of 0.69. These findings suggest that eliminating immune cell infiltration in thyroid tissue may have a beneficial effect in curing HT.

Notably, when we examined the levels of CD27 mRNA and protein in clinical samples, our results showed that CD27 levels were elevated in HT samples, corroborating findings from public databases. In addition, studying the expression levels of CD27 in thyroid cancer revealed that the mRNA level of CD27 in tumor groups was lower than that in normal groups, and there was no difference in the protein level of CD27 in thyroid cancer as compared to the normal group. This further indicates that CD27 could distinguish between HT and thyroid cancer in diagnosis.

Conclusion

In conclusion, the infiltration of immune cells in HT represents a significant aspect of the disease pathology, emphasizing the importance of immune-focused treatment strategies to preserve thyroid gland function in addressing HT. Targeting CD27 as a therapeutic intervention for HT may be beneficial, and this approach should not be limited only to cases where thyroid function is abnormal but should be considered for all patients with HT. However, given the limited number of HT cases in this study, more in-depth research is warranted to validate these findings.

Abbreviations

HT, Hashimoto’s thyroiditis; DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis; PPI, protein–protein interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; DO, Disease Ontology; ssGSEA, Single-sample GSEA; TPOAb, anti-thyroperoxidase antibody; TGAb, anti-thyroglobulin antibody; BP, biological processes; CC, cell component; MF, molecular function.

Data Sharing Statement

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Ethics Approval and Consent to Participate

The study complies with the Declaration of Helsinki.

Acknowledgments

This paper has been uploaded to ResearchSquare as a preprint: https://www.researchsquare.com/article/rs-3389491/v1.

Funding

National Natural Science Foundation of China (No. 81572916).

Disclosure

The authors declare that they have no competing interests.

References

1. Caturegli P, De Remigis A, Rose NR. Hashimoto thyroiditis: clinical and diagnostic criteria. Autoimmun Rev. 2014;13(4–5):391–397. doi:10.1016/j.autrev.2014.01.007

2. Antonelli A, Ferrari SM, Corrado A, Di Domenicantonio A, Fallahi P. Autoimmune thyroid disorders. Autoimmun Rev. 2015;14(2):174–180. doi:10.1016/j.autrev.2014.10.016

3. Hollowell JG, Staehling NW, Flanders WD, et al. Serum TSH, T(4), and thyroid antibodies in the United States population (1988 to 1994): national health and nutrition examination survey (NHANES III). J Clin Endocrinol Metab. 2002;87(2):489–499. doi:10.1210/jcem.87.2.8182

4. McLeod DS, Cooper DS. The incidence and prevalence of thyroid autoimmunity. Endocrine. 2012;42(2):252–265. doi:10.1007/s12020-012-9703-2

5. Ragusa F, Fallahi P, Elia G, et al. Hashimotos’ thyroiditis: epidemiology, pathogenesis, clinic and therapy. Best Pract Res Clin Endocrinol Metab. 2019;33(6):101367. doi:10.1016/j.beem.2019.101367

6. Pyzik A, Grywalska E, Matyjaszek-Matuszek B, Roliński J. Immune disorders in Hashimoto’s thyroiditis: what do we know so far? J Immunol Res. 2015;2015:979167. doi:10.1155/2015/979167

7. Pishdad P, Pishdad GR, Tavanaa S, Pishdad R, Jalli R. Thyroid ultrasonography in differentiation between graves’ disease and Hashimoto’s thyroiditis. J Biomed Phys Eng. 2017;7(1):21–26. doi:10.1245/s10434-014-3556-2

8. Saraf SR, Gadgil NM, Yadav S, Kalgutkar AD. Importance of combined approach of investigations for detection of asymptomatic Hashimoto Thyroiditis in early stage. J Lab Physicians. 2018;10(3):294–298. doi:10.4103/JLP.JLP_72_17

9. Gonzalez-Aguilera B, Betea D, Lutteri L, et al. Conversion to graves disease from Hashimoto thyroiditis: a study of 24 patients. Arch Endocrinol Metab. 2018;62(6):609–614. doi:10.20945/2359-3997000000086

10. Lee MJ, Kim BY, Ma JS, et al. Hashimoto thyroiditis with an unusual presentation of cardiac tamponade in Noonan syndrome. Korean J Pediatr. 2016;59(Suppl 1):S112–s5. doi:10.3345/kjp.2016.59.11.S112

11. Özon A, Tekin N, Şıklar Z, et al. Neonatal effects of thyroid diseases in pregnancy and approach to the infant with increased TSH: Turkish Neonatal and Pediatric Endocrinology and Diabetes Societies consensus report. Turk Pediatri Ars. 2018;53(Suppl 1):S209–s23. doi:10.5152/TurkPediatriArs.2018.01819

12. Liang J, Zeng W, Fang F, et al. Clinical analysis of Hashimoto thyroiditis coexistent with papillary thyroid cancer in 1392 patients. Acta Otorhinolaryngol Ital. 2017;37(5):393–400. doi:10.14639/0392-100X-1709

13. Silva de Morais N, Stuart J, Guan H, et al. The impact of Hashimoto thyroiditis on thyroid nodule cytology and risk of thyroid cancer. J Endocr Soc. 2019;3(4):791–800. doi:10.1210/js.2018-00427

14. Feldt-Rasmussen U. Hashimoto’s thyroiditis as a risk factor for thyroid cancer. Curr Opin Endocrinol Diabetes Obes. 2020;27(5):364–371. doi:10.1097/MED.0000000000000570

15. Ralli M, Angeletti D, Fiore M, et al. Hashimoto’s thyroiditis: an update on pathogenic mechanisms, diagnostic protocols, therapeutic strategies, and potential malignant transformation. Autoimmun Rev. 2020;19(10):102649. doi:10.1016/j.autrev.2020.102649

16. Klubo-Gwiezdzinska J, Wartofsky L. Hashimoto thyroiditis: an evidence-based guide to etiology, diagnosis and treatment. Pol Arch Intern Med. 2022;132(3). doi:10.20452/pamw.16222

17. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–883. doi:10.1093/bioinformatics/bts034

18. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. doi:10.1093/nar/gkv007

19. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9(1):559. doi:10.1186/1471-2105-9-559

20. Fang J, Chen F, Liu D, Gu F, Chen Z, Wang Y. Prognostic value of immune checkpoint molecules in breast cancer. Biosci Rep. 2020;40(7). doi:10.1042/BSR20201054

21. Hader M, Savcigil DP, Rosin A, et al. Differences of the immune phenotype of breast cancer cells after ex vivo hyperthermia by warm-water or microwave radiation in a closed-loop system alone or in combination with radiotherapy. Cancers. 2020;12(5):1082. doi:10.3390/cancers12051082

22. Bertucci F, Boudin L, Finetti P, et al. Immune landscape of inflammatory breast cancer suggests vulnerability to immune checkpoint inhibitors. Oncoimmunology. 2021;10(1):1929724. doi:10.1080/2162402X.2021.1929724

23. Rapoport BL, Steel HC, Benn CA, et al. Dysregulation of systemic soluble immune checkpoints in early breast cancer is attenuated following administration of neoadjuvant chemotherapy and is associated with recovery of CD27, CD28, CD40, CD80, ICOS and GITR and substantially increased levels of PD-L1, LAG-3 and TIM-3. Front Oncol. 2023;13:1097309. doi:10.3389/fonc.2023.1097309

24. Kamphorst AO, Pillai RN, Yang S, et al. Proliferation of PD-1+ CD8 T cells in peripheral blood after PD-1-targeted therapy in lung cancer patients. Proc Natl Acad Sci. 2017;114(19):4993–4998. doi:10.1073/pnas.1705327114

25. Larroquette M, Guegan JP, Besse B, et al. Spatial transcriptomics of macrophage infiltration in non-small cell lung cancer reveals determinants of sensitivity and resistance to anti-PD1/PD-L1 antibodies. J Immunother Cancer. 2022;10(5):e003890. doi:10.1136/jitc-2021-003890

26. Wang Q, He Y, Li W, et al. Soluble immune checkpoint-related proteins in blood are associated with invasion and progression in non-small cell lung cancer. Front Immunol. 2022;13:887916. doi:10.3389/fimmu.2022.887916

27. Sukri A, Hanafiah A, Kosai NR, Mohamed Taher M, Mohamed Rose I. Surface antigen profiling of helicobacter pylori-infected and -uninfected gastric cancer cells using antibody microarray. Helicobacter. 2016;21(5):417–427. doi:10.1111/hel.12295

28. Yamakoshi Y, Tanaka H, Sakimura C, et al. Immunological potential of tertiary lymphoid structures surrounding the primary tumor in gastric cancer. Int J Oncol. 2020;57(1):171–182. doi:10.3892/ijo.2020.5042

29. Yang Y, He W, Wang ZR, et al. Immune cell landscape in gastric cancer. Biomed Res Int. 2021;2021:1930706. doi:10.1155/2021/1930706

30. Ding J, Wang H, Hou R, et al. Total T cell density and expression of t memory stem cell markers are associated with better prognosis in colon cancer. Int J Gen Med. 2023;16:2285–2294. doi:10.2147/IJGM.S411122

31. Liu H, Li L, Hao Y, et al. Identification of two migratory colon ILC2 populations differentially expressing IL-17A and IL-5/IL-13. Sci China Life Sci. 2023;66(1):67–80. doi:10.1007/s11427-022-2127-2

32. Hu B, Yu M, Ma X, et al. IFNα potentiates anti-PD-1 efficacy by remodeling glucose metabolism in the hepatocellular carcinoma microenvironment. Cancer Discov. 2022;12(7):1718–1741. doi:10.1158/2159-8290.CD-21-1022

33. Wajant H. Therapeutic targeting of CD70 and CD27. Expert Opin Ther Targets. 2016;20(8):959–973. doi:10.1517/14728222.2016.1158812

34. Crinier A, Narni-Mancinelli E, Ugolini S, Vivier E. SnapShot: natural killer cells. Cell. 2020;180(6):1280–e1. doi:10.1016/j.cell.2020.02.029

35. Starzer AM, Berghoff AS. New emerging targets in cancer immunotherapy: CD27 (TNFRSF7). ESMO Open. 2020;4(Suppl 3):e000629. doi:10.1136/esmoopen-2019-000629

36. Hu Y, Feng W, Chen H, et al. Effect of selenium on thyroid autoimmunity and regulatory T cells in patients with Hashimoto’s thyroiditis: a prospective randomized-controlled trial. Clin Transl Sci. 2021;14(4):1390–1402. doi:10.1111/cts.12993