Single-cell analysis with childhood and adult systemic lupus erythematosus

Abstract

Patients with systemic lupus erythematosus (SLE), a heterogeneous and chronic autoimmune disease, exhibit unique changes in the complex composition and transcriptional signatures of peripheral blood mononuclear cells (PBMCs). While the mechanism of pathogenesis for both childhood-onset SLE (cSLE) and adult-onset SLE (aSLE) remains unclear, cSLE patients are considered more unpredictable and dangerous than aSLE patients. In this study, we analysed single-cell RNA sequencing data (scRNA-seq) to profile the PBMC clusters of cSLE/aSLE patients and matched healthy donors and compared the PBMC composition and transcriptional variations between the two groups. Our analysis revealed that the PBMC composition and transcriptional variations in cSLE patients were similar to those in aSLE patients. Comparative single-cell transcriptome analysis between healthy donors and SLE patients revealed IFITM3, ISG15, IFI16 and LY6E as potential therapeutic targets for both aSLE and cSLE patients. Additionally, we observed that the percentage of pre-B cells (CD34^-) was increased in cSLE patients, while the percentage of neutrophil cells was upregulated in aSLE patients. Notably, we found decreased expression of TPM2 in cSLE patients, and similarly, TMEM150B, IQSEC2, CHN2, LRP8 and USP46 were significantly downregulated in neutrophil cells from aSLE patients. Overall, our study highlights the differences in complex PBMC composition and transcriptional profiles between cSLE and aSLE patients, providing potential biomarkers that could aid in diagnosing SLE.

Keywords:

• Single-cell analysis
• childhood SLE
• adult SLE
• PBMC composition
• transcription

Introduction

Systemic lupus erythematosus (SLE) is a chronic autoimmune disease characterised by its heterogeneity and its association with cardiac and cardiovascular (CV) comorbidities, which arise due to chronic inflammation [1]. Approximately 20% of SLE cases are diagnosed during childhood (cSLE), while the majority occur in adulthood (aSLE) [2]. The disease course, especially during early adolescence, is unpredictable and can lead to progressive organ damage [3]. The pathogenesis of SLE is complex and involves a multifactorial interaction between genetic and epigenetic factors that contribute to the development of autoimmune responses. Several studies indicate that interferon (IFN)-α plays a significant role in driving SLE [4] and can serve as a potential biomarker for SLE patients [5]. MicroRNAs (miRNAs) also play a crucial role in the pathogenesis of SLE [6,7]. The characteristics of cSLE are similar but different from that of aSLE, exhibiting a higher prevalence of fever, lymphadenopathy, lupus nephritis and anti-dsDNA antibodies [8–10]. Notably, cSLE tends to be more severe and aggressive compared to aSLE, exhibiting manifestations such as malar rash, haemolytic anemias, lupus nephritis, and anti-double stranded DNA [11,12]. The 10-year survival rate of SLE patients treated with corticosteroids and immunosuppressive drugs is nearly 90%. However, the evident toxicity associated with most medications and the diverse manifestations of SLE pose challenges to the design of clinical trials. Therefore, it is essential to further investigate the underlying mechanisms of pathogenesis in both cSLE and aSLE.

Single-cell RNA sequencing (scRNA-seq) not only reveals cell-to-cell variations in gene expression at the single nucleotide level but also enables the identification of gene expression heterogeneity between individual cells [13]. This technology offers an unprecedented resolution for studying the pathogenesis of SLE, a disease characterised by heterogeneous and phenotypically diverse cell populations. However, there have been limited studies comparing the differences between cSLE and aSLE using scRNA-seq. Recent research indicates that the clinical stratification and transcriptional profiles of aSLE resemble those of cSLE [14]. These findings establish a framework for stratifying SLE and highlight specific cell subpopulations that could serve as potential therapeutic targets. Nevertheless, the novel biomarkers specific to cSLE or aSLE have yet to be defined. In this study, we conducted scRNA-seq analysis on peripheral blood mononuclear cells (PBMCs) obtained from cSLE, aSLE and healthy donors and identified distinct cell subpopulations and specific biomarkers, which may have the potential to guide new strategies for SLE therapy.

Results

scRNA-seq reveals the characteristics of PBMC composition and transcription in cSLE patients

To investigate the differences in PBMC composition and transcriptional profiles between cSLE patients and healthy children donors, we obtained 10X scRNA-seq matrix data of PBMCs from cSLE (n = 33) and healthy children donors (n = 11) from the Gene Expression Omnibus (GEO) database (accession number, GSE135779 [14]. We conducted quality control on the scRNA-seq data and retained cells that met specific criteria: (1) detectable gene expression in at least 100 cells, (2) gene copy numbers between 400 and 3,000 detected in cells, and (3) gene copy numbers in mitochondria less than 5% in cells. Consequently, we retained approximately 162,841 cells from cSLE patients and 71,397 cells from healthy donors for further analysis (Figure 1A). Using unsupervised clustering analysis followed by two-dimensional t-distributed stochastic neighbour embedding (t-SNE) [15], we identified 34 molecularly distinct clusters (Figure 1B-C). These clusters were annotated to eight cell types in PBMCs: B cells [16], monocytes (CD16^-), monocytes (CD16⁺), NK cells, CD4⁺ T cells, CD8⁺ T cells, pre-B cells (CD34^-), and common myeloid progenitor (CMP) (Figure 1D-E). Furthermore, when compared with the matched healthy donors, we observed that the percentage of monocytes (CD16^-), monocytes (CD16⁺), and CD8⁺ T cells were upregulated in the cSLE group, while the percentage of NK cells and pre-B cells (CD34^-) was decreased (Figure 1F-G). To further explore the molecular changes associated with cSLE, we analysed the differentially expressed genes (DEGs) in different PBMC clusters between cSLE patients and healthy children donors. The results revealed significant upregulation of IFI27 (Interferon Alpha Inducible Protein 27; p value = 0, Wilcoxon test), ISG15 (p value = 0, Wilcoxon test), and IFITM3 (Interferon Induced Transmembrane Protein 3; p value = 0, Wilcoxon test) in monocytes (CD16^-) of cSLE patients. Additionally, monocytes (CD16⁺) from cSLE patients exhibited significant upregulation of IFI27 (p value = 1.45E-96, Wilcoxon test), ISG15 (p value = 1.29E-258, Wilcoxon test), and IFI6 (Interferon-alpha inducible protein 6; p value = 1.66E-239, Wilcoxon test). Moreover, CD8⁺ T cells of cSLE patients demonstrated significant upregulation of IFI27 (p value = 0, Wilcoxon test), HBA2 (Hemoglobin Subunit Alpha 2; p value = 6.51E-199, Wilcoxon test), and HBB (Hemoglobin Subunit Beta; p value = 3.41E-98, Wilcoxon test) (Figure 1H). These findings indicate changes in the composition of PBMCs in cSLE patients and highlight the transcriptional alterations in these specific genes, which may be associated with the development of cSLE.

Figure 1. Single-cell RNA sequencing reveals the composition and transcriptional characteristics of peripheral blood mononuclear cells (PBMCs) in childhood-onset systemic lupus erythematosus (cSLE) patients. (A-B) Two-dimensional projection of the sample distribution (A) and identification of 34 subclusters (B) across a total of 234,238 PBMCs from 44 individuals, including 33 cSLE patients (case) and 11 healthy donors (control). (C) t-SNE analysis identifies eight distinct PBMC subsets, including B cells, monocytes (CD16^-), monocytes (CD16⁺), natural killer (NK) cells, CD4⁺ T cells, CD8⁺ T cells, pre-B cells (CD34^-), and common myeloid progenitor (CMP). (D) Annotation of the eight clusters identified in C. The dot plot shows the expression values of selected genes (x-axis) across each cluster (y-axis). dot size represents the percentage of cells expressing the marker gene, and colour intensity indicates the mean expression within expressing cells. (E-F) Bar plots illustrating the abundance of cells within each of the eight clusters (B cells: 25,451 cells, monocytes (CD16^-): 30,413 cells, monocytes (CD16⁺): 5,215 cells, NK cells: 14,810 cells, CD4⁺ T cells: 134,874 cells, CD8⁺ T cells: 22,632 cells, pre-B cells (CD34^-): 673 cells, CMP: 170 cells). (G) Differential expression analysis highlights the differentially expressed genes (DEGs) across the eight clusters.

scRNA-seq reveals the characteristics of PBMC composition and transcription in aSLE patients

To investigate the differences in PBMC composition and transcription between aSLE patients and matched healthy donors, we obtained 10X scRNA-seq matrix data of PBMCs from aSLE (n = 10) and healthy adult donors (n = 7) from the GEO database with accession numbers GSE135779 and GSE142016 [14,17]. Following a similar quality control process as before, we retained approximately 32,839 cells from aSLE patients and 36,749 cells from healthy donors for further analysis (Figure 2A). Applying unsupervised clustering analysis followed by two-dimensional t-SNE, we identified 30 molecularly distinct clusters based on gene expression values compared to all other cells (Figure 2B-C). After annotation, these clusters were assigned to seven cell types in PBMCs: B cells, monocytes (CD16^-), monocytes (CD16⁺), NK cells, CD4⁺ T cells, CD8⁺ T cells, and neutrophil cells (Figure 2D-E). Moreover, we observed an upregulation of monocytes (CD16^-), monocytes (CD16⁺), CD8⁺ T cells, and neutrophil cells, while the percentage of NK cells and B cells was lower in the aSLE group compared to the healthy donors’ group (Figure 2F-G). Further analysis of the differential expressed genes (DEGs) in PBMC clusters between aSLE patients and matched healthy donors revealed significant upregulation of HBB, HBA2, and IFITM3 (p value = 7.36E-287, p value = 7.42E-259, and p value = 0, Wilcoxon test) in monocytes (CD16^-) of aSLE patients. Monocytes (CD16⁺) of aSLE patients exhibited significant upregulation of IFI27, ISG15, and IFI6 (p value = 5.56E-61, p value = 2.67E-95, and p value = 1.49E-109, Wilcoxon test). Furthermore, CD8⁺ T cells of aSLE patients showed significant upregulation of IFI27, IFI44L (Interferon Induced Protein 44 Like), and HBB (p value = 1.32E-293, p value = 0, and p value = 0, Wilcoxon test), while neutrophil cells of aSLE patients displayed significant upregulation of IFITM2 (Interferon Induced Transmembrane Protein 2), S100A8 (S100 Calcium Binding Protein A8), and S100A9 (S100 Calcium Binding Protein A9; p value = 0.00065, p value = 0.00055, p value= 0.00056, Wilcoxon test) (Figure 2H). These findings indicate that alterations in PBMC composition and the associated transcriptional changes are closely linked to the development of aSLE.

Figure 2. Single-cell RNA sequencing reveals the composition and transcriptional characteristics of peripheral blood mononuclear cells (PBMCs) in adult-onset systemic lupus erythematosus (aSLE) patients. (A-B) Two-dimensional projection of the sample distribution (A) and identification of 30 subclusters (B) across a total of 69,588 PBMCs from 17 individuals, including 10 aSLE patients (case) and 7 adult healthy donors (control). (C) t-SNE analysis identifies seven distinct PBMC subsets, including B cells, monocytes (CD16^-), monocytes (CD16⁺), natural killer (NK) cells, CD4⁺ T cells, CD8⁺ T cells, and neutrophil cells. (D) Annotation of the seven clusters identified in C. The dot plot shows the expression values of selected genes (x-axis) across each cluster (y-axis). dot size represents the percentage of cells expressing the marker gene, and colour intensity indicates the mean expression within expressing cells. (E-F) Bar plots illustrating the abundance of cells within each of the seven clusters (B cells: 3,632 cells, monocytes (CD16^-): 9,356 cells, monocytes (CD16⁺): 2,571 cells, NK cells: 3,981 cells, CD4⁺ T cells: 39,908 cells, CD8⁺ T cells: 8,751 cells, and neutrophil cells: 1,389 cells). (G) Differential expression analysis shows the differentially expressed genes (DEGs) across the seven clusters.

The similarities and differences in molecular pathogenesis between cSLE and aSLE patients

Upon comparing the PBMC proportions between cSLE and aSLE patients based on the previous data (Figures 1D and 2D), we observed that both cSLE and aSLE patients exhibited similar clusters with the same cell types, such as B cells, monocytes (CD16^-), monocytes (CD16⁺), NK cells, CD4⁺ T cells, and CD8⁺ T cells. Moreover, we identified overlapping DEGs between cSLE and aSLE patients within these cell types (Figure 3A), and despite these similarities, there were also distinct differences in the distribution of PBMC clusters between cSLE and aSLE patients, including in monocytes (CD16⁺), monocytes (CD16^-), CD8⁺ T cells, NK cells, B cells, pre-B cells (CD34^-) and neutrophil cells (Figure 3).

Figure 3. Single-cell RNA sequencing reveals the overlapped differentially expressed genes (DEGs) in peripheral blood mononuclear cells (PBMCs) between childhood-onset (cSLE) and adult-onset (aSLE) systemic lupus erythematosus patients. Venn diagram illustrating the overlapped DEGs between cSLE patients and aSLE patients within the same PBMC cluster.

To gain further insights into the expression variation within PBMC subtypes between aSLE and cSLE patients, Gene Ontology (GO) enrichment analysis for the DEGs was performed in monocytes (CD16⁺), monocytes (CD16^-) and CD8⁺ T cells in both aSLE and cSLE patients compared to matched healthy donors (Figures 1G and 2G). The results revealed three upregulated networks in monocyte (CD16⁺) cells from both cSLE and aSLE patients: the interferon signalling pathway, cell death pathway, and response to virus pathway (Figure 4A-B). Notably, the cell death pathway exhibited significant upregulation in monocyte (CD16⁺) cells from cSLE patients compared to aSLE patients (Figure 4A-B, Table 1). Additionally, we visually presented the top five DEGs in monocytes (CD16⁺) between SLE patients and healthy donors and observed significant upregulation of IFITM3 (aSLE: p value = 1.31E-151, cSLE: p value = 1.16E-284, Wilcoxon test) and LY6E (Lymphocyte Antigen 6 Family Member E; aSLE: p value = 1.04E-122, cSLE: p value = 1.25E-252, Wilcoxon test) in monocyte (CD16⁺) cells of both aSLE and cSLE patients (Figure 4C-D). Previous studies suggested that LY6E could serve as a potential biomarker for SLE patients [18], and IFITM3, an antiviral factor known to inhibit viral replication through cholesterol homeostasis suppression [19], could serve as a novel diagnostic biomarker for cSLE. Furthermore, we observed significant upregulation of ISG15, EEF1A1 (Eukaryotic Translation Elongation Factor 1 Alpha 1) and IFI6 (cSLE: p value = 1.29E-258, p value = 4.75E-256, and p value = 1.66E-239, Wilcoxon test) in cSLE patients but not in aSLE patients (Figure 4C-D). In addition, we observed a significant increase in the expression levels of HLA-DQA2 (Major Histocompatibility Complex, Class II, DQ Alpha 2) [20], HBB and HBA2 (aSLE: p value = 7.78E-140, p value = 2.38E-126, and p value = 6.24E-132, Wilcoxon test), specifically in monocyte (CD16⁺) cells of aSLE patients compared to cSLE patients. Taken together, these results suggest that HLA-DQA2, HBB and HBA2 could be as new potential diagnostic biomarkers for aSLE.

Figure 4. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in monocyte (CD16⁺) cells between systemic lupus erythematosus (SLE) patients and matched healthy donors. (A-B) Gene Ontology (GO) enrichment analysis of DEGs in monocyte (CD16⁺) cells was conducted separately for adults (A) and children (B). (C-D) the top 5 DEGs in monocyte (CD16⁺) cells for both adults (C) and children (D) were further analysed using the Seurat R package.

Table 1. Functional enrichment results of differentially expressed genes between SLE patients and healthy donors in different cell types

term_name	adjusted_p_value
response to biotic stimulus	3.29E-21
immune system process	2.30E-18
regulation of viral process	5.94E-17
defense response	2.74E-16
innate immune response	3.14E-14
response to type I interferon	2.36E-10
response to interferon-beta	4.38E-08
antigen processing and presentation of peptide or polysaccharide antigen via MHC class II	1.37E-07
response to interferon-alpha	3.63E-07
response to interferon-gamma	5.11E-07
antigen processing and presentation of exogenous peptide antigen via MHC class II	1.80229E-06
mitochondrial ATP synthesis coupled electron transport	6.29785E-06
MHC protein complex assembly	1.18522E-05
peptide antigen assembly with MHC protein complex	1.18522E-05
leukocyte cell-cell adhesion	3.32679E-05
cell death	4.4851E-05
cell killing	5.89823E-05
regulation of lymphocyte activation	0.000189073
leukocyte mediated immunity	0.000193137
regulation of cell death	0.000656868

Significantly upregulated MHC-associated DEGs were observed in monocyte (CD16^-) cells from cSLE patients compared to aSLE patients (Figure 5A-B, Table 1). These MHC molecules are involved in various processes such as antigen processing and presentation, T cell differentiation, immune regulation and graft-versus-host disease, indicating their association with the development of cSLE. Additionally, we visualised the top five DEGs in monocyte (CD16^-) cells between SLE patients and matched healthy donors and found that IFITM3 (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test), IFI16 (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test) and ISG15 (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test) were significantly upregulated in both aSLE and cSLE patients (Figure 5C-D). Notably, MT2A (Metallothionein 2 A; cSLE: p value = 0, Wilcoxon test) expression was significantly increased in the monocytes (CD16^-) of cSLE patients (Figure 5C), suggesting its potential as a diagnostic biomarker for cSLE.

Figure 5. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in monocyte (CD16^-) cells between SLE patients and matched healthy donors. (A-B) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in monocyte (CD16^-) cells was conducted separately for adults (A) and children (B). (C-D) the top 5 DEGs in monocyte (CD16^-) cells for both adults (C) and children (D) were further analysed using the Seurat R package.

The DEGs related to immune function in CD8⁺ T cells from cSLE patients were similar to those in CD8⁺ T cells from aSLE patients (Figure 6A-B, Table 1). We also found that IFI6 (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test), IFI144L (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test), 1SG15 (aSLE: p value = 0, cSLE: p value = 0, Wilcoxon test), IFI27 (aSLE: p value = 1.32E-293, cSLE: p value = 0, Wilcoxon test) and EPSTI1 (Epithelial Stromal Interaction 1; aSLE: p value = 1.62E-258, cSLE: p value = 0, Wilcoxon test) [21] were significantly upregulated in CD8⁺ T cells of cSLE or aSLE patients (Figure 6C-D). Moreover, similar to monocytes (CD16⁺), HBB (aSLE: p value = 0, Wilcoxon test) and HBA2 (aSLE: p value = 0, Wilcoxon test) were also significantly upregulated in CD8⁺ T cells (Figure 6D), consistent with the findings in monocytes (CD16⁺) cell of aSLE patients. Collectively, these results indicate the promising potential of these genes as diagnostic biomarkers for aSLE.

Figure 6. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in CD8⁺ T cells between SLE patients and matched healthy donors. (A-B) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in CD8⁺ T cells was conducted separately for adults (A) and children (B). (C-D) the top 5 DEGs in CD8⁺ T cells for both adults (C) and children (D) were further analysed using the Seurat R package.

In both aSLE and cSLE patients, we observed a decrease in the percentage of NK cells compared to healthy donors (Figures 1G and 2G). However, no significant differences were observed in the DEGs associated with the NK-cell cytotoxicity pathway between NK cells of cSLE patients and aSLE patients (Figure 7A-B, Table 1). Despite the decreased percentage of NK cells in SLE patients, our findings revealed significant upregulation of IFI6, IFI144L, ISG15, IFI27 and EPSTI1 (aSLE: p value = 5.11E-76, p value = 8.30E-172, p value = 3.29E-69, p value = 3.91E-75, and p value = 1.71E-95; all five genes in cSLE: p value = 0; Wilcoxon test) in SLE patients (Figure 7C-D). Furthermore, we also observed significant upregulation of HBB (aSLE: p value = 2.19E-162, Wilcoxon test), HBA2 (aSLE: p value = 5.79E-105, Wilcoxon test) and IFITM1 (Interferon Induced Transmembrane Protein 1; aSLE: p value = 1.56E-111, Wilcoxon test) in NK cells of aSLE patients (Figure 7D).

Figure 7. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in NK cells between SLE patients and matched healthy donors. (A-B) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in NK cells was conducted separately for adults (A) and children (B). (C-D) the top 5 DEGs in NK cells for both adults (C) and children (D) were further analysed using the Seurat R package.

Further analysis revealed that the lymphocyte activation pathway and B cell activation pathway were significantly upregulated in B cells from aSLE patients compared to cSLE patients (Figure 8A-B, Table 1). Interestingly, in cSLE patients, the top DEGs in B cells were similar to those in NK cells and CD8⁺ T cells (Figure 8C). Notably, we identified RACK1 (Receptor For Activated C Kinase 1; aSLE: p value = 2.62E-231, Wilcoxon test), ATP5F1E (ATP5F1E gene encodes a subunit of mitochondrial ATP synthase; aSLE: p value = 2.08E-229, Wilcoxon test), NOP53 (NOP53 Ribosome Biogenesis Factor; aSLE: p value = 1.04E-223, Wilcoxon test), IGHM (Immunoglobulin Heavy Constant Mu; aSLE: p value = 4.90E-219, Wilcoxon test) and IGHD (Immunoglobulin Heavy Constant Delta; aSLE: p value = 7.43E-209, Wilcoxon test) as highly upregulated in B cells of aSLE patients (Figure 8D).

Figure 8. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in B cells between SLE patients and matched healthy donors. (A-B) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in B cells was conducted separately for adults (A) and children (B). (C-D) the top 5 DEGs in B cells for both adults (C) and children (D) were further analysed using the Seurat R package.

Additionally, regulated network analysis of neutrophil cells from aSLE and cSLE patients revealed that the immune system pathway was upregulated in neutrophil cells from aSLE patients compared to adult healthy donors (Figure 9A, Table 1). Furthermore, we identified significant downregulation of TMEM150B (Transmembrane Protein 150B; aSLE: p value = 1.28E-152, Wilcoxon test), IQSEC2 (IQ Motif And Sec7 Domain ArfGEF 2; aSLE: p value = 1.28E-152, Wilcoxon test), CHN2 (Chimerin 2; aSLE: p value = 4.81E-102, Wilcoxon test), LRP8 (LDL Receptor Related Protein 8; aSLE: p value = 3.57E-77, Wilcoxon test) and USP46 (Ubiquitin Specific Peptidase 46; aSLE: p value = 3.57E-77, Wilcoxon test) in neutrophil cells from aSLE patients (Figure 9C). For instance, TMEM150, a transmembrane protein of DRAM family, was degraded in neutrophil cells from aSLE patients. Therefore, the significantly lower expression of TMEM150B, IQSEC2, CHN2, LRP8 and USP46 in neutrophil cells may affect the prognosis of aSLE patients.

Figure 9. Functional enrichment analysis was performed on differentially expressed genes (DEGs) in neutrophil and pre-B (CD34^-) cells between SLE patients and matched healthy donors. (A) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in neutrophil between aSLE and matched healthy donors. (B) Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in pre-B (CD34^-) cells between cSLE and matched healthy donors. (C) Analysis of the top 5 DEGs enrichment in neutrophil between aSLE and matched healthy donors using the Seurat R package. (D) Analysis of the top 5 DEGs enrichment in pre-B (CD34^-) cells between cSLE and matched healthy donors using the Seurat R package.

Moreover, we observed a significant increase in the differentially expressed genes (DEGs) associated with the immune system process pathway in pre-B cells (CD34^-) from SLE patients compared to pre-B cells (CD34^-) from healthy donors (Figure 9B, Table 1). Specifically, IFI44L, ISG15 and LY6E (cSLE: p value = 1.60E-39, p value = 1.15E-36, and p value = 2.65E-32, Wilcoxon test) were significantly upregulated in pre-B cells (CD34^-) of cSLE patients, and β2M (Beta-2-Microglobulin, cSLE: p value = 1.25E-33, Wilcoxon test) expression was significantly increased in pre-B cells (CD34^-) from cSLE patients (Figure 9D). Notably, TPM2 (Tropomyosin 2, cSLE: p value = 8.46E-29, Wilcoxon test), a member of the actin filament-binding protein family that encodes-tropomyosin [22,23], was significantly decreased in cSLE patients (Figure 9D). Mutations in TPM2 randomly distributed along the chains of tropomyosin isoforms have been associated with congenital myopathies, such as cap disease, congenital myopathies (CM), distal arthrogryposis (DA), and histologically variable disorders of skeletal muscles [22,24]. Taken together, these findings suggest that the identified DEGs in PBMCs from SLE patients have the potential to serve as diagnostic biomarkers for SLE.

Discussion

Understanding the underlying pathogenesis of SLE is crucial due to its complex and unpredictable nature. In this study, we utilised single-cell RNA sequencing (scRNA-seq) analysis to investigate the variations in PBMC proportions and transcriptional profiles between SLE patients and healthy donors. Our results revealed novel alterations in eight PBMC clusters in cSLE patients and seven PBMC clusters in aSLE patients. Notably, we observed an increased proportion of monocytes (CD16⁺), monocytes (CD16^-), and CD8⁺ T cells in both aSLE and cSLE patients, while the percentages of B cells and NK cells were decreased. Additionally, we identified a decreased proportion of pre-B cells (CD34^-) in cSLE patients and a significant upregulation of neutrophil cells in aSLE patients. These findings indicate the presence of shared PBMC clusters between aSLE and cSLE patients, although with distinct alterations in specific cell populations.

In addition, GO analysis revealed the upregulation of various networks related to the development of SLE within the same cell types, including the type I interferon signalling pathway, neutrophil function, lymphocyte function, T cell function, white blood cell function, B cell function, apoptosis and immune regulation. Moreover, we identified several shared DEGs between aSLE and cSLE patients, such as IFITM3, ISG15, IFI16, and LY6E, which exhibited significant upregulation in both patient groups, indicating their potential as therapeutic targets for SLE and aligning with the findings of previous research [17,19,25].

Furthermore, HBA2, HBB, EEF1A1, RACK1, ATP5F1E, NOP53, IGHM and IGHD were identified as potential diagnostic biomarkers for aSLE, while MT2A and TPM2 were identified as potential diagnostic biomarkers for cSLE. Notably, the expression of RACK1, ATP5F1E, NOP53, IGHM and IGHD were specifically upregulated in B cells of aSLE patients, indicating differential transcriptional signatures between aSLE and cSLE. On the other hand, the expression of TPM2, which encodes β-tropomyosin, was significantly downregulated in cSLE patients but not in aSLE patients. Several studies have demonstrated that mutated TPM2 randomly distributes along the chains of tropomyosin isoforms, leading to congenital myopathies such as cap disease, congenital myopathies (CM), and distal arthrogryposis (DA), and histologically variable disorders of skeletal muscle [22–24]. Moreover, it is important to note that SLE can also impact the musculoskeletal system, leading to various symptoms such as arthralgia, arthritis, osteonecrosis (avascular necrosis of bone) and myopathy. Therefore, it is likely that low expression of TPM2 is associated with symptoms of musculoskeletal system involvement in cSLE patients and may as a potential biomaker for cSLE patients.

Neutrophils, as key components of the innate immune system, play a crucial role in defending against infections through various pro-inflammatory functions [26], for instance, by performing chemotactic, phagocytic and bactericidal activities. Their ability to form swarms is commonly observed in acute or suppurative infections caused by diverse pathogens, including abscesses, pneumonia, appendicitis, visceral perforation, scarlet fever, and others. Furthermore, neutrophil swarming can occur in various conditions such as poisonings (acidosis, uremia), tissue damage, malignancy, acute haemorrhage and acute hemolysis [27,28]. In aSLE, the lower expression of TMEM150B, IQSEC2, CHN2, LRP8 and USP46 in neutrophil cells may be associated with a poorer prognosis. Altogether, these genes, identified within neutrophil cells, provide valuable insights into the underlying pathogenesis of aSLE and may serve as potential therapeutic targets for the disease.

In conclusion, our study reveals significant differences in the composition and transcriptional profiles of PBMCs between aSLE and cSLE patients. These findings provide insights into potential diagnostic biomarkers for SLE and offer new opportunities for disease diagnosis and therapeutic interventions.

Methods

scRNA-seq data analysis

The R (version 4.0.2) software and R Seurat package (version 3.2.1) were used for the scRNA-seq data analysis [29]. Samples from adults and children with SLE were separately processed.

Quality control and cell filtering

The matrix of SLE patients and healthy donors were loaded into the Seurat package (version 3.2.1) [29]. To ensure high-quality data for downstream analysis, we applied the following filtering parameters. Firstly, genes that were not detected in at least 100 cells were excluded. Secondly, cells with less than 400 total unique transcripts were removed. Additionally, cells with mitochondrial gene transcripts accounting for more than 5% of the transcripts were considered poor quality and were filtered out. Lastly, cells displaying unusually high unique gene counts (>2,500 genes) were considered outliers and excluded from further analysis.

Normalised and integration

After filtering, the scaled data were first normalised using ‘LogNormalize’ of the Seurat R package. Then, ‘FindVariableFeatures’ of the Seurat R package was used to identify highly variable genes. Following this, anchors were identified using the ‘FindIntegrationAnchors’ function and integrated the patients and healthy donors datasets with ‘IntegrateData’.

Clustering

To identify distinct clusters within the integrated dataset, we utilised the ‘FindClusters’ function in Seurat. This function was applied with the top 15 principal components (PCs) and a resolution parameter of 1.2. Subsequently, marker genes specific to each cluster were identified using the ‘FindAllMarkers’ function in Seurat. The marker genes were tested for differential expression using the Wilcoxon test.

Cell type identification

The R ‘SingleR’ package [30] was used to annotate each cell cluster.

Differential expression analysis

Based on the filtered gene expression matrix by Seurat [29], differential expression analysis between the patients and healthy donors was performed using the ‘FindAllMarkers’ function to identify significantly differentially expressed genes.

Enrichment analysis of differentially expressed genes

To perform enrichment analysis, we utilised the online software g:Profiler (https://biit.cs.ut.ee/gprofiler/) [31]. Then, we applied the Benjamini-Hochberg correction method to adjust the p-values for multiple testing and control the FDR.

Statistical analysis

Statistical analysis was performed using the R software (version 4.0.2). Tests were used to determine data distribution, and depending on the normality of the data, comparisons were performed using the Wilcoxon signed rank test (for two groups, paired) with two-tailed P values unless otherwise stated. Differences were considered to be significant for P values <0.05 (*), < 0.01 (**), < 0.001 (***) and p < 0.0001 (****).

Author’ contributions

J.C. and Y.Z. conceived and designed the project; Y.L., X.Y. and J.W. analysed data; X.J. and Y.L. collected the data; Y.L., X.Y. and J.W. wrote the manuscript; and all authors revised and approved the manuscript.

Disclosure statement

The authors declare no competing interests.

Availability of data and materials

We downloaded the PBMC’ 10X scRNA-seq matrix data of 43 SLE patients (33 cSLE and 10 aSLE patients) and 18 health donors (11 cHD and 7 aHD) from Gene Expression Omnibus (GEO) database under the accession number GSE135779 and GSE142016 [14,17].

Additional information

Funding

This work was supported by the Medical Science and Technology Reserved Talent of Kunming City Health Science and Technology cultivation program of China (grant no. 2022-SW (Reserve) − 00).