Exploring differential gene expression and biomarker potential in systemic lupus erythematosus: a retrospective study

Research participants and sample collection

This study is a retrospective study. A total of 45 SLE patients and 40 healthy controls (HCs) were collected from Nanfang Hospital of Southern Medical University from January 2024 to December 2024. All SLE patients met the 2012 ACR SLE disease classification criteria (Petri et al., 2012; Tan et al., 1982). All the SLE patients with concurrent infection, acute concurrent illness, and use of probiotics or antibiotics within 1 month before admission were excluded. The collected HCs also needed to have no history of known autoimmune diseases and be of comparable age to SLE patients. All the participants were collected two mL EDTA-K₂ blood sample.

Ethics approval was granted by the Ethics Committee of Nanfang Hospital, under the surveillance of ethical number: NFEC-2025-026. All the methods used were in accordance with the approved guidelines. Written informed consent was required from all patients and healthy volunteers in the study.

The SLE patients and HCs in this study were comparable in age and gender (Table 1). Average age of the SLE and HC groups was 38.09 ± 11.97 and 42.20 ± 10.34, respectively (p = 0.099).

Datasets

Two raw datasets, GSE13887 and GSE10325, were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE13887 was an analysis of CD3-positive T cells isolated from 10 SLE patients and nine HCs. GSE10325 was an analysis of freshly isolated lymphocytes (CD4⁺ T cells and CD19⁺ B cells) and CD33⁺ myeloid subsets from the blood of 28 HCs and 39 SLE patients.

Identification of differentially expressed genes

Differentially expressed genes (DEGs) between SLE and HCs were identified using the limma package (version 3.64.1) in R (version 4.2.1). The GSE13887 and GSE10325 datasets were first converted into expression matrices, appropriately grouped, and normalized. Normalization and differential expression analysis were performed using the limma package. Genes with an adjusted p-value (false discovery rate, FDR) < 0.05 and an absolute log2 fold change (|log2FC|) ≥ 0.5 were considered statistically significant. DEGs were classified as upregulated or downregulated based on whether their log2FC values were above 0.5 or below −0.5, respectively. An online Venn diagram tool was employed to identify overlapping DEGs between SLE and HC. The results were visualized using ggplot2 (version 3.4.4) for volcano plots (with thresholds set at p ≤ 0.05 and |log2FC| ≥ 1) and box plots to assess normalization across datasets. Heatmaps were generated using the ComplexHeatmap package (version 2.13.1), while unique and shared DEGs were further illustrated using the VennDiagram package (version 1.7.3) (Gu, Eils & Schlesner, 2016).

Gene Ontology term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

The R language packages clusterProfiler (version 4.4.4) and ggplot2 were used to perform Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on common DEGs (co-DEGs) and visualize the results, with a threshold of p value < 0.05, and bar charts, bubble charts, and chord plots were drawn (Yu et al. , 2012).

Construction of protein-protein interaction network and identification of hub genes

The PPI network of DEGs that may play an important role in the progression of SLE was constructed using the STRING online database (https://www.string-db.org/). The protein–protein interaction (PPI) network was performed using Cytoscape (version 3.10.2) (http://www.cytoscape.org/) to better visualize the interaction information. Hub genes with connectivity ≥10 were selected using the CytoHubba (version 0.1).

Hub gene validation

Receiver operating characteristic (ROC) analysis was performed on the data using the R language package pROC (version 1.18.0). The results were visualized using ggplot2 to plot the ROC curve of the hub gene. The area under the curve (AUC) of the ROC curve corresponding to the hub gene was used to evaluate the ability to distinguish between SLE and HC. The expression profiles of the hub genes in the two data sets were used as variables for principal component analysis (PCA) to obtain PC1 and PC2, and ggplot2 was used to visualize the results.

CD3⁺ T cell isolation from SLE and HC periphal blood samples

Human primary PBMC were isolated from periphral blood with EDTA-K2 by density gradient centrifugation using Ficoll. 1 X 10⁶ CD3⁺ T cells were enriched by positive selection with an CD3⁺ T cells isolation kit (130-097-043, Miltenyi Biotec, Bergisch Gladbach, Germany). The purity of T cell were analyzed by Flow cytometry (CantoII, BD Bioscience, Franklin Lakes, NJ, USA) staining with anti-CD3 -PE antibody (562310, clone: HIT3A, BD Bioscience, Franklin Lakes, NJ, USA).

Real-time quantitative polymerase chain reaction

The sorted CD3⁺ T cells were lyzed with TRIzol reagent (Cat. #G3013 Servicebio, Wuhan, China). After chloroform substitute (Cat. #G3014 Servicebio, Wuhan, China) isolation, RNA were percipitated by isopropanol by centrifugation at 4 °C, 15,000 g, 15 min. After that, RNA were rinsed with 2 times 75% ethanol. The cDNA were synthesized with SweScript RT I First Strand cDNA Synthesis Kit (Cat. # G331-1, Servicebio, Wuhan, China) by Random Hexamer Primer (Cat. #G331-4, Servicebio, Wuhan, China). The RNA template were added one µg.

The qPCR were conducted by using Servicebio^® 2 ×SYBR Green qPCR Master Mix (Cat. #G3326-1, Servicebio, Wuhan, China). The reation volume were used 20 µL in each wells, cDNA templated were diluted in 10 times with DEPC water and 2 µL, forward and reverse primers were added 0.4 µL. The qPCR were performed with Roche Cobas 4800 PCR machine (Roche, Basel, Switzerland). This experiment used a two-step reaction procedure: pre-denaturation: 95 °C, 2–5 min; cyclic amplification (35–45 cycles): 95 °C denaturation 5–15 s; 60–68 °C annealing/extension 20–60 s (simultaneous fluorescence collection, SYBR Green collection at the end point, TaqMan collection in real time). Melting curve (SYBR Green): gradually increase the temperature from 65 °C to 95 °C and monitor the fluorescence. Relative gene expression was calculated using the 2^−ΔΔCt method, and the results were normalized using GAPDH. Primer sequences are shown in Table S1.

Flow cytometry

The FCER1A and RGS1 protein expression levels of CD3⁺ T cells in plasma collected from HCs and SLE patients and the levels of inflammatory factors in plasma were detected by flow cytometry. Samples were analyzed using a BD FACS Fortassa flow cytometer (BD Biosciences, USA). Data acquisition was performed using BD FACSDiva software, and at least 10,000 events were recorded for each sample. Polyclonal rabbit anti-human RGS1 antibody (Cat# LS-C162570-400; LSBio, Seattle, WA, USA) and FITC-conjugated goat-anti-rabbit secondary antibodies were used for RGS1 analysis (Jiang et al., 2024). Polyclonal rabbit anti-human RGS1 antibody (Cat# LS-C162570-400; LSBio, Seattle, WA, USA) and FITC-conjugated goat-anti-rabbit secondary antibodies were used for RGS1 analysis. A rabbit polyclonal FCER1 alpha Monoclonal Antibody (Cat# 16-5899-82, Functional Grade, eBioscience, San Diego, CA, USA) was used for FCER1A analysis (Greer et al., 2014).

Statistics

All data processing and analysis were performed in GraphPad Prism 9.0. To compare two groups of variables, the Mann–Whitney U test was applied. Fisher’s exact test was performed to analyze the statistical significance between two variable data sets. P < 0.05 was considered statistically significant. Statistical significance is reported as follows: p < 0.05(*), p < 0.01(**), p < 0.001(***), p < 0.0001(****), ns(not significant).

Identification of DEGs

Utilizing bioinformatics analysis through R software (version 4.2.1), a total of 1,912 DEGs were identified in the GSE13887 dataset, while 52 DEGs were detected in the GSE10325 dataset. Subsequent hierarchical clustering analysis was performed on these DEGs, prioritizing the top 20 upregulated and downregulated genes from GSE13887 (Table S3) and GSE10325 (Table S4), respectively. Volcano plots and heatmaps were generated to visualize the clustering results for GSE13887 (Figs. 1A, 1C) and GSE10325 (Figs. 1B, 1D). The heatmaps distinctly illustrated a high degree of consistency in sample clustering.

Following data normalization and comparative evaluation across datasets, box plots were constructed to confirm the normalization quality for GSE13887 (Fig. 1E) and GSE10325 (Fig. 1F). The results demonstrated that the data distribution across both datasets conformed to established quality standards, affirming the robustness and cross-comparability of the microarray data. A Venn diagram (Fig. 1G) was subsequently generated to illustrate the overlap of DEGs between the two datasets, identifying eight shared DEGs.

Enrichment analysis of DEGs

GO and KEGG pathway enrichment analyses were conducted to elucidate the functional roles of the eight common differentially expressed genes (co-DEGs). GO terms were categorized into biological process (BP), cellular component (CC), and molecular function (MF), with the top five most significantly enriched terms in each category selected based on the lowest p-values (Table S5). These results were visualized via bar charts (Fig. 2A) and bubble plots (Fig. S1A). The co-DEGs were predominantly associated with processes such as cell killing, leukocyte-mediated cytotoxicity, natural killer cell-mediated immunity, and natural killer cell-mediated cytotoxicity.

Subsequently, the five most significantly enriched KEGG pathways were identified (Table S6) and represented through bar plots (Fig. 2B) and bubble charts (Fig. S1B). The co-DEGs were chiefly involved in pathways related to transcriptional dysregulation in cancer, apoptosis, thyroid cancer, African trypanosomiasis, and asthma.

PPI network and hub genes

A chord diagram was generated to illustrate the enrichment of genes within the top five GO categories across BP, CC, and MF, identifying seven key genes: GZMB, LAG3, APOL1, CXCL13, RGS1, FCER1A, and DPEP2. Notably, RGS1, CXCL13, and GZMB were enriched in at least two GO categories (Figs. 2C–2E). The DEGs were uploaded to the STRING database to construct the PPI network (Fig. S2). Using the MCODE plugin, a tightly interconnected protein cluster consisting of CXCL13, GZMB, LAG3, and FCER1A was identified (Fig. 3A).

Overall, seven hub genes—GZMB, LAG3, APOL1, CXCL13, RGS1, FCER1A, and DPEP2—were selected based on both the PPI network and chord diagram analyses. Volcano plot assessments of their expression levels in the GSE13887 and GSE10325 datasets revealed that only FCER1A and RGS1 exhibited consistent expression trends across both datasets, with FCER1A being downregulated and RGS1 upregulated (Fig. S3). These two genes were further analyzed using the STRING database to construct an additional PPI network (Fig. S4), visualized in Cytoscape. MCODE analysis identified a highly interconnected protein cluster comprising ten genes (Figs. 3B, 3C).

Prognostic value of hub genes

ROC curves were generated to evaluate the diagnostic efficacy of FCER1A and RGS1, both individually and in combination, using their expression levels in the GSE13887 and GSE10325 datasets (Fig. 4). Both genes demonstrated strong diagnostic potential, with area under the curve (AUC) values exceeding 0.7. The combined analysis of these two genes exhibited superior diagnostic performance.

PCA was performed to further validate their discriminatory power. The PCA plot, derived from dimensionality reduction of gene expression data in both HCs and SLE patients, demonstrated clear separation between the two groups based on the expression profiles of FCER1A and RGS1 (Fig. 5). The first two principal components (PC1 and PC2) accounted for 100% of the explained variance. SLE samples clustered predominantly on the negative axis of PC1, while HCs were localized on the positive axis, indicating significant differences in gene expression patterns between the groups and highlighting the robust discriminatory capability of these two genes.

Verification by other datasets and clinical samples

To further substantiate these findings, the expression of FCER1A and RGS1 was validated using the GSE61635 dataset, which corroborated their diagnostic efficacy. The expression trends observed in GSE61635 were consistent with those from GSE13887, GSE10325, and clinical samples, with FCER1A downregulated and RGS1 upregulated in SLE patients relative to HCs (Fig. 6).

Peripheral blood samples were collected, and CD3⁺ T cells were isolated using magnetic bead separation. RNA extraction followed by RT-qPCR revealed significant downregulation of FCER1A (p < 0.0001) and upregulation of RGS1 (p < 0.0001) in SLE patient samples compared to HCs (Fig. 7A). These findings were corroborated by flow cytometry, which demonstrated a marked reduction in FCER1A protein levels (p < 0.01) and a significant increase in RGS1 protein expression (p < 0.001) in the plasma of SLE patients (Figs. 7B, 7C).

Additionally, flow cytometric analysis of inflammatory cytokines revealed elevated plasma levels of IL-6 and TNF-α (p < 0.01), alongside a significant reduction in IL-10 (p < 0.001) in SLE patients compared to HCs (Fig. 7D). These results further support the potential of FCER1A and RGS1 as robust biomarkers for the diagnosis and pathophysiological understanding of SLE.

Conclusion

In this study, we employed integrative bioinformatics analysis and clinical validation to identify and verify FCER1A and RGS1 as potential biomarkers for SLE. These genes exhibited consistent differential expression patterns across multiple independent datasets and patient-derived clinical specimens, with FCER1A markedly downregulated and RGS1 significantly upregulated in individuals with SLE. Functional enrichment and PPI analyses implicated these genes in critical immunological processes, including immune regulation, cell-mediated cytotoxicity, and apoptosis—hallmarks of SLE pathogenesis. Moreover, dysregulated levels of pro- and anti-inflammatory cytokines in patient plasma further substantiate the immunopathological relevance of these findings. Collectively, our results provide a foundation for the potential application of FCER1A and RGS1 in SLE diagnosis and suggest their promise as targets for future therapeutic strategies. Further mechanistic studies are warranted to elucidate their functional roles in disease progression.