Introduction

Rheumatoid arthritis is a chronic inflammatory autoimmune disease characterized by persistent synovitis and can destroy bone and cartilage tissue, leading to joint damage, chronic disability, and a shortened life span.1–3 The global prevalence of RA is approximately 1%.4 The underlying molecular mechanism of its pathogenesis is still unknown, and the onset and progression of the disease are closely related to a variety of factors (genetic, environmental, and immunologic), among which an abnormal infiltration of immune cells has a significant impact on RA.5,6 Although significant progress has been made with several targeted immunotherapies,3 the limited therapeutic efficacy prompts us to investigate alternative pathogenic mechanisms and pathways, among which programmed cell death (PCD) is an innovative and major breakthrough. The imbalance of cell death in FLS and immune cell populations and the release of inflammatory cytokines ultimately promote synovitis, contributing to the RA’s deterioration.7 Apoptosis,8 autophagy,9 NETosis,10 necroptosis,11 ferroptosis,12 and pyroptosis13 are PCDs affecting the initiation and development of RA. RA may be an outcome of the interaction of multiple PCDs, and few studies have examined how ferroptosis and pyroptosis interact to exacerbate the course of RA. Thus, our study is looking for common disease targets.

Ferroptosis is an iron-dependent PCD14 characterized by the interaction between iron and lipid peroxides to produce reactive oxygen species (ROS).15 Studies have shown that RA disrupts iron metabolism and has iron deposits in the synovial tissue absent in healthy synovium.16 Oxidative stress and lipid peroxidation reduce Fe3+ to Fe2+ and the production of free radicals (-OH, RO-), leading to cellular volume expansion and plasma membrane rupture.17 It has been demonstrated that there is a significant positive correlation between the ROS level and the severity of RA.18 Additionally, immune cells that undergo ferroptosis may release proinflammatory cytokines that eventually trigger and worsen inflammatory disorders.19 Therefore, ferroptosis may be a key driver of autoimmune and inflammatory diseases.20 Pyroptosis is also a newly discovered type of PCD,21 characterized by pore formation, cell enlargement, and plasma membrane rupture.22 Pyroptosis exacerbates the disorders by inducing plasma membrane rupture and release of inflammatory cytokines in cells with RA via the inflammasome signaling pathway,23 in which caspase proteins play a crucial role.24

Quercetin is a natural flavonoid with potent antioxidant and anti-inflammatory activity.25 In both human and animal models of rheumatoid arthritis, studies have demonstrated that quercetin reduces symptoms and delays disease progression.26 Additionally, quercetin can alleviate seizures,27 acute kidney injury,28 liver injury,29 and intestinal injury30 by inhibiting ferroptosis. Quercetin also attenuates atherosclerosis,31 cardiomyopathy,32 lung injury,33 and kidney injury34 by inhibiting pyroptosis. Therefore, quercetin might be an ideal drug for RA by targeting ferroptosis and pyroptosis. Our study identified a common biomarker of ferroptosis and pyroptosis in RA through bioinformatics, in vitro experimental validation, and a potential drug targeting ferroptosis and pyroptosis by network pharmacology and computer-aided drug design (CADD), which will provide new insights for experimental study and clinical diagnosis and treatment.

Methodology Data Collection and Processing

Two datasets for RA were included, GSE5523535 and GSE77298,36 which were screened using the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) as shown in Table 1. The GSE55235 dataset contains synovial tissue samples from 10 patients with rheumatoid arthritis and 10 normal subjects for deg identification. The GSE77298 dataset contains synovial tissue samples from 16 rheumatoid arthritis patients and 7 normal subjects to validate the hub genes. From the FerrDb database (http://www.zhounan.org/ferrdb/current/),37 382 ferroptosis-related genes (FRGs) were obtained, including 150 driver genes, 109 suppressor genes, and 123 marker genes. In addition (Table S1), 442 genes associated with ferroptosis were searched in the genecards database (https://www.genecards.org/) (Table S2).38 The integration of these two gene banks screened 592 FRGs. The GeneCards database was applied to acquire the 254 pyroptosis-related genes (PRGs) (Table S3).

Table 1 Basic Information of Selected Datasets

Identification of DEGs, DEFRGs, and DEPRGs

The empirical Bayesian approach of the limma R package (http://www.bioconductor.org/packages/release/bioc/html/limma.html)39 was used to normalize GSE55235 and GSE77298 after standardization and setting |log2 FC| > 0.5 and P<0.05 to obtain DEGs of synovial tissues from RA in the 2 datasets. The ggplot2 package was used to plot the volcano map for both datasets.40 DEFRGs were created by intersecting RA-DEGs and FRGs, and DEPRGs were created by intersecting RA-DEGs and PRGs. DEFRGs and DEPRGs in GSE55235 were expressed using a clustered heat map. The data were analyzed by three methods: Principal Component Analysis (PCA), Uniform Manifold Approximation and Projection (UMAP), and t-distributed Stochastic Neighbour Embedding (t-SNE) dimensionality reduction using unsupervised clustering methods of the stats package, UMAP package, and Rtsne package in R software.

GO, KEGG, and ClueGO Enrichment Analyses

The “clusterProfiler” R package was used to run GO and KEGG enrichment analyses for DEFRGs and DEPRGs to investigate the found genes’ route and function (P < 0.05 and FDR < 0.25).41 ClueGO, a commonly used enrichment tool in Cytoscape software, employs the kappa method to visualize functional pathways in gene sets, enabling the construction of interactive gene network maps and the analysis of the potential functions of DEFRGs and DEPRGs.

PPI Network and CytoHubba Analysis

The STRING database (https://string-db.org/)42 was utilized to develop a PPI network, and the Cytoscape platform43 was used to construct and view it. The MCODE plug-in was applied for the PPI network module analysis. The DEFRGs and DEPRGs were computed using the Degree algorithm in the CytoHubba tool to find the top 10 high-scoring genes and visualize them by Cytoscape mapping.

Identification of Hub Genes and Machine Learning Validation

The top 10 genes of each of the DEFRGs and DEPRGs were taken to intersect to get the hub genes, and the expression of the hub genes was verified in GSE55235 using the joy plot and violin plot. In the ANN training model, the expression data of hub genes were transformed into “gene scores” using the min-max normalization method, and the neural network was constructed by the “NeuralNetTools” and “neuralnet” packages in R software.44 The ANN consists of an input layer (gene scores of candidate genes), a hidden layer (creating 5 hidden nodes and using modified linear units as activation functions), and an output layer (2 nodes: HC and RA). The gene weights in the hidden layer are multiplied by the corresponding gene expression to obtain the sum of the disease classification scores. Constructing ANN models for hub genes to determine phenotypic accuracy.

Construction of microRNA/TF-Hub Genes Regulatory Network

MicroRNA/TF-hub genes networks were constructed utilizing NetworkAnalyst (https://www.networkanalyst.ca/).45 The computational background of microRNA for hub genes was performed using the miRTarBase v8.0 and TarBase v8.0 databases, while the ENCODE and JASPAR databases were used for TF. The output is processed and plotted using Cytoscape for first-order and minimum networks.

Pathway Enrichment Analysis of Hub Genes

Metascape (https://metascape.org/gp/index.html#/main/step1)46 is a gene function enrichment platform that brings together 40 databases and sets P < 0.01, a minimum count of 3, and an enrichment factor > 1.5 to group genes into clusters to capture potential relationships between clusters further. GeneMANIA (http://www.genemania.org)47 integrates multiple databases from BioGRID and Pathway Commons and can predict the functional pathways by which hub genes interconnect with the most relevant genes. SPEED2 (https://speed2.sys-bio.net/)48 was able to probe the upstream relevant signaling pathways of hub genes; the genetic background was from human cell line pathway perturbation experiments, and the crosstalk of hub genes with 16 classical signaling pathways was assessed based on the magnitude of the P-value.

Validation of Hub Genes and Identification of Key Gene

The key gene was obtained by intersecting DEGs (GSE77298), DEFRGs, and DEPRGs, and the expression of the key gene in the validation set is shown in a box plot. ROC curves were used to validate diagnostic efficacy and verify the sensitivity and specificity of key genes, with AUC > 0.70 being diagnostically significant.

In vitro Experiments and Treatments

Immortalization of Synovial Fibroblasts in Human Rheumatic Joints (iCell-008a) was purchased from iCell Bioscience Inc., Shanghai. The iCell primary fibroblast culture system (priMed-iCell-003) was used as the medium for culturing FLS in vitro. The culture flasks were placed in a humidified incubator at 37°C, 5% CO2 for static incubation, and the medium was changed every 2 days until the cells reached the desired adherent growth state. We took logarithmic growth phase cells, adjusted the cell density, inoculated them into 6-well cell culture plates, and cultured them in a normal incubator with 5% CO2 for 24 h. The cells were treated in the following groupings. Control group: cultivation with complete medium; Ferroptosis group: pre-treatment with 1 μm RSL3 (MCE, HY-100218A) was added for 8 h, and incubation was continued for 24 h. Pyroptosis group: incubation was performed with 4 μm Polyphyllin VI (MCE, HY-N0816) for 24 h. Ferroptosis+pyroptosis group: After adding 1 μm RSL3 for 8 h of pre-treatment, 4 μm Polyphyllin VI was added for 24 h of incubation. Ferroptosis+quercetin group: after adding 1 μm RSL3 pre-treatment for 8h, continue to add 100 μm quercetin (MCE, HY-18085) and incubate for 24h. Pyroptosis+quercetin group: incubation was performed with 4 μm Polyphyllin VI and 100 μm quercetin for 24 h. Ferroptosis+pyroptosis+quercetin group: After adding 1 μm RSL3 for 8 h of pre-treatment, 4 μm Polyphyllin VI and 100 μm quercetin were added for 24 h of incubation. To ensure the reliability of the results, we tested the Control group and Ferroptosis+pyroptosis group in corresponding concurrent periods according to the time of phenotype formation in the Ferroptosis group and Pyroptosis group. In the quercetin-related studies, we re-experimented with the seven groups mentioned above: Control, Ferroptosis, Ferroptosis+quercetin, Pyroptosis, Pyroptosis+quercetin, Ferroptosis+pyroptosis and Ferroptosis+pyroptosis+quercetin groups.

Reverse Transcription-Quantitative PCR (RT-qPCR)

(1) RNA extraction. The four groups were rinsed with 1 mL, 4°C of pre-cooled PBS solution, added 1 mL of TRIpure Total RNA Extraction Reagent (EP013) and a series of other treatments, and then added 10 μL of RNase-Free Water to dissolve the RNA fully. (2) First-strand cDNA synthesis. The M-MLV Reverse Transcriptase kit (ELK Biotechnology, EQ002) was used. (3) RT-qPCR. Primers H-GAPDH and H-caspase-8 are shown in Table 2. A StepOne™ Real-Time PCR instrument was used with the QuFast SYBR Green PCR Master Mix kit. A total reaction volume of 10 µL was prepared prior to amplification, which included 5.0 μL 2× Master Mix, 1.0 μL Primer Working Solution (2.5 μm), 1.0 μL Template, 2.0 μL ddH2O, and 1.0 μL Rox. After pre-denaturation at 95°C for 1 min, 40 thermal cycles were performed, and melting curves were recorded. Quantification was performed using the 2−ΔΔCT method.49

Table 2 Primer Sequences Used for Reverse Transcription-Quantitative PCR

Analysis of Immune Cell Infiltration

The CIBERSORT algorithm (http://CIBERSORT.stanford.edu/)50 is a linear support vector regression-based method for assessing the composition and number of immune cells in the RA and HC groups. The ggplot2 package was used to illustrate the distribution of LM22, while the corrplot package was used to exhibit the correlation of LM22. Person correlation coefficients were used to analyze the relationship between key genes and immune cell abundance.

Quercetin-Related Targets

We used The Comparative Toxicogenomics Database (CTD, http://ctdbase.org/),51 the TCMSP database (https://old.tcmsp-e.com/tcmsp.php),52 the GeneCards database, and the CHEMBL database (https://www.ebi.ac.uk/chembl/)53 to find potentially relevant genes for quercetin (Supplementary Material has been uploaded to GitHub). The quercetin and key genes were further analyzed using molecular docking and molecular dynamics simulations.

Molecular Docking

The crystal structure of quercetin was downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/)54 as a docking ligand, and the crystal structure of the key gene-expressed protein was downloaded from the PDB database (https://www.rcsb.org/)55 as a docking acceptor. The Maestro software from the Schrödinger 2022–3 software package was used for this molecular docking. The charge state and the energy minimum conformation of quercetin in the physiological environment were calculated by the LigPrep module. The target proteins were optimized using the Protein Preparation Wizard in Maestro, and the eutectic complexes and the complex systems of target molecules and proteins were finely prepared using the protein-preparation workflow module. The lattice files defining the docking pockets were generated by the Receptor Grid Generation module in Maestro, and the center of mass of the eutectic ligand of the crystal structure was selected as the center of the active pocket and calculated to generate the lattice files. Finally, molecular docking was performed by the SP precision algorithm of the Glide docking module. The complexes obtained by docking were used as initial conformations for subsequent molecular dynamics simulations.

Molecular Dynamics Simulation

MDS of quercetin-target protein complexes was performed using the Desmond program. First, each complex system was dissolved separately in an SPC solvent box with a period boundary size of 10*10*10Å using the System Builder module to form the complex-solvent system. Second, adding a certain amount of Na+/Cl− balances the system charge. Finally, OPLS4 was chosen to simulate the force field. After completing the setup of the simulation system, we used the Molecular Dynamics module for formal simulations. Before the simulation, the simulation temperature was set to 300 K, and the pressure to 1.01325 bar, and a 100 ns NPT simulation was performed (energy was recorded every 1.0 ps, orbital atomic coordinates were monitored every 20 ps, and a total of 5000 frames of conformation were output). The behavior of quercetin and target proteins was analyzed using the Simulation Interaction Diagram tool in Desmond software, and root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) calculations throughout the simulation were output.

The flowchart of this study is shown in Figure 1, and all raw data has been uploaded to GitHub, which is HTTPS (https://github.com/zheng5862/zheng5862.git).

Figure 1 Schematic diagram of the workflow of this study.

Results Processing of Datasets and Identification of DEGs

The GSE55235 dataset includes 2230 DEGs, comprising 1279 upregulated genes and 951 downregulated genes. The GSE77298 dataset contains 499 DEGs, comprising 233 upregulated genes and 266 downregulated genes. The DEGs of the 2 datasets were presented using volcano plots (Figure 2). The RA-DEGs intersected with 592 FRGs to obtain 89 DEFRGs (Figure 3A). The RA-DEGs intersected with 254 PRGs to obtain 53 DEPRGs (Figure 3B). In the DEFRGs of GSE55235, 27 genes were upregulated, and 62 genes were downregulated. In the DEPRGs of GSE55235, 35 genes were upregulated, and 18 were downregulated. Using the clustering heat map, the distribution of up and down-regulated genes for 89 DEFRGs and 53 DEPRGs was readily visible (Figure 3C and D). PCA analysis of DEFRGs in GSE55235 showed PC1 (58.19%) and PC2 (7.02%), UMAP and tSNE plots showed significant classification of RA and control groups; PCA analysis of DEPRGs in GSE55235 showed PC1 (58.42%) and PC2 (8.27%), UMAP and tSNE plots showed significant categorization of RA and control groups (Figure 4A and B).

Figure 2 Identification of RA-DEGs. Yellow dots indicate genes expressed in synovial tissues of RA and HC groups that were not significantly different (P > 0.05); red triangles represent genes upregulated in RA (P < 0.05); and blue triangles represent genes down-regulated in RA (P < 0.05).

Figure 3 Screening of DEFRGs and DEPRGs. (A) Venn diagram of RA-DEGs of GSE55235 with ferroptosis-related genes; (B) Venn diagram of RA-DEGs of GSE55235 with pyroptosis-related genes; (C and D) Aggregated heat map distribution of 89 DEFRGs and 53 DEFRGs in GSE55235. The RA group is in red, the HC group is in blue, and the expression profiles from green to red indicate elevated gene expression, with green being down-regulated genes and red being upregulated genes.

Figure 4 (A and B) shows the PCA, UMAP, and tSNE 3 data downscaling analyses for DEFRGs and DEPRGs, respectively.

Functional Enrichment Analyses of DEFRGs and DEPRGs

In the GO analysis of DEFRGs, the biological process class (BP) was enriched in response to oxidative stress (Figure 5A); the cell formation class (CC) was enriched in a vesicle (Figure 5B); and the molecular function class (MF) was enriched in enzyme binding, oxidoreductase activity (Figure 5C). In the GO analysis of DEPRGs, BP was enriched in response to stress, immune system processes, and cytokine production (Figure 5D); CC was enriched in inflammasome complex, NLRP3, AIM2 inflammasome complex, and death-inducing signaling complex (Figure 5E); and MF was enriched in enzyme binding, cytokine receptor binding, death receptor binding, and cytokine-type endopeptidase regulator involved in apoptotic activities (Figure 5F). In the KEGG analysis, DEFRGs were enriched in ferroptosis, autophagy, mitophagy, apoptosis, and the FoxO signaling pathway (Figure 5G); DEPRGs were enriched in the NOD-like receptor signaling pathway, the IL-17 signaling pathway, the NF-kB signaling pathway, and the Toll-like receptor signaling pathway (Figure 5H). From the ClueGO analysis, it can be seen that DEFRGs were enriched in rheumatoid arthritis, ferroptosis, autophagy, mitophagy, the TNF signaling pathway, and the IL-17 signaling pathway (Figure 5I); DEPRGs were enriched in programmed cell death, pyroptosis, inflammasomes, NLR signaling pathways, the NLRP3 inflammasome, the AIM2 inflammasome, NF-κB activation through the FADD/RIP-1 pathway mediated by caspase-8 and −10, caspase activation via death receptors in the presence of ligand, caspase activation via extrinsic apoptotic signaling pathway and interleukin signaling (Figure 5J).

Figure 5 Pathway enrichment analysis of DEFRGs and DEPRGs. (A-C) show BP, CC, and MF analyses in GO for DEFRGs, respectively; (D-F) show BP, CC, and MF analyses in GO for DEPRGs, respectively; (G and H) shows KEGG analyses for DEFRGs and DEPRGs, respectively; and (I and J) below are ClueGO analyses for DEFRGs and DEPRGs, respectively.

Identification of Hub Genes

The PPI network of DEFRGs has 88 nodes and 356 edges (Figure 6A). CytoHubba identified the top 10 genes, including JUN, MYC, PTEN, EGFR, VEGFA, MAPK8, CASP8, MAPK1, PTGS2, and ATM (Figure 6B, Table 3). 52 nodes and 332 edges made up the PPI network of DEPRGs (Figure 6C). CytoHubba identified the top 10 genes as IL1B, CASP8, JUN, MYD88, CXCL8, CASP1, PTGS2, NLRP3, CASP3, and IL18 (Figure 6D, Table 4). These 2 sets of top 10 genes intersect as hub genes (CASP8, PTGS2, and JUN) (Figure 7A). We visualized the expression distribution and expression level of hub genes by joy plot and split-face violin plot (Figure 7B and C). CASP8 was substantially more expressed in RA than in the HC group, but the opposite was true for JUN and PTGS2 (P < 0.05). Our validation of hub genes in GSE55235 using ANN found that the model test was excellent (Figure 7D).

Table 3 Information of the 10 Hub Genes (DEFRGs)

Table 4 Information of the 10 Hub Genes (DEPRGs)

Figure 6 Screening of candidate hub genes. (A) The PPI network of DEFRGs has 88 nodes and 356 edges. (B) DEFRGs screened the top 10 candidate hub genes by Degree algorithm. (C) The PPI network of DEPRGs has 52 nodes and 332 edges. (D) DEPRGs screened the top 10 candidate hub genes by Degree algorithm.

Figure 7 (A) Screening of hub genes; (B) Expression distribution of hub genes in the RA group of GSE55235; (C) Differential expression of hub genes in the RA and HC groups of GSE55235 is shown by split-face violin plot, with red representing the RA group and blue representing the HC group; (D) ANN model showing the input layer, hidden layer, and output layer separately.

MicroRNA/TF-Hub Genes Regulatory Network

The miRNA-hub genes network graph obtained from the miRTarBase database includes 101 nodes, 103 edges, and 3 seeds, and the streamlined minimum network includes 5 nodes, 5 edges, and 3 seeds (Figure 8A), which shows that hsa-mir-26b-5p is closely related to CASP8, PTGS2, and JUN. The miRNA-hub genes network graph obtained from the TarBase database includes 184 nodes, 230 edges, and 3 seeds, and the minimum network includes 4 nodes, 3 edges, and 3 seeds (Figure 8B), which shows that hsa-mir-16-5p is tightly linked to CASP8, PTGS2, and JUN. The network diagram of TF-hub genes obtained from the ENCODE database includes 90 nodes, 93 edges, and 3 seeds, and the minimum network includes 5 nodes, 4 edges, and 3 seeds (Figure 8C), which shows that SP3 is closely related to CASP8, JUN, and SAP30 is closely related to PTGS2, JUN. The network diagram of TF-hub genes obtained from the JASPAR database includes 30 nodes, 32 edges, and 3 seeds, and the minimum network includes 4 nodes, 3 edges, and 3 seeds (Figure 8D), which shows that GATA2 is associated with CASP8, PTGS2. Thus, 2 microRNAs (hsa-mir-26b-5p and hsa-mir-16-5p) and 3 TFs (SP3, SAP30, and GATA2) may play important roles in the expression of hub genes.

Figure 8 MicroRNA/TF regulatory network of hub genes. (A and B) MicroRNA: first-order network obtained from hub genes in the miRTarBase and TarBase databases, respectively, with the minimum network in the lower right corner; (C and D) TF: first-order network obtained from hub genes in the ENCODE and JASPAR databases, respectively, with the minimum network in the lower right corner.

Pathway Enrichment Analysis of Hub Genes

Metascape analysis showed that hub genes were enriched in the IL-17 signaling pathway and visualized using network diagrams (Figure 9A). In GeneMANIA, the pathways between hub genes and the 20 most relevant genes were enriched in the regulation of apoptotic signaling pathway, positive regulation of cysteine-type endopeptidase activity, programmed necrotic cell death, mitochondrial outer membrane permeabilization involved in programmed cell death, IκB kinase/ NF-κB signaling, positive regulation of proteolysis, and regulation of cytokine-mediated signaling pathway (Figure 9B). Positive regulation of hub genes by Hypoxia and the MAPK+PI3K signaling pathways can be seen in SPEED2 (Figure 9C).

Figure 9 Enrichment analysis of hub genes. (A) Pathway and process richness analysis by Metascape and the corresponding network diagram; (B) Gene-gene interaction network of hub genes and the 20 most relevant genes detected by GeneMANIA database; (C) SPEED2 analysis of hub genes, with gray to red in the scale indicating elevated correlation.

Identification of Key Gene

CASP8 was identified as a key gene by intersecting RA-DEGs of GSE77298 with the top 10 genes of DEFRGs/DEPRGs (Figure 10A). CASP8 was significantly higher expressed in the RA group compared to healthy controls in the GSE77298 dataset (P < 0.01) (Figure 10B). With excellent sensitivity and specificity, the AUC values of CASP8 in the GSE55235 and GSE77298 datasets were 0.94 (95% CI 0.84–1.00) and 0.86 (95% CI 0.69–1.00), respectively (Figure 10C).

Figure 10 Identification and validation of key genes. (A) Further screening of key gene by GSE77298; (B) Expression of CASP8 in the validation set of GSE77298; red indicates the RA group, blue indicates the HC group; (C) Diagnostic validity of CASP8 was verified by ROC curve in GSE55235 and GSE77298; (D-H) In vitro experiments and RT-qPCR validation. In these bar plots, (D) the horizontal coordinates are control group, ferroptosis group, and ferroptosis+pyroptosis group; (E) the horizontal coordinates are control group, pyroptosis group, and ferroptosis+pyroptosis group; (F) the horizontal coordinates are control group, ferroptosis group and ferroptosis+quercetin group; (G) the horizontal coordinates are control group, pyroptosis group and pyroptosis+quercetin group; (H) the horizontal coordinates are control group, ferroptosis+pyroptosis group and ferroptosis+pyroptosis+quercetin group. The vertical coordinate is the relative mRNA levels of caspase-8/GAPDH.

In vitro Validation of Key Gene

PCR-Amplification curves were attached to Figures S1-S5. The results of Figure 10D showed that in the ferroptosis phenotype constructed by adding RSL3 to FLS, caspase-8 level was significantly higher than that in the control group (P < 0.0001); in the ferroptosis+pyroptosis phenotypes constructed by adding RSL3 and Polyphyllin VI to FLS, caspase-8 level was significantly higher than that in the control group (P < 0.0001) and ferroptosis group (P < 0.0001) (Figure S1). Figure 10E shows that the level of caspase-8 in the pyroptosis phenotype constructed by adding Polyphyllin VI to FLS was significantly higher than that in the control group (P < 0.0001); in the ferroptosis+pyroptosis phenotypes constructed by adding RSL3 and Polyphyllin VI to FLS, caspase-8 level was significantly higher than that in the control group (P < 0.0001) and pyroptosis group (P < 0.0001) (Figure S2). These results suggest that caspase-8 is significantly more expressed than normal FLS in both the ferroptosis and pyroptosis phenotypes present. When both ferroptosis and pyroptosis phenotypes co-occur, caspase-8 is also more highly expressed than the ferroptosis alone or pyroptosis alone phenotypes. The possible reason for this is that caspase-8 is the common hub center of ferroptosis and pyroptosis, and the superposition of the two phenotypes is positively correlated with the expression of caspase-8. Figure 10F shows that the caspase-8 level was significantly reduced in the ferroptosis treatment model constructed by adding RSL3 and quercetin to FLS compared to the ferroptosis group (P < 0.0001) (Figure S3). Figure 10G shows that caspase-8 levels were significantly reduced in the pyroptosis treatment model constructed by adding Polyphyllin VI and quercetin to FLS compared to the pyroptosis group (P = 0.0005) (Figure S4). Figure 10H shows that caspase-8 levels were also significantly reduced in the ferroptosis+pyroptosis treatment model constructed with the addition of RSL3, Polyphyllin VI, and quercetin to FLS compared to ferroptosis+pyroptosis group (P < 0.0001) (Figure S5). Thus, quercetin can reduce caspase-8 levels in ferroptosis, pyroptosis, and ferroptosis+pyroptosis phenotypes constructed in FLS. Quercetin is a potential therapy for RA via targeting caspase-8 through ferroptosis and pyroptosis.

Immune Infiltration Analysis

We used the CIBERSORT method to map the proportions of 22 immune cells in RA samples in GSE55235 (Figure 11A). Differences in immune cell infiltration between the RA and HC groups were analyzed using a box plot (Figure 11B). The results demonstrated that 4 kinds of immune cells were enriched in RA: plasma cells, CD8+ T cells, memory-activated CD4+ T cells, and follicular helper T cells (Tfh) (P < 0.05). The relationship between the degree of LM22 infiltration in RA samples was further investigated by correlation matrix analysis (Figure 11C), and it was found that all 4 types of immune cells (plasma cells, CD8+ T cells, memory-activated CD4+ T cells, and Tfh) were positively correlated in RA (P < 0.05). According to the analysis of the Person correlation coefficient, only memory-activated CD4+ T cells (R=0.85, P=0.0018) and Tfh (R=0.69, P=0.03) among the above 4 types of immune cells were found to be correlated and positively correlated with CASP8 (Figure 11D). Memory-activated CD4+ T cells and Tfh were significantly higher expressed in RA with statistically significant results (P < 0.05). These results imply that CASP8 may be implicated in the pathogenesis of RA via modulating immune cell infiltration and that memory-activated CD4+ T cells and Tfh may play a significant role in this regard.

Figure 11 Analysis of immune cell infiltration. (A) The plot of the percentage of LM22 in RA samples using the CIBERSORT algorithm, where the horizontal coordinate is the sample and the vertical coordinate is the percentage of immune cells, and the legend shows the LM22 species; (B) Box plots of immuno-infiltrating cells in the RA versus HC groups, with the RA group in red and the HC group in blue (*P < 0.05,**P < 0.01,***P < 0.001, ****P < 0.0001); (C) Correlation matrix analysis between immune cells in RA groups. Both horizontal and vertical coordinates are LM22, with positive correlations in red and negative correlations in blue (*P < 0.05,**P < 0.01,***P < 0.001, ****P < 0.0001); (D) Correlation analysis of CASP8 expression with memory-activated CD4+ T cells and Tfh.

Quercetin-Related Targets

Our results of extracting quercetin-related target genes from CTD, TCMSP, GeneCards, and CHEMBL databases and taking intersections with hub genes were only CASP8 and were in the top 34, 66, 76, and 122 of the four databases, respectively (Table S4).

Molecular Docking

Molecular docking predicts the conformation of quercetin within the caspase-8 target binding site and assesses binding affinity. We obtained quercetin’s 2D and 3D structures and the protein structure caspase-8 (PDB ID: 4JJ7) (Figure 12A-C). As seen in the interaction diagram of quercetin with caspase-8 (Figure 12D), quercetin binds in the internal cavity of the caspase-8 protein. We can see from the detailed diagram on the right that quercetin interacts with SER-316, ARG-260, ASP-455, ARG-258, ARG-413, GLY-318, HIS-317, GLN-358 on the caspase-8 in a hydrogen-bonding-related manner and with TYR-412, TYR-365 on the caspase-8 in a hydrophobic manner. These interactions underlie the maintenance of quercetin and caspase-8 binding. In addition, docking scoring showed −7.17 kcal/mol, implying that quercetin binds closely to caspase-8 protein.

Figure 12 (A) 2D structure of quercetin molecule; (B) 3D structure of quercetin molecule; (C) 3D structure of caspase-8 protein; (D) Binding mode of quercetin and caspase-8 protein obtained based on molecular docking. The left figure shows the overall view, and the right figure shows the partial view, in which the Orange stick is the small molecule, the blue cartoon is the protein, the blue solid line indicates the hydrogen bonding interaction, and the gray dashed line indicates the hydrophobic interaction.

Molecular Dynamics Simulation (100 Ns)

Considering that the complexes obtained by molecular docking can only observe the binding static information, MDS can further explore the dynamic nature, stability, and protein structural flexibility of quercetin-caspase-8 complexes, a deepening experiment in molecular docking technology. From the RMSD (Figure 13A), it can be observed that caspase-8 gradually converges and stabilizes within 15 ns. It is worth mentioning that the RMSD of the protein after convergence is around 3.5–4 Å. The high RMSD is due to the excessive loop regions in the structure of caspase-8 itself. The RMSD of quercetin was within 2 Å, indicating that quercetin did not break away from the active pocket of the protein, and the fluctuation of quercetin implied that there was a conformational change of the small molecule in the active pocket, but the fluctuation was within 1.5 Å. Therefore, we believe this fluctuation did not affect the overall binding conformation. From Figure 13B, it can be observed that the overall RMSF of caspase-8 is within 2 Å, implying that the protein is low in flexibility, and although the protein has a high number of loop regions (white areas), quercetin binds to the loop-enriched region, and this binding region is also within 1.5 Å. This result implies that the proteins in the binding region of the small molecule did not undergo large conformational changes and could exhibit dynamic stable binding. The intrinsic RMSF fluctuation of quercetin can be seen in Figure 13C. It is noteworthy that the RMSF of the dihydroxybenzene region is larger, implying that the flexibility of this part is greater, probably due to the presence of large variations in the dihedral angles 12-10-11-19 during the molecular dynamics simulations.

Figure 13 (A) RMSD of quercetin and caspase-8 during the simulation; (B) RMSF of caspase-8, with light red indicating the helix region, blue indicating the beta-sheet region, white indicating the loop region, and dark green indicating the small-molecule binding site (amino acid numbers in the horizontal coordinate plus 222 is the actual amino acid number); (C) RMSF plot of heavy atoms on Quercetin with a distribution of dihedral angles 12–10-11-19.

Small molecules and proteins are motile in physiological environments, and we used MDS to simulate this dynamic pattern and observe amino acid interactions in the active pockets of quercetin and caspase-8 during movement. High-frequency water-bridging interactions between quercetin and GLY-318 and CYS-360 on caspase-8, as well as hydrogen-bonding interactions with ARG-260, ARG-413, TYR-365, and ASP-455, and cation-Pi interactions with ARG-413, can be observed from Figure 14A. In addition, we counted more detailed interaction forces and the amino acids with which the interactions occur. Quercetin interacts with SER-256, ARG-258, ASP-259, ARG-260, ASN-261, ASP-266, SER-316, HIS-317, GLY-318, GLN-358, ALA-359, CYS-360, GLY-362, ASP-363, ASN-364, TYR-365, GLN-366, VAL-410, SER-411, TYR-412, ARG-413, ASN-414, TRP-420, ASP-455, ASN-458 on caspase-8, forming around 8 stable interactions during the simulation, with a high frequency of interactions with amino acids ARG-258, TRP-420, ARG-260, TYR-365, ARG-413, and ASP-455 (above 0.4) (Figure 14B and C). Hydrogen bonding and water bridging are the main interaction forces between small molecules and amino acids in the pocket. These interactions suggest that quercetin and caspase-8 can bind dynamically and stably.

Figure 14 (A) 3D/2D plot of quercetin and caspase-8 interaction during MDS. Purple lines indicate hydrogen-bonding interactions and red lines indicate water-bridging interactions; (B) Histogram of the percentage of quercetin and caspase-8 interactions during MDS. The cyan color indicates hydrogen bonding, the blue color indicates water-bridging, and the lavender color indicates hydrophobicity; (C) Heat map analysis of quercetin and caspase-8 interactions during MDS.

We calculated the radius of gyration (rGyr), intramolecular HBs (intraHB), molecular surface area (MolSA), solvent accessible surface area (SASA), and polar surface area (PSA) of quercetin with time using kinetic simulation trajectories (Figure 15). The fluctuations in these properties are relatively stable, indicating that quercetin can maintain its basic physicochemical properties in the caspase-8 protein binding pocket and undergoes its active action with relatively consistent properties. In addition, the intraHB data show that small molecules can form an intramolecular hydrogen bonding interaction during the 100 ns simulation, contributing to the conformational stabilization of small molecules.

Figure 15 RGyr, intraHB, MolSA, SASA, PSA of quercetin/caspase-8 system over time during MDS.

Discussion The Relationship Between RA and Ferroptosis

Ferroptosis has a crucial regulatory function in autoimmune and chronic inflammatory disorders.19,20 Luo et al found that RSL3 (ferroptosis activator) induced synoviocyte ferroptosis and exacerbated synovitis.56 First, high iron ions inhibit the growth and activity of osteoblasts57 and activate the differentiation of osteoclasts, resulting in bone destruction.58 High iron ions also promote the aberrant proliferation of synovial cells and the establishment and progression of pannus and vascular synovitis in the synovium,59 consistent with iron deposits in RA synovial tissue.60 Second, ROS, as a product of oxidative stress, plays an important role in the pathogenesis of RA,61 and there is a significant positive association between ROS levels and RA severity.18 The ROS-activated NF-κB signaling pathway produces TNF-α, further accelerating RA via the p38/JNK pathway.62 ROS can also activate NLRP3, crosstalking with the pyroptosis pathway. Third, ferroptosis exacerbates inflammatory responses by stimulating the immune system to release inflammatory cytokines.29 Iron accumulation, ferroptosis activation, and their associated signaling pathways cause a vicious cycle of hemorrhage-inflammation-re-hemorrhage-re-inflammation in the synovium, destroying bone and cartilage in the joints of RA patients (Figure 16).

Figure 16 Diagram of the proposed mechanism for the involvement of caspase-8 in the ferroptosis, pyroptosis pathway in RA.

The Relationship Between RA and Pyroptosis

It is now known that pyroptosis primarily works through the NLRP3 inflammasome, which consists of NLRP3, ASC, caspase-1, caspase-8, GSDMD, IL-1β, and IL-18 effector molecules, to exacerbate the progression of RA.63,64 The inflammasome is an intracellular multiprotein complex mediating immune responses to cellular damage,65 where important pyroptosis effector molecules include NLRP3 inflammasomes, which are divided into two phases of initiation and activation.66 The ROS-NF-κB-NLRP3 axis is the initiation phase.67 The activation phase is the NLRP3/ASC/procaspase-1/caspase-1 axis. NLRP3 binds to the adaptor ASC and recruits procaspase-1 to form the classical inflammasome, ultimately driving caspase-1 activation.68 The following effect phase is divided into two independent routes. First, the cleavage of GSDMD by activated caspase-169 results in a GSDMD-N structural domain, triggering a breach in the cell membrane and initiating pyroptosis.70,71 Second, the cleavage of pro-IL-1β and pro-IL-18 by activated caspase-1 releases IL-1β and IL-18 to initiate inflammatory responses (Figure 16).72

Caspase-8 is a Biomarker of Ferroptosis and Pyroptosis in RA

Our study identified a key gene, CASP8, associated with ferroptosis and pyroptosis in RA, which may be an important biomarker and a therapeutic target. Caspase is a class of cysteine-dependent endonucleases that cleave and activate themselves and other caspases by hydrolyzing their targets.73 Effector caspases must be cleaved by initiating caspases to activate pyroptosis and release inflammatory factors.74,75 Caspase-1 was originally assumed to be involved in pyroptosis, but as research progressed, other caspases (such as caspase-3, −4, −5, and −8) were also found to be involved.76–78 Within the family of caspases, caspase-8 functions as both an initiator and an apical activator, comprising two N-terminal death effector domains (DED), a large (p18) and a small (p10) subunit at its C-terminus, and was originally recognized for its role in apoptosis.79 A subsequent genetic investigation showed that caspase-8 is a molecular switch for apoptosis, necroptosis, and pyroptosis.80 Several investigations also demonstrated that while apoptosis and necroptosis are suppressed, NLRP3/caspase-1 axis-mediated pyroptosis continues to be induced by caspase-8.81,82 Caspase-8 can function as both an initiating and an effector caspase and play essential roles in distinct signaling pathways (such as cell death, inflammasome, and inflammatory cytokine production).83–85

First, caspase-8 was an upstream regulator of classical and non-classical NLRP3 inflammasome initiation and activation.86,87 Second, caspase-8 acts as an initiating caspase that activates caspase-1 and its downstream signaling pathways.88,89 Third, caspase-8 can substitute caspase-1 as an effector caspase to induce pyroptosis and IL-1β production.78 In the inflammasome pathway, the PYD region of the ASC binds to the DED of caspase-8 and thus recruits and activates caspase-8,89 which completes the pyroptosis process by cleaving the substrate at the same amino acid D276 as caspase-1 in the GSDMD.24 Fourth, GSDMD-driven pyroptosis may be triggered by caspase-8 via regulating direct channels of GSDMD.90 Fifth, caspase-8 enhances IL-1β production by directly regulating the NF-κB signaling pathway.91 Non-apoptotic caspase-8 is involved in pyroptosis and inflammatory cytokine release and is a facilitator and reinforcer of the NLRP3 inflammasome signaling pathway.78,92 Thus, caspase-8 is essential for pyroptosis-related pathways (Figure 16).

Another investigation has demonstrated that caspase-8 activation occurs concurrently with ROS generation.93 During ferroptosis, ROS generated by excessive intracellular Fe2+ and H2O2 via the Fenton reaction activates the NF-κB pathway to generate TNF-α, which in turn can trigger caspase-8 via the p38/JNK pathway and further accumulate ROS.12,17,62,94 In the NLRP3 inflammasome pathway activated via the TLR4/MyD88/NF-κB axis, caspase-8 is an upstream regulator of NF-κB and promotes TNF-α production, indirectly promoting ROS synthesis and ferroptosis (Figure 16).95,96 Caspase-8 promotes ferroptosis and pyroptosis, acts as an upstream positive regulator of the production and secretion of inflammatory cytokines such as TNF-α, IL-1β, and IL-18, and plays an important role in the positive feedback loop of RA progression. Therefore, Caspase-8 is a potentially important therapeutic target.

Cells in RA

(1) Fibroblast-like synoviocytes. In RA, FLS exhibited a distinct aggressive tendency.97 FLS can induce chronic synovitis, pannus development, bone degradation, and cartilage damage through the production of inflammatory factors, proangiogenic factors, and matrix-degrading enzymes.98 FLS plays an important role in the development and progression of RA, in which TNF-α is a key trigger in the proinflammatory process of FLS.99 The NF-κB/TNF-α pathway generates ROS in response to NOX activation100 and induces caspase-8 activation,101 then ROS enhances TNF-mediated activation of the NF-κB pathway,102 generating a positive loop of accumulating ROS that drives ferroptosis and pyroptosis. (2) Immune cells. The immune-infiltrating environment of the synovium consists mainly of immune cells such as T cells, B cells, and macrophages.103 Studies have shown that memory-activated CD4+ T cells accumulate in RA joints and are important immune cells promoting inflammation.104,105 Tfh can further assist B cells in promoting antibody production.106 Th17 cells can secrete IL-17