Single Cell and Transcriptomic Analysis of Regulatory Mechanisms of Key Genes in Systemic Lupus Erythematosus

Introduction

Systemic lupus erythematosus (SLE) is a chronic autoimmune disease primarily caused by immune system dysfunction, leading to damage in multiple organs such as the skin, kidneys, joints, nervous system, and cardiovascular system.¹ The global incidence and prevalence of SLE are steadily increasing, with incidence rates ranging from 1.4 to 15.13 per 100,000 individuals and an annual addition of approximately 400,000 new cases. The prevalence ranges from 15.87 to 108.92 per 100,000 individuals, posing a serious threat to human health.² SLE is characterized by heterogeneous clinical manifestations and complex molecular mechanisms, with abnormal immune regulation at the cellular level playing a central role in its pathogenesis.³ Among immune cells, B cells are particularly crucial in driving disease development. Dysregulated B cell function contributes to SLE through mechanisms such as the production of autoantibodies, excessive secretion of inflammatory cytokines, and antigen presentation, all of which promote chronic inflammation and tissue damage.⁴ Despite the availability of various treatment options, such as immunosuppressants, antimalarials, and biologics, the cure rate for SLE remains low, and treatment outcomes vary among individuals.⁵ Thus, the identification of effective biomarkers is essential for early diagnosis and the prediction of treatment efficacy.

In recent years, bioinformatics has emerged as a powerful approach for analyzing large-scale biological data to uncover disease-related mechanisms and identify key molecular drivers. By leveraging transcriptomic, genomic, and proteomic datasets, bioinformatics enables systematic screening for disease-associated biomarkers and therapeutic targets. Additionally, single-cell RNA sequencing (scRNA-seq) has revolutionized our ability to profile cellular heterogeneity, revealing previously unrecognized immune cell subpopulations and gene expression patterns at single-cell resolution. Numerous studies have applied these methods to investigate immune mechanisms in SLE. However, despite these advances, a comprehensive screening of key genes associated with B cell dysfunction in SLE remains limited.

Therefore, our study integrates single-cell transcriptomics with conventional bulk transcriptomic data to enable high-resolution mapping of immune cell alterations in SLE. This integrative approach allows us to identify key pathogenic genes more robustly, particularly those associated with dysfunctional B cells, thus advancing our understanding of SLE pathogenesis and providing a theoretical foundation for novel therapeutic strategies.

Materials and Methods

Data Sources

This study was approved by the Ethics Committee of the Second Affiliated Hospital of Kunming Medical University, and written informed consent was obtained from all participants. A total of 6 SLE patients and 6 age- and sex-matched healthy controls were recruited. Peripheral blood (8 mL) was collected from each participant, and peripheral blood mononuclear cells (PBMCs) were isolated for single-cell transcriptomic sequencing analysis. The GSE82221 dataset was sourced from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), which included transcriptome data (GPL10558) and methylation data (GPL13534) of peripheral blood from 25 control (CTL) samples and 30 systemic lupus erythematosus (SLE) samples.

Differential Gene and Immune Cell Analysis

Differentially expressed genes (bulk DEGs) between the SLE and CTL groups were identified in the GSE82221 dataset using the limma package,⁶ with thresholds set at P<0.05 and |log2FC|>1. The results were visualized using volcano and heat maps generated by the ggplot2 (v 3.3.2)⁷ and heatmap (v 4.1.0)⁸ packages, respectively, to display DEG distributions. Additionally, to investigate immune function differences between the SLE and CTL groups, differential immune cells (DICs) were identified using single-sample gene set enrichment analysis (ssGSEA) within the GSVA package (v 1.70–3),⁹ based on 28 immune cell gene sets.⁶

Single-Cell Sequencing Data Analysis

Single-cell data were filtered using the R package Seurat (v 4.0.6)¹⁰ by excluding cells with ≥5% mitochondrial gene expression, low-quality cells (fewer than 200 or more than 6000 genes), and genes detected in 3 or fewer cells. The remaining cells were integrated using the Harmony method, followed by gene variance analysis on batch-effect-corrected cells. A set of 1500 highly variable genes was selected for subsequent analysis. Principal component analysis (PCA) was then performed (p < 0.05), and sample clustering was conducted using the UMAP algorithm based on the principal components (PCs). Cluster annotations were performed using Azimuth (https://github.com/satijalab/azimuth),¹¹ with annotation scores ranging from 0 to 5 (0–2: low annotation confidence; 2–3.75: possible annotation issues; 3.75–5: high annotation confidence). Finally, the cell proportion differences between the SLE and CTL groups (cell type count/total cell count per patient) were calculated and analyzed. Cells exhibiting the same differential trends as those observed in the CTLs from section 2.2 were defined as key cells.

Key Gene Selection and Gene Set Enrichment Analysis (GSEA)

Differentially expressed genes (single-cell DEGs) in key cells were identified using the FindMarkers function in Seurat (SLE vs CTL), with criteria of p < 0.05 and |log2FC| > 0.5 (v 4.0.6).¹⁰ The bulk DEGs and single-cell DEGs were taken to intersect to get candidate genes, in which the genes with area under curve (AUC) values of receiver operating characteristic (ROC) curves greater than 0.7 were designated as key genes. In order to explore whether there were interactions between key genes, protein-protein interaction (PPI) network was constructed for them using GeneMANIA website (http://genemania.org/). The role played by key genes in SLE was explored by GSEA, where the expression of all genes in the single-cell dataset with key genes was first calculated and ranked, and then enriched using the clusterProfiler package (P value < 0.05) (v 6.1.0).¹² Additionally, methylation analysis of key genes was performed in the GSE82218 dataset using the ChAMP package (v 2.80).¹³

Construction of Regulatory Networks

The transcription factors (TFs) targeting key genes were predicted using NetworkAnalyst database (https://www.networkanalyst.ca/). In addition, the miRWalk (http://mirwalk.umm.uni-heidelberg.de/) and starBase (https://ngdc.cncb.ac.cn/databasecommons/database/id/169) databases were used to predict miRNAs targeting key genes. The miRNAs that were predicted in both databases were used as the co-rmiRNAs in this study. Next, starBase (clipExpNum > 0) and miRNet (https://www.mirnet.ca/) databases were utilized to predict target co-lncRNAs targeting the co-miRNAs. Finally, the key gene-TF and lncRNA-miRNA-mRNA network were constructed by Cytoscape (v 3.8.2).¹⁴

Analysis of Single Cell Expression Level

To explore the expression differences of key genes between different B cell subpopulations as well as between SLE and CTL, differential analyses were performed using kruskal.test. At the same time, monocle (v 2.22.0)¹⁵ was used for analyzing the cell differentiation trajectory of differentially expressed cells. Finally, the cellular time-series analysis was combined with the expression of key genes and presented.

Statistical Analysis

Differences were analyzed via the Wilcox. Test and kruskal.test between two groups. P value < 0.05 was considered statistically significant.

Results

Identification of Bulk DEGs and DICs

A total of 2041 bulk DEGs were identified between the SLE and CTL groups, including 847 up-regulated genes and 1194 down-regulated genes (Figure 1A and B). Meanwhile, the results of ssGSEA algorithm showed that the scores of 14 immune cells such as activated CD4 T cell, CD56 dim natural killer cell, central memory CD8 T cell, effector memory CD4 T cell and memory B cell were significantly different between the SLE and CTL groups (Figure 1C).

Figure 1 Identification of differentially expressed genes and immune cells between SLE and normal groups. (A) Heatmap of gene expression density illustrating the expression levels of 2041 differentially expressed genes (DEGs) between the SLE and control groups. (B) Volcano plot depicting the distribution of DEGs between the SLE and control groups. (C) Violin plot showing the significant differences in 14 immune cell types between the SLE and control groups, as calculated using the ssGSEA algorithm(ns: not significant, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

Identification of Key Cell Clusters at the Single-Cell Level

To further analyze key cell clusters at the single-cell level, single-cell analysis was performed after cell filtering (Figure S1A). A total of 1500 highly variable genes were identified (Figure S1B). PCA plots revealed that single cell samples in SLE and CTL groups were dispersed, indicating that there was no significant batch effect of sequencing results (Figure S1C). Subsequently, 35 principal components (p < 0.05) were selected for further analysis (Figure S1D–E). Clustering results revealed 24 cell clusters (Figure 2A and S1F), annotated into 30 different cell types, such as CD4 TCM cells, CD8 TEM cells, and NK cells (Figure 2B). Analysis of the cellular proportion differences between subpopulations showed significant differences in ASDC cells, memory B cells, and CD8. Proliferating cells (Figure 2C). The trend of memory B cell differences was consistent with the differential immune cells at the transcriptional level, being higher in the CTL group and lower in the SLE group (Figure 2D). Therefore, memory B cells were selected as the key cell type for further analysis.

Figure 2 Clustering and Cellular Proportions in SLE and Control Groups. (A) UMAP plot displaying the clustering results of single-cell samples from the SLE and control groups, with cells classified into 24 clusters. (B) UMAP plot annotated with 30 distinct cell types. (C) Stacked bar plot showing the proportions of different immune cell subsets in each sample from the SLE and control groups. (D) Comparison of the proportions of different immune cell subsets between the SLE and control groups. (ns: not significant; *p < 0.05; **p < 0.01).

Key Genes Involved in Virus-Related Pathways

In memory B cells, 19 differentially expressed genes (single-cell DEGs) were identified between SLE and CTL groups. The top four single-cell DEGs included IFI44L, IGHV5-51, IFI27, and MX1 (Figure 3). Six candidate genes (MX1, IFI44L, ISG15, OAS1, IFI6, and IFI27) were selected by intersecting bulk DEGs and single-cell DEGs (Figure 4A). ROC analysis revealed that these genes had diagnostic value for SLE (AUC > 0.7) (Figure 4B). Thus, they were defined as key genes for this study. PPI network analysis showed that the key genes are involved in regulating antiviral responses. For example, the OAS1 gene co-expressed with OAS2, participating in pathways such as antiviral response regulation, negative regulation of viral processes, and viral genome replication (Figure 4C). Additionally, methylation analysis indicated that IFI27, IFI44L, IFI6, ISG15, MX1, and OAS1 had 6, 4, 12, 6, 14, and 1 methylation sites, respectively, with significant differences in methylation levels (Figure 4D). For instance, significant methylation differences were observed in cg10778971, cg12177001, cg14352715, and cg23893332 in IFI27. To further explore the potential roles of MX1, IFI44L, ISG15, OAS1, IFI6, and IFI27, single-gene GSEA analysis was performed. KEGG analysis revealed that most key genes were significantly enriched in pathways such as cytoplasmic translation, viral defense response, and symbiont defense response (Figure 5 and S2).

Figure 3 UMAP plot illustrating the top four DEGs in memory B cells.

Figure 4 Identification and analysis of key genes in SLE. (A) Venn diagram showing the overlap of differentially expressed genes (DEGs) between memory B cells and the whole transcriptome dataset. (B) Receiver operating characteristic (ROC) curve analysis demonstrating the diagnostic value of these key genes in SLE. (C) Protein-protein interaction (PPI) network of key genes. (D) Methylation analysis showing the number of significantly methylated sites for each key gene (ns: not significant; *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001; red text indicates statistically significant differences).

Figure 5 Single-gene GSEA analysis displaying KEGG enrichment results for MX1, IFI6, IFI27, IFI44L, ISG15 and OAS1.

ceRNA Regulation of Genes and miRNA-mRNA-TF Networks

To explore the regulatory mechanisms of MX1, IFI44L, ISG15, OAS1, IFI6, and IFI27, a key gene-TF network was constructed, consisting of 41 nodes and 48 edges (Figure 6A). In this network, BRCA1 was found to potentially influence the expression of both IFI6 and MX1. Moreover, an mRNA-miRNA-lncRNA network comprising 53 nodes and 73 edges was constructed, consisting of four miRNAs (hsa-miR-485-5p, hsa-miR-28-5p, hsa-miR-185-5p, and hsa-miR-326) (Figure 6B), two mRNAs (IFI44L and OAS1), and 47 lncRNAs (such as GAS5, FGD5-AS1) (Figure 6C). This network included pathways such as hsa-miR-485-5p-OAS1-MEG3 and IFI44L-hsa-miR-326-LINC00894 (Figure 6D).

Figure 6 Construction of Regulatory Networks Involving ceRNA, miRNA-mRNA, and Transcription Factor (TF) Interactions. (A) Key gene-transcription factor (TF) network illustrating the regulatory relationships between key genes and transcription factors. (B) Venn diagram showing the overlapping miRNA-mRNA interactions identified by MirWalk and StarBase. (C) Venn diagram displaying the overlapping lncRNAs identified by MirWalk and StarBase. (D) mRNA-miRNA-lncRNA regulatory network.

Expression of Key Genes in Cell Pseudotime Analysis

MX1, ISG15, and OAS1 were expressed in memory B cells, naïve B cells, and intermediate B cells, while IFI6 was highly expressed in intermediate B cells, and IFI44L was highly expressed in naïve B cells and intermediate B cells (Figure 7A). Additionally, the expression levels of the six key genes were higher in the B cell subpopulations of the SLE group compared to the CTL group (Figure 7B). Further exploration of biomarker expression during B cell differentiation revealed three distinct states (Figure 8A and B), with state 1 primarily comprising memory B cells, and naïve B cells and intermediate B cells predominantly in state 2 (Figure 8C and D). Ultimately, combining the cell pseudotime analysis with key gene expression showed that IFI44L and MX1 initially increased and then decreased over pseudotime (Figure 8E).

Figure 7 Expression of Key Genes in Different B Cell Subpopulations. (A) Gene expression analysis across various B cell subpopulations. (B) Comparison of key gene expression levels in B cell subpopulations between the SLE and control groups (****p < 0.0001).

Figure 8 Pseudotime analysis of B cell subpopulations and the expression of key genes in these subpopulations. (A and B) The pseudotime trajectory plot shows the differentiation process of B cell subpopulations, with three distinct states observed in the B cell subpopulations. (C) Pseudotime analysis combined with cell state classification, displaying the distribution of B cell subpopulations across different states. (D) Pseudotime analysis in the control and SLE groups. (E) Temporal dynamics of key genes such as IFI27, IFI44L, IFI6, ISG15, MX1, and OAS1 in each state based on pseudotime analysis.

Discussion

Systemic lupus erythematosus is a globally prevalent autoimmune disease, and elucidating its underlying mechanisms is critical for developing novel diagnostic and therapeutic strategies. In this study, single-cell and transcriptomic analyses identified six key genes associated with SLE (MX1, IFI44L, ISG15, OAS1, IFI6, IFI27), which were significantly differentially expressed in immune cells, particularly memory B cells, from SLE patients. Further analysis suggested that these genes contribute to immune dysregulation in SLE through pathways such as viral response and T-cell receptor signaling. Additionally, we explored their epigenetic regulation, ceRNA network interactions, and associations with transcription factors and miRNAs, offering new insights into potential therapeutic targets for SLE.

IFI44L, ISG15, IFI6, and IFI27, as interferon-stimulated genes, play crucial roles in immune responses, with methylation modifications potentially influencing SLE susceptibility. OAS1 is involved in antiviral immunity, while MX1 is upregulated following viral infection, and its high expression may contribute to immune dysregulation.¹⁶ Methylation of MX1 and IFI44L has been proposed as biomarkers for SLE,^17,18 with MX1 demonstrating strong diagnostic potential in proteomics and single-cell sequencing studies.¹⁹ Additionally, MX1 and IFI16 have been identified as potential biomarkers for lupus nephritis (LN), suggesting applications in LN diagnosis and treatment.^20,21 IFI44L promoter methylation has been established as a blood biomarker for SLE,²² influencing SLE pathogenesis through the type I interferon signaling pathway.²³ Our study confirmed that IFI44L is highly expressed in SLE patients, particularly in memory B cells, further supporting its critical role in SLE pathogenesis. Furthermore, IFI44L, IFI27, and OAS1 are closely associated with SLE progression,²⁴ with IFI44L methylation having been suggested as a highly sensitive diagnostic biomarker.²⁵ Free-form ISG15 may exacerbate disease severity by enhancing pro-inflammatory activity,^26,27 and its high expression in B cells of SLE patients may contribute to immune response imbalance.

Significant differences exist between naive and memory B cells in SLE patients,²⁷ with abnormal memory B cells disrupting immune homeostasis and being associated with SLE relapse.^28,29 Despite exhibiting reduced responsiveness to B cell receptor (BCR) signaling, SLE memory B cells retain responsiveness to toll-like receptor and type I interferon signaling,³⁰ which may reduce the efficacy of B cell-targeted therapies such as Bruton’s Tyrosine Kinase (BTK) inhibitors and anti-CD22 antibodies.³¹ These findings suggest that aberrant memory B cells may play a role in SLE pathogenesis and could potentially serve as therapeutic targets. However, further functional studies and clinical investigations are needed to confirm their pathogenic role and assess the feasibility of such therapeutic strategies.³²

Functional enrichment analysis revealed that these key genes are primarily involved in antiviral defense, cell cycle regulation, and immune-related pathways, aligning with the immunopathological characteristics of SLE. MX1, IFI44L, OAS1, and ISG15 are closely associated with antiviral responses, with OAS1 implicated in Epstein–Barr virus infection and MX1 involved in NK cell-mediated cytotoxicity.^33,34 Additionally, IFI6 and IFI44L participate in T cell receptor signaling and Th17 cell differentiation, both of which are critical in SLE pathogenesis.^35,36 Methylation analysis identified multiple significant methylation sites in MX1, IFI44L, ISG15, and OAS1, highlighting the essential role of epigenetic modifications in SLE immune dysfunction.³⁷

Pseudotime analysis revealed the temporal expression dynamics of key genes in B-cell subsets. MX1, ISG15, and OAS1 were linked to immune dysregulation in SLE, while IFI6 and IFI44L played roles in B-cell differentiation and regulation. These genes were expressed in memory, naïve and intermediate B cells, suggesting viral influence on B-cell immunity. Notably, IFI6 was enriched in intermediate B cells, whereas IFI44L was upregulated in naïve and intermediate B cells, highlighting their functional significance. Their widespread upregulation indicates dysregulated B-cell activation in SLE. Pseudotime analysis also identified three B-cell differentiation states, with IFI44L and MX1 peaking early before declining, underscoring their roles in B-cell development. These findings link gene dynamics to immune regulation, cell cycle control, and viral responses, providing insights into SLE pathogenesis and potential therapies.

Additionally, this study constructed a key gene–transcription factor regulatory network, revealing that BRCA1 may regulate IFI6 and MX1, supporting its involvement in immune responses.³⁸ Furthermore, mRNA–miRNA–lncRNA network analysis identified hsa-miR-485-5p as a potential regulator of OAS1, while the IFI44L–hsa-miR-326–LINC00894 axis may be implicated in SLE immune mechanisms. Further experimental validation is needed to explore the clinical significance of these regulatory interactions However, it should be noted that these conclusions are primarily based on computational predictions, which represent a limitation of this study. Future research should incorporate functional validation through knockdown or overexpression experiments, as well as chromatin immunoprecipitation (ChIP) assays, to confirm the biological relevance of these regulatory interactions.

Conclusions

This study systematically identified key genes implicated in SLE pathogenesis (MX1, IFI44L, ISG15, OAS1, IFI6, and IFI27) and their potential regulatory mechanisms in memory B cells using single-cell sequencing and transcriptomic analysis. These genes are likely to contribute to SLE pathogenesis by modulating antiviral responses, immune signaling pathways, and epigenetic regulation, thereby providing novel insights into SLE diagnosis and targeted therapy. These findings suggest that the identified genes may serve as potential biomarkers for early screening and as molecular targets for precision therapies in clinical settings. However, this study is limited by its reliance on computational predictions, the possibility of batch effects, and sample heterogeneity. Future studies should incorporate functional validation approaches, including gene knockout or overexpression, CRISPR-Cas9-based editing, and clinical validation in larger SLE cohorts, to further verify the biological significance and clinical applicability of these findings.