Identification of potential immunotherapeutic targets and prognostic biomarkers in Graves' disease using weighted gene co-expression network analysis

Abstract

Graves' disease (GD) is an autoimmune disorder characterized by hyperthyroidism resulting from autoantibody-induced stimulation of the thyroid gland. Despite recent advancements in understanding GD's pathogenesis, the molecular processes driving disease progression and treatment response remain poorly understood. In this study, we aimed to identify crucial immunogenic factors associated with GD prognosis and immunotherapeutic response. To achieve this, we implemented a comprehensive screening strategy that combined computational immunogenicity-potential scoring with multi-parametric cluster analysis to assess the immunomodulatory genes in GD-related subtypes involving stromal and immune cells. Utilizing weighted gene co-expression network analysis (WGCNA), we identified co-expressed gene modules linked to cellular senescence and immune infiltration in CD4⁺ and CD8⁺ GD samples. Additionally, gene set enrichment analysis enabled the identification of hallmark pathways distinguishing high- and low-immune subtypes. Our WGCNA analysis revealed 21 gene co-expression modules comprising 1,541 genes associated with immune infiltration components in various stages of GD, including T cells, M1 and M2 macrophages, NK cells, and Tregs. These genes primarily participated in T cell proliferation through purinergic signaling pathways, particularly neuroactive ligand-receptor interactions, and DNA binding transcription factor activity. Three genes, namely PRSS1, HCRTR1, and P2RY4, exhibited robustness in GD patients across multiple stages and were involved in immune cell infiltration during the late stage of GD (p < 0.05). Importantly, HCRTR1 and P2RY4 emerged as potential prognostic signatures for predicting overall survival in high-immunocore GD patients (p < 0.05). Overall, our study provides novel insights into the molecular mechanisms driving GD progression and highlights potential key immunogens for further investigation. These findings underscore the significance of immune infiltration-related cellular senescence in GD therapy and present promising targets for the development of new immunotherapeutic strategies.

Keywords: Gene co-expression modules; Graves' disease; Immune infiltration; Immunotherapeutic biomarkers; Molecular mechanisms; WGCNA.

Fig. 1

Flowchart describing the schematic overview of the study design for investigating immune infiltration characteristics in patients with GD.

Fig. 2

Differences in immune infiltration characteristics among patients with different stages of GD (A) and different ages (B). (C) The nine immune cell scores differed among the three stages based on IPS, MCP, and Xcell algorithms. (D) Immune cell fractions between different stages of GD. Values are mean ± SEM, *p < 0.05, **p < 0.001.

Fig. 3

Immune screening of gene expression in CD4⁺ and CD8⁺ GD samples from the TCGA database was evaluated based on the IPS (A) and MCP (B) algorithms. The vertical gray line denotes a subset of high immune score expression outliers (right) determined by the 1.5 × interquartile range. Box plots (far right) display aggregate expression of the genes in high and low immune score groups. Statistical analysis was performed using the Mann-Whitney U test.

Fig. 4

Weighted gene correlation network analysis (WGCNA) of CD4⁺ and CD8⁺ GD samples based on stage (mild, moderate, and severe) and age range. (A) Clustering dendrogram and module-trait analysis. Different colors of the column indicate different hub modules. (B) Analysis of the scale-free fit index and mean connectivity for various soft-threshold powers. (C) Clustering dendrogram of all genes with dissimilarity based on the topological overlap and assigned module colors. (D) Heatmap for the correlation between immune modules and traits. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5

Functional analyses of the selected immune module. (A) Pathway enrichment and main pathway networks of top-selected immune module dark-red module. (B) The PPI bipartite network analysis for the dark-red module. (C) Cognitive impairment analysis of the selected immune module. (D) Integration of the selected biological processes of purinergic signaling pathways, neuroactive ligand-receptor interaction, and DNA binding transcription factor activity. (E) Signaling enriched KEGG pathway analysis of purinergic signaling pathways. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6

Diagnostic accuracy analysis of targeted mRNA. (A) Boxplots of the targeted potential key mRNA drivers T cell proliferation, PRSS1, HCRTR1, and P2RY4 between high- and low-immune score samples (A), between CD4⁺ and CD8⁺ T cell samples (B), and different stages of GD samples (C) groups. These three genes were overexpressed in the metastatic groups. Horizontal bars represent mean ± SD. (D) Receiver operating characteristic (ROC) curve for PRSS1, HCRTR1, and P2RY4 expression. (E) Correlation of PRSS1, HCRTR1, and P2RY4 with other clinical trials of the samples.

Fig. S1

Immune cell subsets distribution by three algorithms - IPS, MCP, and Xcell.

Fig. S2

Significance of Biological Pathways in Purinergic Signaling in GD. This Figure shows the significance of the most prominent biological pathways related to neuroactive ligand-receptor interactions and DNA binding transcription factor activity in the purinergic signaling pathways, which have been linked to the regulation of immune response and inflammation in GD.