GEO (https://www.ncbi.nlm.nih.gov/geo/) is a high-throughput gene expression database that contains microarray, second-generation sequencing, and other high throughput sequencing data. GSE58331, which was established on the platform of GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) was obtained by searching for thyroid associated ophthalmopathy in GEO. The raw data and annotation files of GSE58331 were acquired from the GEO database. Of the 175 samples in this dataset, 27 and 22 samples of orbital tissue were from TAO patients and healthy people, respectively. These samples were used as the train group to identify DEGs and differential expressed FRGs. Seven TAO and 8 healthy samples of lacrimal gland tissue were selected as the test group to plot ROC curves. We collected the gender information of the samples used in our work and presented in the Supplementary Table S1. The flow chart of the whole research was shown in Figure 1.

ID conversion of raw data was accomplished using Perl based on the annotation files. The Limma package was used to screen out DEGs, setting the selection criterion at |log₂Fold Change|(log₂FC)>0.25 and adjusted p-value<0.05. FerrDb (http://www.zhounan.org/ferrdb/current/) is the first manually compiled database which is dedicated to ferroptosis regulators and ferroptosis-disease associations. A set of ferroptosis drivers, suppressors, and marker genes were downloaded from FerrDb. To identify differentially expressed FRGs, we intersected genes in the obtained sets with DEGs. A volcano plot and heatmap were employed to visualize the DEGs and differentially expressed FRGs, respectively.

The ROC curve was plotted in the R studio using the R package “qROC” to evaluate the diagnosis value of the screened FRGs in TAO. The area under curve (AUC) was calculated to quantitatively assess the power of prediction accuracy. An AUC value more than 0.7 indicated high diagnostic value for the disease.

We performed immune cells infiltration analysis using the CIBERSORT algorithm. This algorithm was developed by Newman et al. and is used to estimate the abundance of 22 types of immune cells in gene expression profile data. The source code of CIBERSORT and LM22. txt which includes signature genes and their expression profile in 22 types of immune cell were downloaded from CIBERSORTx (https://cibersortx.stanford.edu/). We calculated the p-value for each sample to estimate their heterogeneity of immune infiltration by using the permutation test inside the CIBERSORTx. One sample with the p-value greater than 0.05 indicates that there is no statistical difference in immune cell infiltration and would be excluded in the follow-up analysis. Spearman correlation analysis was applied to analyze the correlation between differentially expressed FRGs and immune cells in TAO samples with the filter of Spearman’s R > 0.4 or < -0.04 and p-values <0.05. The analysis was also performed using the R packages (“corplot,” “vioplot,” “ggplot2,” and “glment”).

The human genome annotation file GRCh38 was obtained from Genecode (Frankish et al., 2021) and used to distinguish lncRNAs from mRNAs in the gene matrix. The Limma R package was used to identify differential expressed lncRNAs in TAO groups. The Pearson correlation analysis was utilized to investigate the co-expression of differential expressed lncRNAs and FRGs in the gene matrix. The absolute value of Person’s R > 0.4 and p-values <0.05 were considered to indicate a significant co-expression relationship. Differential expressed lncRNAs which showed a significant co-expression relationship with differential expressed FRGs were regarded as ferroptosis-related lncRNAs (FR-lncRNAs).

The DGIdb (Drug-Gene Interaction database) was utilized to predict drugs targeted by the FRGs identified in our study. The DGIdb is a drug-gene interaction database that provides information on the association of genes with their known or potential drugs and has been updated to version 4.0 (Freshour et al., 2021). More detailed information about these drugs were obtained from the DrugBank database (Wishart et al., 2018).

TRRUST was used to predict the transcription factors that regulate related FRGs. TRRUST (https://www.grnpedia.org/trrust/) is an artificially annotated database of transcriptional regulatory networks for human and mouse which not only contains the target genes corresponding to transcription factors, but also the regulatory relationships between transcription factors.

Orbital tissues were obtained from TAO patients who underwent orbital decompression surgery (n = 3) and healthy individuals who received oculoplastic surgery (n = 3) in the Shanghai Changzhen Hospital. This program has been reviewed and approved by Committee on Ethics of Biomedicine, Second Military Medical University. Each sample was collected with acquiring informed consent and strictly subjected to 1964 Helsinki declaration. Primary OFs were isolated from orbital tissues as mentioned in the previous study (Diao et al., 2020). Briefly, the orbital tissue obtained during surgery was separated into small pieces and neatly placed in a Petri dish containing small amount of DMEM (10% FBS, L-glutamine, 110 mg/L sodium pyruvate, 1% penicillin/streptomycin). Long spindle-shaped fibroblasts could be observed migrating out of the tissue mass after approximately 7 days. Passage and phenotypic validation were performed when the primary fibroblasts grew to confluence. The human OFs at 3-8 passage was selected for the experiment. Total RNA of OFs was extracted using the RNA-Quick Purification Kit (RN001, Yishan, Shanghai, China) and reverse transcribed to cDNA by HiScript^® Q RT SuperMix for qPCR (R122-01, Vazyme, Nanjing, Jiangsu, China). The primer sequence of six FRGs, two lncRNAs and GAPDH used as the housekeeping gene are presented in the Supplementary Table S2. After amplification of cDNA by ChamQ Universal SYBR qPCR Master Mix (Q711-02, Vazyme, Nanjing, Jiangsu, China), the fluorescence signal during the exponential amplification phase was collected and analyzed. All results were analyzed as three replicates from three independent experiments.

A total of 162 DEGs (32 upregulated and 130 downregulated) were identified based on the criteria of log₂|FC|>0.25 and adjust p-value<0.05 ( Figure 2A and Supplementary Table S3). Six differentially expressed ferroptosis driver genes including CYBB, CTSB, SLC38A1, TLR4, PEX3, and ABCC1 were identified through the intersection of DEGs and ferroptosis-related genes (Table 1). The heatmap showing their expression profile in each sample and the data variation between the two groups is shown in Figure 2B. Except for SLC38A1, the other five genes were downregulated in the orbit tissue of TAO patients.

The ROC curve of six FRGs based on GSE58331 is presented in Figure 2C.

To further determinate the diagnose power of these genes in disease, we plotted the ROC curve in both 49 orbit tissue samples (27 TAO and 22 healthy) and another 15 lacrimal gland samples (7 TAO and 8 healthy) from GSE58331. The AUC value of differentially expressed FRGs was calculated to perform quantitative analysis. The AUC values for six genes in the orbital tissue group were all higher than 80 (Figure 2C). In lacrimal gland tissues groups, the AUC value of CTSB (AUC = 0.875; CI: 0.625–1.000), SLC38A1 (AUC = 0.893; CI:0.679–1.000), and PEX2 (AUC = 1.000; CI:1.000–1.000) was greater than 80, which indicated high clinical diagnostic power (Figure 3). Meanwhile, we plotted the ROC of TSHR and IGF1R which is the existing TAO molecular markers in both orbit tissue and lacrimal gland based on GSE58331 and compared their AUC value with FRGs in this study. The AUC value of TSHR and IGF1R were not greater than six FRGs in tissues at transcriptional levels (Supplementary Figure S1).

The violin plot (Figure 4A) of immune cell infiltration indicated high infiltration level of monocytes (p < 0.001), macrophages M0(p = 0.039), mast cells activated (p = 0.008), and neutrophils (p = 0.045) in TAO samples. Meanwhile, mast cells resting (p = 0.043) and macrophages M2 (p = 0.02) which can suppress T cell-mediated immune responses exhibited low infiltration level in TAO samples (Table 2). No gender differences regarding the immune microenvironment were observed in the TAO samples (Figure 4B). The bar plot (Figure 5A) revealed the relative percentage of 22 types of immune cells in each sample. The correlations among 19 types of immune cells (T cells CD4 memory activated, T cells gamma delta, and dendritic cells activated were not included in the analysis because they did not show infiltration in all samples) are presented in the correlation heatmap (Figure 5B). Early inflammatory response in TAO is characterized by monocytes infiltration. In our results, monocytes showed a significant negative correlation with macrophages M0 and mast cells activated but showed a positive correlation with macrophages M0, mast cells activated, neutrophils, and T cells CD4 naïve. Besides, macrophages M2 and monocytes (R = −0.74), mast cells resting and plasma cell (R = −0.7), mast cells resting, and mast cells activated (R = −0.66), T cells follicular helper and mast cells resting (R = −0.58), natural killer (NK) cells activated and NK cells resting (R = −0.55) showed negative relationships. There was a positive relationship between mast cells resting and macrophages M2 (R = 0.5), mast cells activated and plasma cell (R = 0.65), NK cells resting and neutrophils (R = 0.47).

To explore the potential function of differential expressed FRGs (CYBB, CTSB, SLC38A1, TLR4, PEX2 and ABCC1) in the pathogenesis of TAO, we further analyzed the relationship between these genes and 18 immune cells (T cells CD4 memory activated, T cells gamma delta, eosinophils, and dendritic cells activated were not included in the analysis because they did not show infiltration in TAO samples) in TAO groups. The results were presented in a correlation heatmap and were selected based on Spearman’s R > 0.4 and p-value<0.05 (Figure 6). Noticeably, ABCC1 exhibited a significant strong positive correlation with mast cells resting (R = 0.71) and neutrophils (R = 0.73), but a negative correlation with plasma cells (R = −0.87). CYBB showed a significant strong positive correlation with macrophages M2 (R = 0.75) and mast cells resting (R = 0.83), but a negative correlation with monocytes (R = −0.7) and macrophages M0 (R = 0.77). A significant strong positive correlation was observed between PEX2 and mast cells resting (R = 0.72).

By applying the criteria of Pearson’s R > 0.4, two lncRNAs (LINC01140 and ZFHX4-AS1) were identified from the DEGs. LINC01140 exhibited a significant positive correlation with CYBB (R = 0.57, p < 0.01) and TLR4 (R = 0.46, p < 0.05), but a negative correlation with SLC38A1 (R = −0.64, p < 0.001). ZFHX4-AS1 showed a significant positive correlation with CTSB (R = 0.64, p < 0.001) and CYBB (R = 0.53, p < 0.01) (Figure 7; Table 3). CYBB-LINC01140-TLR4, CYBB- LINC01140- SLC38A1, TLR4- LINC01140- SLC38A1, and CTSB- ZFHX4-AS1- CYBB were considered to be potential RNA regulatory pathways involved in the pathogenesis of TAO.

The DGIdb database was used to predict drugs that interact with the six ferroptosis related genes. Fifty targeted drugs for ABCC1, 11 targeted drugs for TLR4, 3 targeted drugs for CTSB and 3 targeted drugs for CYBB were identified in the database. However, drugs that were targeted by SLC38A1 and PEX2 were not obtained. These results were visualized by the Cytoscape software as a gene-targeted drugs network shown in Figure 8A. The detailed information, including types, chemical formulas, and structure diagrams of these drugs is presented in the (Supplementary Table S4). Notably, 21 transcription factors for four FRGs (CYBB, CTSB, TLR4, ABCC1) were screened out in the TRRUST database. Nine transcription factors for CYBB, 8 transcription factors for CTSB, 5 transcription factors for TLR4 and 2 transcription factors for ABCC1 were also identified (Figure 8B). Among these, CYBB and TLR4 share two common transcription factors (SPI1 and IRF8) whereas CTSB and ABCC1 share one common transcription factor (SP1). Transcription factors that regulate SLC38A1 and PEX2 were not obtained.

We further validated the expression of the identified differentially expressed FRGs (CYBB, CTSB, SLC38A1, TLR4, PEX2 and ABCC1) and two FR-lncRNAs (LINC01140 and ZFHX4-AS1) in cultured OFs from TAO patients and healthy controls (Figure 9). The results indicated that the expression of CTSB, PEX2, ABCC1 and ZFHX4-AS1were significantly down-regulated in TAO groups which is consistent with our analysis results. Though not significant, the expression of CYBB and TLR4 tended to be decreased in OFs from TAO patients.