Exploring Prognostic Gene Expression Signature for Thyroid Eye Disease

1. Materials and Methods

1.1. Dataset and Data Processing

The TED datasets GSE58331[15,16,17,18,19,20] (35 TED samples and 29 Control samples), GSE105149[21] (4 TED samples and 7 Control samples) were obtained by downloading them from the GEO database[22] (https://www.ncbi.nlm.nih.gov/geo/), separately. The GSE105149 dataset analyzes gene expression in lacrimal gland biospecimens and surgical waste tissues, while GSE58331 focuses on gene expression in orbital biopsies. Both datasets utilized the GPL570 microarray platform, ensuring technical consistency. MEMRGs were sourced from the GeneCards database[23] (https://www.genecards.org/). The set of MEMRGs from the published literature[24] was obtained from the PubMed website (https://pubmed.ncbi.nlm.nih.gov/), and after merging and de-duplication, a total of 386 MEMRGs. Datasets GSE58331 and GSE105149 were de-batched using R package sva[25], which were normalised by R package limma[26], annotated probes were standardized and normalized. Distribution boxplots indicate successful removal of batch effects, reducing the impact of sample type variations on results. All chosen samples from GSE58331 were part of the test set, while those from GSE105149 were the validation set for further analyses.

1.2. Analysis of Differentially Expressed Genes

Genes in the TED and Control groups were differentially analyzed using the R package limma. A threshold of |logFC| > 0.30 and p-value < 0.05 was applied for identifying Differentially Expressed Genes (DEGs). Among these, genes with logFC > 0.30 and p-value < 0.05 were classified as Up-regulated Genes, while genes with logFC < -0.30 and p-value < 0.05 were categorized as Down-regulated Genes. Mitochondrial Energy Metabolism-Related Differentially Expressed Genes (MEMRDEGs) were identified by intersecting genes from all TED samples in the DEGs dataset with MEMRGs.The cutoff (|Log2FC| > 0.3) chosen for the DEG analysis was carefully considered. We aimed to capture genes potentially linked to TED's pathological processes, ensuring inclusion of candidates with biological relevance. Additionally, similar or more permissive Log2FC thresholds have been used in previous studies (e.g., PMID: 39021997), supporting our decision.

1.3. Functional Enrichment Analysis of Differentially Expressed MMRGs

For Gene Ontology (GO)[27] and KEGG[28] enrichment analysis, the “clusterProfiler”[29]software package was used to explore the enrichment degree of biological process (BP), cellular component (CC), molecular function (MF) , and biological pathway (KEGG) of the differentially expressed MMRGs. P value < 0.05 and FDR value (q value) < 0.25 was considered statistically significant. Data visualization was performed.

1.4. GSEA

Gene Set Enrichment Analysis (GSEA)[30]was conducted on all genes in dataset GSE58331 using the "clusterProfiler" package in R, sorted based on logFC values. TED samples from GSE58331 were then divided into high-risk and low-risk groups using the median RiskScore value. Genes meeting the criteria of |logFC| > 0.30 and p-value < 0.05 in these groups were selected for variance analysis with the limma package. Subsequently, GSEA was repeated on all genes within the different risk groups using the "clusterProfiler" package..

The analysis utilized the seed 2020, 1000 permutations, a minimum of 10 genes per set, and a maximum of 500 genes per set. The gene set c2.cp.all.v2022.1.Hs.symbols.gmt [All Canonical Pathways] (3050 sets) was sourced from the Molecular Signatures Database (MSigDB)[31] for GSEA, and the gene set GSEA screening criteria were adj. p < 0.05 and FDR value (q value) < 0.25, and the p value correction method was Benjamini-Hochberg (BH).

1.5. Screening of Key Genes

To identify Key Genes and build diagnostic models, we initially conducted one-way logistic regression on MEMRDEGs, considering a significance threshold of p < 0.05 for screening. The selected MEMRDEGs were then included in the subsequent random forest analysis.Model construction was performed using the Random Forest package[32] with the expression of the MEMRDEGs obtained from the screening of the one-way logistic regression in the expression matrix of the dataset GSE58331 with the parameters set.seed (234), ntree = 300. MeanDecreaseGini is the average decrease in the Gini coefficient. the Gini coefficient represents the impurity of a node, the larger the Gini coefficient the less pure it is and the more impurities it contains. We then trade-off the number of variables by performing five ten-fold cross-validations in conjunction with the cross-validation curves. We used the training set itself for cross-validation, retaining a relatively small number of variables with relatively small errors, and then selected our significant variables for subsequent analyses based on MeanDecreaseGini. Least Absolute Shrinkage and Selection Operator(LASSO)[33] Regression analyses were conducted following random forest screening using the "glmnet" R package, with set.seed(300) to set the seed and a run period of 200 to prevent overfitting. The outcomes of the LASSO regression analyses were visualized using diagnostic model plots and variable trajectory plots. MEMRDEGs incorporated in the ultimate LASSO regression model were identified as the Key Genes for further analyses.

We selected LASSO regression for its robust feature selection, effectively managing high-dimensional data while minimizing overfitting. Its penalty term inherently eliminates less relevant variables, making it ideal for identifying key genes. Compared to elastic net, LASSO offers greater simplicity and interpretability, which is crucial when prioritizing model clarity.

1.6. Key Genes Constructs Diagnostic Logistic Regression Models

We developed a multifactorial logistic model to derive coefficients for each Key Gene. These coefficients were multiplied by their respective expressions to calculate the Risk Score for each sample. Based on the median Risk Score, the dataset's disease groups were categorized into high and low-risk groups (High, Low). The Risk Score was computed using the following formula:

$$\mathrm{Risk}Score=\sum _{i}Coefficient\left({gene}_i\right)\mathrm{*}mRNAExpression\left({gene}_i\right)$$

Nomogram[35] is a graphical representation of the functional relationships among multiple independent variables in a two-dimensional rectangular coordinate system, typically showing disjointed line segments. The results of the multifactorial logistic model were visualized using the "rms" R package. Plotting the columns illustrated the interconnections among the Key Genes featured in the multifactorial logistic model.

Calibration Curve to evaluate the accuracy and discriminatory power of the model based on the results of the multifactorial logistic regression analysis through Calibration analysis. Decision Curve Analysis (DCA) was plotted through Decision Curve Analysis (DCA)[36] using the R package “ggDCA” based on Key Genes from the dataset GSE58331.

Finally, the R package pRO was used to plot the ROC curves of the Logistic Risk Score (Risk Score) in datasets GSE58331 and GSE105149 and to calculate the Area Under the Curve (AUC) to assess the expression of Logistic Risk Score (Risk Score) volume on the diagnostic effect of TED occurrence. The functional correlation of Key Genes was calculated by R package “GOSemSim”[37], and the functional correlation between Key Genes was analysed for Friends. Finally, the chromosomal localisation map was drawn to show the position of Key Genes on chromosomes by R package “RCircos”[38].

1.7. Validation of Differential Expression of Key Genes

To further explore the expression differences of Key Genes in the datasets GSE58331 and GSE105149 in the TED group and on the Control group, group comparison plots were drawn based on the expression of Key Genes. Finally, the R package “pROC” was used to plot the ROC curve of Key Genes and calculate the Area Under the Curve (AUC) to assess the diagnostic effect of Key Genes expression on the occurrence of TED The AUC of the ROC curve is generally between 0.5 and 1. The closer the AUC is to 1, the better the diagnostic effect. The AUC of 0.5 to 0.7 has low accuracy, the AUC of 0.7 to 0.9 has some accuracy, and the AUC of 0.9 or more has high accuracy.

1.8. Key Genes' ssGSEA Algorithm for Immune Infiltration Analysis

Single-Sample Gene-Set Enrichment Analysis (ssGSEA)[39] quantifies the relative abundance of each immune cell infiltration. First, each infiltrating immune cell type is labeled, such as Activated CD8 T cell, Activated dendritic cell, Gamma delta T cell, Natural killer cell, Regulatory T cell, and many other human immune cell subtypes. Next, the enrichment scores calculated using ssGSEA analysis were used to represent the relative abundance of each immune cell infiltration in each sample to derive the immune cell infiltration matrix. Subsequently, immune cells with significant differences in the two groups were screened for subsequent analysis, the correlation between immune cells was calculated based on the Spearman algorithm, and the correlation heatmap was drawn using the R package pheatmap to demonstrate the results of the correlation analysis of the immune cells themselves. The correlation between Key Genes and immune cells was calculated based on the Spearman algorithm, and the result of p value < 0.05 was retained, and the correlation bubble map was plotted using the R package “ggplot2” to show the results of the correlation analysis between Key Genes and immune cells.

1.9. mRNA-miRNA, mRNA-TF, mRNA-RBP Interaction Network Analysis of Key Genes

ENCORI Database[40] is the version 3.0 of starBase database. miRNA-ncRNA, miRNA-mRNA, ncRNA-RNA, RNA-RNA, RBP-ncRNA and RBP-mRNA interactions in ENCORI database are based on the data mining of CLIP-seq and degradome sequencing (for plants), and we utilized ENCORI database to predict miRNAs interacting with Key Genes and screen mRNA-miRNA interaction pairs using pancancerNum > 7 as a screening criterion, and visualize the mRNA-miRNA interaction network using Cytoscape software.

CHIPBase database (version 3.0)[41] (https://rna.sysu.edu.cn/chipbase/) identifies thousands of binding motif matrices and their binding sites from ChIP-seq data of DNA-binding proteins and predicts transcriptional regulatory relationships between millions of transcription factors (TFs) and genes. We searched for TFs binding to Key Genes through the CHIPBase database (version 3.0) and used the sum of Number of samples found (upstream) and Number of samples found (downstream) greater than 10 as a filtering criterion for mRNA-TF interaction pairs. The mRNA-TF interaction pairs were screened and the mRNA-TF interaction network was visualized using Cytoscape software.

RNA-Binding Protein (RBP)[42] plays a key role in gene regulation, and is critical for the regulation of life activities such as RNA synthesis, selective splicing, modification, translocation and translation. Based on the StarBase v3.0 database (https://starbase.sysu.edu.cn/), the target RBPs of Key Genes were predicted and the mRNA-RBP interaction pairs were screened using clusterNum > 12 as the screening criterion. mRNA-RBP interaction networks were visualized by Cytoscape software ( mRNA-RBP Interaction Network).

1.10. Statistical Analysis

All data processing and analyses in this study were conducted using R software (Version 4.2.2). For comparing two groups of continuous variables, statistical significance of normally distributed variables was assessed using the independent Student's T-Test, while non-normally distributed variables were analyzed with the Mann-Whitney U-Test (Wilcoxon Rank Sum Test). Comparisons involving three or more groups utilized the Kruskal-Wallis test. Spearman correlation analysis was employed to calculate correlation coefficients between different molecules. Unless specified otherwise, all p values were two-sided, with p < 0.05 considered statistically significant.

2. Results

2.1. MEMRGs Landscape in TED

After using R package “sva” to remove the batch effect on TED dataset GSE58331, GSE105149, the difference in the expression value of the dataset before and after the removal of the batch effect was compared by distribution box line plots (Fig.1A-B), (Fig.1C-D) , The results showed that the batch effect of the samples in TED dataset is basically eliminated after removal of the batch treatment.

2.2. TED-related Differentially Expressed MMRG-Mediated Signaling Pathways

The GSE58331 dataset had a total of 3994 DEGs; up-expressed genes had a total of 1972, and down-expressed genes had a total of 2022, and the GSE105149 dataset had a total of 3952 DEGs; the up-expressed genes had a total of 2027, and the down-expressed genes had a total of 1925, and the results of the differential analysis of the dataset (Fig.2A, Fig.2B).

The intersection (Fig.2C) was taken for all TED samples of genes and MEMRGs in all GSE58331-DEGs, GSE105149-DEGs obtained. A total of 26 MEMRDEGs were obtained, namely UQCRFS1, ALDH1A1, GPI, IMMT, GCK, NDUFS3, IDH1, ETFA, UQCRC2, PGK1, MFF, COX6C, BOLA3, SOAT2, MDH1, C1QBP, ALDH3A1, NDUFV2, NDUFS1, AKR1A1, ALDH3A2, DHTKD1, FXN, ACADS, SLC25A47, EPO. Based on the intersection results, the expression differences of MEMRDEGs between different sample subgroups in GSE58331 and GSE105149 datasets were analyzed separately and heatmaps were plotted to show the results of the analysis using the R package pheatmap (Fig.2D-E).

BOLA3 displays contrasting expression patterns in the GSE58331 and GSE105149 datasets, possibly due to various factors. Variations in sample origins and processing methods, with GSE58331 utilizing anterior orbital tissues and GSE105149 using lacrimal glands, may impact gene expression levels. The smaller sample size in GSE105149 could result in greater result instability and variability. Furthermore, clinical sample heterogeneity, including diverse biological backgrounds and disease stages among patients, may contribute to expression disparities of the same gene across datasets.

2.3. Analysis of GO and KEGG

A total of 26 MEMRDEGs were subjected to GO and KEGG enrichment analyses as shown in Table 1. The results showed that 26 MEMRDEGs were predominantly enriched in TED. The results of the analyses were visualized by bubble plots (Fig.3A). Meanwhile, the network diagrams of biological process (BP), cellular component (CC), molecular function (MF) and biological pathway (KEGG) were plotted based on the GO and KEGG enrichment analysis (Fig.3B-E). The connecting lines illustrated the corresponding molecules and the annotations of the corresponding entries; the larger the node, the more molecules the entry contained.

2.4. GSEA in TED

The association between the expression of all genes in the dataset GSE58331 and the biological processes involved, the cellular components affected, and the molecular functions performed was investigated by GSEA (Fig.4 A), and the specific results were shown in Table 2. The results showed that all genes in the dataset GSE58331 were significantly enriched in Fatty Acids (Fig.4 B), the Overview of Proinflammatory and Profibrotic Mediators (Fig.4 C), Dilated Cardiomyopathy (Fig.4 D), ADORA2B Mediated Anti Inflammatory Cytokines Production (Fig.4 E) and other biologically relevant functions and signaling pathways.

2.5. Screening of Key Genes

To determine the diagnostic value of MEMRDEGs in TED, 26 MEMRDEGs were screened for significance using one-way logistic regression (p value<0.05). The RF algorithm was used to analyze the expression of these 26 MEMRDEGs in the TED/Control disease subgroup of the dataset GSE58331. Setting the seed to 234 and the number of decision trees to 300, we plotted the decision tree error curve (Fig.5A), which showed that the error reached its minimum and leveled off when the number of decision trees was around 50. We then plotted the MeanDecreaseGini scatterplot (Fig.5B) for important gene screening. The number of genes was rounded off by performing five ten-fold cross-validation and plotting the cross-validation error plot (Fig.5C). The graph showed that the model error was minimized when the number of genes was 12, which was then combined with MeanDecreaseGini to select specific genes for subsequent analysis. The results showed (Fig.5B-C) that the 12 MEMRDEGs screened by this algorithm that had a significant impact on the diagnosis of TED were, in descending order of importance, GPI, COX6C, NDUFS3, MDH1, ALDH1A1, IDH1, FXN, IMMT, MFF, GCK, UQCRFS1, NDUFV2.

LASSO risk models were constructed by LASSO regression analysis for these 12 MEMRDEGs. The LASSO regression model was plotted (Fig.5D) and the trajectory of LASSO variables (Fig.5E). The results showed that the LASSO risk model contained a total of 10 MEMRDEGs, which were IDH1, IMMT, MDH1, NDUFS3, NDUFV2, GPI, UQCRFS1, GCK, ALDH1A1, and MFF. These genes were taken as Key Genes for our subsequent study, and the forest plot of Key Genes was drawn ( Fig.5F).

2.6. Key Genes Constructing Diagnostic Logistic Regression Models

The Risk Score formula for the dataset GSE58331 was as follows:

$$\begin{array}{cc}\mathrm{Risk}Score=\hfill & -1.9687\mathrm{*}IDH1-2.0103\mathrm{*}IMMT-0.6120\mathrm{*}MDH1-2.3690\mathrm{*}NDUFS3+4.1636\hfill \\ & \mathrm{*}NDUFV2-2.2142\mathrm{*}GPI-3.5576\mathrm{*}UQCRFS1+3.3602\mathrm{*}GCK-0.5756\hfill \\ & \mathrm{*}ALDH1A1+4.3804\mathrm{*}MFF\hfill \end{array}$$

The TED patients were classified into high-risk and low-risk groups using the risk score. Then we plotted the column line graph (Fig. 6A) to show the connection between the 10 Key Genes, and we could see that the expression of MFF contributes the most to the multifactorial logistic model.

In order to judge the accuracy and discriminative power of the multifactor logistic model, the Calibration Curve was plotted by Calibration analysis, and the predictive effect of the model on the actual results was assessed based on the fit between the actual probabilities and the probabilities predicted by the model in different cases plotted in the graph (Fig.6B). The calibration curve plot of the multifactor logistic model showed that the calibration line shown by the dashed line basically coincided with the diagonal line of the ideal model, which indicated that the model had high accuracy and discriminatory power.

The diagnostic effectiveness of the diagnostic multifactorial logistic model for TED diagnosis was evaluated by decision curve analysis (DCA) based on the dataset GSE58331(Fig.6C). The results showed that the line of the model was stably higher than All and None in a certain range, and the net gain of the model was higher, indicating that the model was more effective in diagnosis.

Finally, ROC curves were plotted using the R package “pROC” based on the Risk Score in the datasets GSE58331 and GSE105149 to validate the diagnostic utility of TED with the diagnostic multifactorial logistic model. showing that the multifactorial logistic model presented a high accuracy for the diagnosis of TED (Fig.6D-E. AUC = 0.910, AUC = 1.000).

The high or low score of Friends analysis was used to determine the genes that TED played an important role in the biological process (Fig.6F). The results showed that NDUFS3 played an important role in TED as the gene closest to the cut-off value (cut-off value = 0.60).

The positions of 10 Key Genes on human chromosomes were analyzed by R package “RCircos” (Fig. 6G). The chromosome localization map showed that more Key Genes were located on chromosome 2, namely, MDH1, IMMT, IDH1, and MFF were located on chromosome 2(Fig. 6G).

2.7. Expression Differential Validation Analysis and Functional Similarity Analysis of Key Genes

In order to explore the expression differences of 10 Key Genes, the results of the difference analysis of in dataset GSE58331 in TED group and Control group with high and low expression were demonstrated by group comparison plot (Fig. 7A). The difference results showed (Fig. 7A) that six Key Genes had highly statistically significant (p value < 0.001) expression in the TED group and Control group of dataset GSE58331, namely, IMMT, NDUFS3, GPI, UQCRFS1, GCK, and ALDH1A1. Finally, using the R package “pROC” to draw ROC curves based on the expression of Key Genes in the dataset GSE58331. The ROC curves (Fig. 7B-F) showed that the expression levels of all Key Genes were very similar to those in the classification of TED group and Control group presented some accuracy (0.7 < AUC < 0.9).

Using the same methodology, the analysis of differences in expression levels of 10 Key Genes between the TED and Control groups in dataset GSE105149 revealed significant findings. A group comparison plot (Fig. 7G) displayed that IDH1, IMMT, and NDUFS3 exhibited highly significant expression variances (p < 0.01) between the TED and Control groups. Subsequently, the R package "pROC" was utilized to generate ROC curves based on dataset GSE105149. The ROC curves (Fig. 7H-L) indicated that IDH1, IMMT, MDH1, NDUFS3, NDUFV2, ALDH1A1, and MFF expression levels among the Key Genes were highly accurate in distinguishing the TED and Control groups (AUC > 0.9). Additionally, GPI, UQCRFS1, and GCK expression levels showed moderate accuracy in classifying the TED and Control groups (0.7 < AUC < 0.9).

2.8. GSEA Based on Logistic Risk Score High/Low Grouping

In order to determine the effect of the expression levels of all genes in dataset GSE58331 on TED between high and low risk score subgroups, the association between the expression of all 21655 genes in dataset GSE58331 of TED between different high and low risk score subgroups (High/Low) and the biological processes involved, cellular components affected, and molecular functions performed were investigated by GSEA . (Fig.8A), the specific results were shown in Table 3. The results showed that all genes in the dataset GSE58331 were significantly enriched in Activation of the Tfap2 Ap 2 Family of Transcription Factors (Fig.8B), Negative Regulation of Activity of the Tfap2 Ap 2 Family of Transcription Factors (Fig.8C), Signaling By Interleukins (Fig.8D), IL18 Signaling Pathway (Fig.8E), and other biologically relevant functions and signaling pathways.

2.9. ssGSEA Algorithm Based on Logistic Risk Score High/Low Grouping for Immune Infiltration Analysis

The expression matrix of TED samples from dataset GSE58331 was applied to calculate the immune infiltration abundance of 28 immune cells in TED the high and low risk groups by the ssGSEA algorithm. First, the differences in the expression of immune cell infiltration abundance in different subgroups were demonstrated by group comparison plots (Fig.9A) , showing that activated CD8 T cell, Immature B cell and T follicular helper cell were statistically significant (p value < 0.05).Finally, the correlation of Key Genes with the abundance of immune cell infiltration was demonstrated by correlation bubble plots (Fig.9B-C). The results of correlation bubble plots showed that TED Key Genes showed a strong positive correlation with immune cells in the Low Risk (LR) and High Risk (HR) groups. Among them, UQCRFS1 showed the strongest positive correlation with T follicular helper cells in the Low Risk group (r value = 0.529, p value < 0.05); MDH1 showed the strongest positive correlation with T follicular helper cells in the High Risk group (r value = 0.733, p value < 0.05).

2.10. Interaction Networks Analysis of mRNA-miRNA, mRNA-TF, and mRNA-RBP at Key Genes

Firstly, miRNAs associated with 10 Key Genes were obtained from ENCORI database, and the mRNA-miRNA Interaction Network was constructed and visualized using Cytoscape software (Fig.10A). Among them, 6 Key Genes (ALDH1A1, GPI, IDH1, MFF, NDUFS3, UQCRFS1) and 44 miRNAs were included.

Next, the TFs combined with Key Genes were obtained from ChIPBase database, and the mRNA-TF Interaction Network was constructed and visualized using Cytoscape software (Fig.10B). Among them, 8 Key Genes (IDH1, ALDH1A1, GCK, IMMT, MDH1, MFF, NDUFV2, UQCRFS1) and 33 TFs were included.

Finally, RBPs associated with Key Genes were predicted by the StarBase database, and the mRNA-RBP Interaction Network was constructed and visualized using Cytoscape software (Fig.10C). Among them, 8 Key Genes (ALDH1A1, GPI, IDH1, IMMT, MDH1, MFF, NDUFS3, NDUFV2) and 36 RBPs were included.

3. Discussion

The prevalence of Thyroid Eye Disease (TED) among Graves' disease patients is approximately 25%, making it the most prevalent orbital condition. Existing treatments such as medications, surgery, and radiation have not consistently yielded satisfactory outcomes and can potentially result in complications and adverse effects. While numerous studies have explored the pathogenesis of TED and potential treatment targets, the precise mechanisms remain elusive. Contrasting gene expression profiles between TED and other orbital inflammatory conditions have shown weaker associations with inflammatory markers in TED, indicating a relatively subdued inflammatory response in TED[19], This implies that targeting metabolic or alternative pathways could be more effective in treating Thyroid Eye Disease (TED). Dysregulation in mitochondrial energy metabolism pathways is linked to oxidative stress. Inflammation is implicated in a range of degenerative conditions, acute illnesses, and the aging process. Mitochondria are key players in amplifying inflammatory signals, leading to heightened mitochondrial oxidative stress and perpetuating a detrimental cycle of inflammation[43]. There is currently no definitive target or treatment mechanism in the study of mitochondrial energy metabolism-related genes in TED, which holds potential and guiding significance for the diagnosis and treatment of TED.

Following GO and KEGG analyses, 26 differentially expressed mitochondria-related genes (MEMRDEGs) were identified as enriched in various biological processes, cellular components, molecular functions, and pathways in Thyroid Eye Disease (TED). The primary aim of this study is to evaluate the diagnostic potential of these 26 MEMRDEGs in TED. While the specific regulation status (upregulation or downregulation) of these genes in TED was not detailed, their enrichment in the disease context suggests diagnostic relevance. Our analysis focuses on the collective impact of these genes in TED pathogenesis. Gene Set Enrichment Analysis (GSEA) revealed significant enrichment of gene expression levels related to TED in dataset GSE58331, particularly in pathways involving fatty acids, proinflammatory and profibrotic mediators, and ADORA2B-mediated anti-inflammatory cytokine production. Mitochondrial energy metabolism pathways encompass glycolysis, the tricarboxylic acid cycle, oxidative phosphorylation, ketone body metabolism, lipid metabolism, and gluconeogenesis. Oxidative phosphorylation, crucial for ATP synthesis, stands out as a central pathway in energy metabolism. The activation of the ADORA2B receptor triggers the cyclic adenosine monophosphate (cAMP) pathway and associated downstream pathways, linked to immunosuppression and potential applications in asthma and tumor immunotherapy[44]. GSEA indicated significant enrichment related to oxidative phosphorylation and ATP energy production within mitochondria. By categorizing high and low-risk groups using Logistic risk scoring, GSEA demonstrated significant enrichment of TED genes across pathways such as Signaling By Interleukins, IL18 Signaling Pathway, and other relevant biological functions. In the context of intracytoplasmic mitochondrial DNA and ROS, the activation of inflammatory vesicles can lead to IL-1β and IL-18 secretion upon release from dysfunctional mitochondria[45]. We suggest that the mitochondrial energy metabolism pathway in TED involves inflammasome signaling pathways. These findings suggested that MEMRGs were involved in the immune aspects of TED and were associated with the pathogenesis of TED. Therefore, we further compared the immune infiltration abundance in TED grouped by high and low risk score.

Using the ssGSEA algorithm to assess immune cell infiltration abundance in Thyroid Eye Disease (TED) high and low-risk groups, significant findings (p < 0.05) were observed for Activated CD8 T cells, Immature B cells, and T follicular helper cells. Key Genes displayed strong positive correlations with immune cells in both Low Risk and High Risk TED groups. Notably, UQCRFS1 exhibited the strongest positive correlation with T follicular helper cells in the Low Risk group (r = 0.529, p < 0.05), while MDH1 showed the strongest positive correlation with T follicular helper cells in the High Risk group (r = 0.733, p < 0.05). Previous immunohistochemical studies have linked cytokines such as interferon-γ (IFN-γ), tumor necrosis factor-α (TNF-α), and interleukin-1α (IL-1α) with T-cell infiltration in TED connective tissues. During the Active phase of TED, CD4+ and CD45RO+ T cells are abundant in early orbital infiltration, secreting cytokines and chemokines to enhance the immune response. Follicular helper T cells (Tfh) play a crucial role as regulators of germinal center (GC) B cells and in sustaining long-term protective humoral responses[46,47]. Tfh cells aided in immunoglobulin class switch reorganization and supported germinal center responses, thereby promoting immunoglobulin affinity maturation and the generation of humoral immune memory[48]. While their primary function was to promote B-cell responses, Tfh cells also exhibited phenotypic and functional diversity depending on the immunologic and spatial environment they encounter[49]. UQCRFS1 was highly expressed in gastric and breast cancers. Knocking out UQCRFS1 cells can reduce cell proliferation, induce cell cycle arrest in G1 phase, increase the proportion of apoptosis, ROS production and DNA damage gene expression, and inhibit the AKT/mTOR pathway[50]. Dysregulation of Tfhs and follicular regulation T cells had been reported to be positively associated with the pathogenic process of autoimmune diseases[51].

The LASSO regression model identified 10 key genes in Thyroid Eye Disease (TED): IDH1, IMMT, MDH1, NDUFS3, NDUFV2, GPI, UQCRFS1, GCK, ALDH1A1, and MFF. ROC curves revealed that IDH1, IMMT, MDH1, NDUFS3, NDUFV2, ALDH1A1, and MFF expression levels are highly accurate in distinguishing between the TED and Control groups (AUC > 0.9). NDUFS3, highlighted by Friends analysis, was deemed crucial in TED's biological processes (cut-off value = 0.60). In invasive ductal breast carcinoma samples, the mitochondrial complex I subunit NDUFS3 exhibited a distinct signature as a marker for hypoxia/necrosis. NDUFS3 expression levels were significantly and positively associated with nuclear grading, enabling improved diagnosis and survival rate predictions for breast cancer patients[52]. Furthermore, Complex I (CI) is the largest enzyme in the mitochondrial respiratory chain, and its defects were a major cause of mitochondrial disease. CI disassembly processes were identified as occurring in a hierarchical and modular manner when the amount of NDUFS3 gradually decreases[53]. Inner membrane mitochondrial protein (IMMT) is crucial for preserving mitochondrial structural integrity, regulating mitochondrial dynamics, and contributing to various cellular functions. It plays a role in cell cycle advancement and mitochondrial antioxidant protection. Elevated IMMT levels are linked to lung adenocarcinoma and independently predict a poor prognosis in patients with this condition. IMMT has the potential to serve as a novel prognostic marker and therapeutic target in lung adenocarcinoma[54]. IMMT played a crucial role in shaping the inner mitochondrial membrane, regulating metabolism and innate immunity. Low expression of IMMT in tumors was associated with mitochondrial inhibition and activation of angiogenesis, indicating a poor prognosis for patients with kidney renal clear cell carcinoma (KIRC) and being linked to the progression of KIRC[55]. NDUFS3 and IMMT exhibited high expression levels in Thyroid Eye Disease (TED) datasets and were linked to mitochondrial energy metabolism pathways, potentially influencing TED's immune responses. These genes are considered promising research directions and potential therapeutic targets for TED. However, our study has limitations. While we identified associations between 26 differentially expressed genes and TED, lack of biological validation hinders a comprehensive understanding of their roles in the disease's pathophysiology. Resource constraints prevent further wet lab experiments, and the absence of clinical relevance studies limits analysis with clinical data. Future research aims to address these gaps through experimental validation. Despite combining data from GSE105149 and GSE58331 to enhance statistical power, the relatively small sample size may impact result reliability. Relying on publicly available data introduces potential bias, necessitating future studies with independent datasets and experimental confirmation to bolster our findings.

4. Conclusions

The study examined differential expression of mitochondria-related genes (MEMRGs) in Thyroid Eye Disease (TED) and Control groups, identifying 10 Key Genes for a TED diagnostic model. Immunoinfiltration analysis demonstrated a notable positive correlation between these key genes and immune cells in TED, offering valuable insights for TED diagnosis and prognosis. Nevertheless, additional validation is necessary to uncover precise pathogenic mechanisms and target effects, necessitating thorough investigation.