A single cell dataset of GSE134420 was obtained from the Gene Expression Omnibus (GEO) public database for bioinformatic analysis., A total of 23 transcriptomic data sets in GSE77298, including 7 control samples and 16 RA samples, were downloaded using GPL570 platform. Differential genes (DEGs) between controls and RA were identified using R (version 3.5.1), and DEGs that were P < 0.05 and |logFC| >1 was selected for subsequent analysis. A dataset of 619 ferroptosis-related genes was retrieved from the GeneCards database (https://www.genecards.org/).

The use of the Metascape database (www.metascape.org ) for annotation and visualization, as well as the GO analysis and KEGG pathway analysis for differentially expressed genes, is a common approach to understand the biological functions associated with disease progression. Setting a minimum overlap of 3 and a p-value threshold of 0.01 is a reasonable criterion for determining statistical significance.

The PPI network of mapped genes was constructed by using the Retrieval of Interacting Genes (STRING) database (http://string-db.org/). Setting a threshold interaction score of 0.4 is a reasonable criterion for determining the strength of the interactions between genes. with a threshold interaction score of 0.4, which was visualized by the Cytoscape software. To understand the complex interactions between genes and identify key nodes in the network, Cytoscape software (version 3.7.2) was used to visualize the PPI network.

CIBERSORT algorithm is widely utilized to infer the relative degree of the 22 immune infiltrating cells in the microenvironment. The method is based on the principle of support vector regression (SVM), and the expression matrix of immune cell subtypes is analyzed by deconvolution. In this study, patient data was analyzed using the CIBERSORT algorithm, and the result was then used to perform Spearman correlation analysis of gene expression as well as immune cell content. By correlating gene expression with immune cell content, this analysis can provide insights into the role of different immune cell types in disease development and progression.

The Gene Atlas database(http://geneatlas.roslin.ed.ac.uk/) is a comprehensive resource that provides annotations against genes and proteins using the UK Biobank cohort. The database is a valuable tool for understanding the genetic basis of complex diseases and traits. The associations recorded in the Gene Atlas database have been calculated using data from 452,264 individuals from the UK Biobank. This large dataset allows for the identification of genetic variants that are associated with specific traits or diseases, providing insights into the underlying biological mechanisms.

GSEA analysis provides insight into the biological processes that are dysregulated in a particular phenotype or disease state by ranking genes based on their differential expression using a set of pre-defined genes, such as a pathway or gene ontology term. The method then tests whether the predefined set of genes is enriched at the bottom or top of this ranking table, indicating whether the pathway or gene ontology term is associated with the phenotype being studied. Here, we use GSEA to compare the differential signaling pathways in both the high- and low-expression group in order to investigate the underlying mechanism of key genes in both groups. Number of substitutions were 1000 and the substitution was phenotype, which means that phenotype labels were randomly permuted during the analysis.

Here, we use the RcisTarget in R to predict TFs. RcisTarget calculates the AUC for each pair of motif-motif. This is the first stage to estimate the over-expression of each motif. This is done based on the calculation of recovery curve according to the order of the motifs. The normalized enrichment score (NES) is computed from AUC distributions The gene-motif ranking database used by RcisTarget is rcistarget. hg19. motifdb. cisbpont.500bp.

Data was primarily processed by the Seurat package, and position relationships among clusters were determined using the tSNE algorithm. The clusters were annotated by Celldex package, and some cells with relationship to RA progression were separately annotated. Finally, marker genes were extracted from single cell expression profiles with the Find All Markers logfc. threshold parameter set at 1. For each cell subtype, gene was selected to be a specific marker with |avg log₂FC| > 1 and p val adj< 0.05.

The CMap is an interventional gene expression profiling database s that will be primarily used for the discovery of functional associations between small molecules, gene and disease condition. This study predicts disease-targeted therapeutics through differentially expressed genes in disease.

Primary RA synovial fibroblasts were isolated from the synovium of male Sprague-Dawley rats weighing 160-180 g purchased from the Beijing Vital River Laboratory Animal Technology Co., Ltd. This study was approved by the Animal Ethics Committee of Peking University First Hospital. The rat model of adjuvant-induced arthritis (AIA) was established according to a previously published protocol (15). The knees of rat models were embedded in paraffin, sectioned at 5-μm thickness sections, and then stained with hematoxylin and eosin.

The methods of primary synovial fibroblast extraction and culture were based on a previous publication (16). In order to obtain pure synovial fibroblasts, cells used herein were at passage 3-5. Synovial fibroblasts were cultured in DMEM high-glucose medium (Gibico, America) containing 10% fetal bovine serum (Pan, Germany) and 1% penicillin-streptomycin (Beyotime, China) at 37 °CC in a 5% CO2 incubator.

Synovial specimens were isolated from 3 RA patients who satisfied the diagnostic criteria of 2010 ACR and underwent total knee replacement surgery, and3 OA (osteoarthritis) patients who fulfilled the diagnostic criteria of 2010 ACR/EULAR and were undergoing surgery at the Peking University First Hospital. This study was approved by the Medical Ethics Committee of Peking University First Hospital. Detailed information on the clinical details is given in Table S1 .

Paraffin-embedded samples were sliced into 8 μm sections. Antigen retrieval was performed using the target retrieval solution Tris-EDTA buffer with a microwave oven for 10-15 min. Tissues were immersed in 3% H₂O₂ for 30 min. 3% BSA was added to evenly cover the tissues for 30 min. Sections were immersed with antibodies against PTGS2(1:200, K001561P, Solarbio), ENO1(1:200, ab227978, Abcam), GRN (1:200, DF7997, Affinity Biosciences) overnight at 4°C.

Human primary FLSs were isolated from patients’ synovial tissues, minced and digested in type I collagenase (2 mg/ml, Solarbio, China) and 0.25%Trypsin-EDTA Solution (Solarbio, China) for 4 h. The resulting cell suspension was centrifuged and the cells were resuspended in medium. After 24 h of incubation, PBS washes removed nonadherent cells.

According to the GenBank database, the siRNA sequence( Table S2) of rat and human ENO1 gene designed and obtained from JTSBIO Co,Ltd (Wuhan, China). ENO1-siRNA and negative control (NC) vectors were transfected into primary cultured FLS cells of rat and human.

The knees of rat models were embedded in paraffin, sectioned at a thickness of 5-μm, and then stained with hematoxylin and eosin. Cells were lysed using RIPA buffer (Beyotime, China). Protein was electrophoresed with SDS-PAGE gel and transferred to PVDF membrane, blocked using QuickBlock™ Primary Antibody Dilution Buffer (Beyotime, China), and then incubated with the ENO1 antibodies (ab227978, Abcam), GRN antibodies (DF7997, Affinity Biosciences), PTGS2 antibodies (K001561P, Solarbio) and ACSL4 (ab155282, Abcam) at 4°C overnight. Membrane was subsequently incubated with HRP-conjugated secondary antibodies. BeyoECL Moon (Beyotime, China) was used to visualize the protein bands. Detection was performed using an automated chemiluminescence image analysis system (SyngeneG Box chemi xx6, British) and the quantification of the bands was performed using Image J software.

Cells were cultured in 96-well plates in 5% CO2 incubator. The medium containing the drug was replaced and add 10 microliters of CCK-8 solution (C0038, Beyotime) to each well. A blank control can be used with the corresponding amount of medium and CCK-8 added but no cells. The cells were then incubated in the cell incubator for 1.5 hours. Absorbance was measured at 450nm through Multifunctional Enzyme Labeler (SpectraMax i3).

Cells were inoculated in fluorescent culture dishes and incubated overnight in 5% CO2 incubator. The medium containing the drug was replaced, and FerroOrange working liquid (F374- DOJINDO) with a concentration of 1 μmol/l was added and incubated at 37°C in 5% CO2 incubator for 30 minutes. Fluorescent was observed under a fluorescence microscope.

To detect the level of reactive oxygen species (ROS), the medium was removed and DCFH-DA (10μmol/l, S0033S, Beyotime) was added and incubated in the incubator at 37°C for 20 minutes. The cells were washed to fully remove the excess DCFH-DA. ROS level was detected through flow cytometry (BECKMAN COULTER- CytoFLEX SRT).

A total of 23 patients, including 7 in the control group and 16 in the disease group, were downloaded from the NCBI GEO public database for the GSE77298 dataset. The limma package was used to count the differential genes between both sample groups with the differential gene screening conditions of P.value<0.05 and |logFC|>1. In total, 1886 differential genes (990 up- and 896 down-regulated genes) were finally screened ( Figure 2 ). We mapped the differential genes to the ferroptosis gene set and obtained 49 overlapping genes, including 28 up- and 21 down-regulated genes.

We then conducted pathway analysis of the 49 intersecting genes, and found that these genes were primally involved in the pathways of small molecule biosynthetic process, metal ion response, and oxidoreductase activity ( Figure 3 ). Protein interaction relationship pairs related to the 49 intersecting genes were retrieved from STRING database and visualized using Cytoscape, and the top 3 genes of degree were selected as key genes, which were PTGS2, ENO1, GRN ( Figure 4 ).

The immune microenvironment plays a crucial role in the pathogenesis of RA and influences its diagnosis, prognosis, and clinical response. By analyzing the relationship between PTGS2, ENO1, GRN and immune infiltration in the disease dataset, the mechanism by which these hub genes influence the progression of RA were explored. The level of immune cell in each sample was shown in Figure 5A . There were several pairs of significant correlation between the levels of immune infiltration ( Figure 5B ). And the levels of naive B cells, resting CD4 memory T cells, and activated NK cells were significantly lower in the disease group samples compared to normal patients ( Figure 5C ). We further explored the association between key genes and immune cells and found that several key genes had strong correlation with immune cells, and ENO1 had strong correlation with M0 macrophages, T follicular helper cells, and activated mast cells, and with resting mast cells, T regulatory cells (Tregs), resting dendritic cells, etc. ( Figure 5D ); GRN had a positive correlation with T gamma delta cells, plasma cells, activated CD4 memory T cells, etc., and negatively correlated with resting dendritic cells, resting CD4 memory T cells, and resting mast cells ( Figure 5E ); PTGS2 had a positive correlation with M0 macrophages and activated mast cells, and showed a negative correlation with naive B cells, resting mast cells, and Tregs ( Figure 5F ).

We retrieved 4881 RA-related pathogenic genes from the GeneCards database (https://www.genecards.org/ ), and analyzed the expression levels of the top 20 genes and found that the expression of HLA-B, MIF, PSTPIP, TLR1 differed in the two groups ( Figure 6A ). We then carried out a correlation analysis of key genes and genes regulating RA. and found that the expression levels of hub genes were significantly correlated with several RA-related genes. Specifically, ENO1correlated significantly and positively with MIF (Pearson r=0.82) and GRN correlated significantly and positively with HLA-B (Pearson r=0.8) ( Figure 6B ).

Next, to confirm the pathogenic regions of three hub genes in RA, we analyzed the RA GWAS data. We also showed the pathogenic regions of the SNPs corresponding to PTGS2, ENO1, and GRN, with PTGS2 and ENO1 located in chromosome 1 pathogenic region, and GRN in chromosome 17 pathogenic region ( Figure 7 ). Table S3 (GWAS data.xls) shows the significant SNP loci corresponding to the 3 genes.

Next, to explore the underlying mechanisms of how these key genes affect disease development, we investigated the specific pathways associated with the three key genes. GSEA results revealed that the pathways enriched by the highly expressed PTGS2 gene are Nod-like receptor signaling pathway, FcγR mediated phagocytosis, Chemokine signaling pathway and other pathways ( Figure 8A ); high expression of GRN gene was enriched in Lysosome, Leishmania infection, B cell receptor signaling pathway and other pathways ( Figure 8B ); high expression of ENO1 gene was enriched in Lysosome, Pathogenic Escherichia coli infection, Leishmania infection and other pathways ( Figure 8C ).

The 3 hub genes used for the analysis were found to be under the control of multiple TFs, as well as other common mechanisms. Therefore, an enrichment analysis was conducted for these TFs using cumulative recovery curves ( Figure 9 ), motif TF annotation and key gene screening. The result revealed that the motif cisbp_M5493 had the highest normalized enrichment score (NES)of 6.17. Genes of ENO1 and PTGS2 were enriched in this motif. We showed the enriched motifs and the corresponding TFs ( Supplement Figure 1 ). We also reverse predicted the 3 key genes using the miRcode database and obtained 56 miRNAs, a total of 80 mRNA-miRNA relationship pairs, and visualized them using Cytoscape ( Figure 10 ). In addition, we utilized the miRcode database to reverse-predicet the three key genes and identify 56 miRNAs. We then visualized the resulting 80 mRNA-miRNA relationship pairs using Cytoscape, as shown in Figure 10 .

Cells were clustered using the tSNE algorithm, and the expression of three key genes, prostaglandin-endoperoxide synthase 2 (PTGS2), enolase 1(ENO1), and granulin (GRN), in the nine cell clusters is shown in Figure 11 .

We utilized the Connectivity Map database to predict potential drugs for the disease. Our analysis revealed that the expression profiles of CID_11676786 (AZ628), mycophenolic acid, sulmazole, and indatraline drug perturbations had the most significant negative correlations with disease perturbation expression profiles. These findings suggest that these drugs could alleviate or even reverse the disease state ( Figure 12 ).

The adjuvant-induced arthritis (AIA) rat model is a well-established and widely used method for establishing preclinical model of RA. After 25 days of exposure to Complete Freund’s adjuvant (CFA), the establishment of the model was evaluated by calculating the polyarthritis index (PI), joint swelling thickness and joint histopathology. The results demonstrated that model rats exhibited a significant increase in PI and joint swelling thickness, compared to control rats ( Figure 13A ). Furthermore, joint histopathological examination revealed that the synovial tissue in the model group was markedly hyperplastic, with thickened cells in the synovial lining layer, obvious inflammatory cell infiltration, and pannus formation ( Figure 13B ). We isolated and cultured primary FLS from rats and detected their protein expression levels using WB. As shown in Figure 13C , the expression of GRN, ENO1 and PTGS2 were significantly upregulated in the model group, which was consistent with the results in the RA patient datasets. To validate our findings in human samples, we further investigated the expression of GRN, ENO1 and PTGS2 through IHC in patient samples and found that these genes were dramatically elevated in RA compared to OA as a control ( Figure 13D ).

Since PTGS2 has been reported as a specific biomarker for ferroptosis (10), we focus on ENO1 to validate its potential function. ENO1, as an RNA-binding protein, can accelerate the messenger RNA decay of ACO1 in cancer cells, leading to repression of ferroptosis (15). ACO1 plays a critical role in the regulation of ferroptosis, which is involved in the activation of ferroptosis by ROS and the upregulation of genes that promote iron uptake and storage (16). We examined the expression of ENO1 and ACO1 in the synovium of RA patients, finding that ENO1 was upregulated while ACO1 was downregulated ( Figure 14B ). To further explore this pathway, we knocked down ENO1 in primary FLS of a rat model ( Figure 14A ), and found that the expression of ACO1 and the signature protein for ferroptosis, ACSL4, were increased. This effect was more pronounced in the presence of iron, and was reversed by treatment with ferrostatin-1(Fer-1) and Liproxstatin-1(Lip-1), both well-known ferroptosis inhibitors ( Figure 14C ). We also knocked down ENO1 in primary FLS of RA patients, leading to the upregulation of ACO1 and ACSL4. However, this effect was not observed to be more pronounced in the presence of iron. Treatment with Fer-1 and Lip-1 was still found to eliminate the increased ACSL4( Figure 14D ). Moreover, consistent with previous reports that increased ACO1 led to excess accumulation of ROS and cell death in cancer cells (16, 17), knockdown of ENO1 led to increased cell death ( FigureS 14E, F ), enhanced lipid ROS generation ( Figures 15A, B ), and increased accumulation of intracellular ferrous ion ( Figure 15C ), particularly in the presence of iron. Treatment with Fer-1 and Lip-1 abolished ENO1-ACO1-mediated intracellular ferrous ion, ROS accumulation, and cell death ( Figures 15A–C ).