Identification of Key Genes in Systemic Lupus Erythematosus through integrated bioinformatics

doi:10.21203/rs.3.rs-3734708/v1

Download PDF

Research Article

Identification of Key Genes in Systemic Lupus Erythematosus through integrated bioinformatics

https://doi.org/10.21203/rs.3.rs-3734708/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background：Systemic lupus erythematosus (SLE) is a chronic autoimmune disease with unclear mechanisms, limiting treatment options. Our study identifies potential core genes of SLE and their clinical applicability.

Method：Using differential expression and weighted gene co-expression network analysis (WGCNA), we identified novel susceptibility modules and associated core genes. Examination of these genes through KEGG and GO analyses revealed their roles. Diagnostic performance of the core genes was evaluated using column line plot models and Receiver Operating Characteristic (ROC) curves. We also assessed the correlation between core genes and immune infiltration and used Mendelian randomization to determine the causal effect of GYPB on SLE.

Results：The gene co-expression network constructed through WGCNA identified 144 key genes associated with SLE. The column line graph model demonstrated strong predictive power for SLE risk, with its diagnostic effectiveness validated by the ROC curve. A causal relationship was established between the top five core genes and immune cell infiltration in SLE. A negative correlation was observed between the gene GYPB and SLE, suggesting that GYPB might serve as a potential therapeutic target.

Conclusion：This investigation provides new insights into SLE molecular mechanisms and potential therapeutic avenues.

SLE

WGCNA

core genes

GYPB

Mendelian randomization

Systemic lupus erythematosus (SLE) is an autoimmune disease typified by chronic systemic inflammation, presenting a high incidence and mortality rate. Predominantly affecting women of reproductive age, the female-to-male ratio is approximately 9:1 [1]. The pathogenesis of SLE is notably complex, with various mechanisms potentially leading to diverse clinical manifestations. The disease's pathogenesis primarily involves a malignant cycle of autoantibodies, immune complexes, inflammatory processes, and tissue damage [2–4]. Current treatment strategies mainly encompass conventional anti-inflammatory medications, immunomodulatory drugs, and corticosteroids, although the treatment efficacy is generally unsatisfactory. Thus, comprehensive research on SLE-related biomarkers and pathogenic mechanisms is crucial for advancing diagnostic and therapeutic strategies.

With the rapid advancement of bioinformatics, it is now possible to utilize extensive microarray data to identify the core genes, interaction networks, and pathways associated with SLE, thereby aiding in improving the disease's diagnosis and treatment. Alongside differential gene expression analysis, weighted gene co-expression network analysis (WGCNA) is a widely utilized bioinformatics technique that effectively explores the correlation between gene expression and clinical attributes [5].

Mendelian Randomization (MR) is another robust method that has gained traction in epidemiological studies [6]. It uses genetic variants associated with exposure factors as instrumental variables to infer causality between these factors and the outcomes, offering a unique approach for examining potential causal links in complex epidemiological scenarios. In the context of SLE, MR can provide invaluable insights into the causal relationship between certain genes and the disease, thus bringing us a step closer to understanding its genetic and immunological aspects.

In this study, we identified differentially expressed genes (DEGs) by analyzing data from both the human SLE group and the control group. We then employed the WGCNA method to identify the gene modules most closely associated with SLE, thereby considerably narrowing the scope of gene selection. Following this, we identified shared genes within the DEGs and the most correlated modules based on WGCNA, and utilized CytoHubba to determine five core genes: SLC4A1, EPB42, FECH, GYPB, and ALAS2. These genes are viewed as potential biomarkers that may aid in the diagnosis of SLE and further enhance understanding of its underlying mechanisms. Ultimately, GYPB was selected for MR analysis to explore its causal relationship with SLE.

2.1 Data source

We downloaded the dataset [GSE61635 (GPL570-55999)] from the GEO database (http://www.ncbi.nlm.nih.gov/geo), which is an open-source platform for retrieving gene expression data [7]. The GSE61635 dataset includes mRNA data from blood samples of 79 patients in the SLE study cohort (99 arrays in total, with some repeated sampling) and 30 healthy volunteers (one array per volunteer, 30 arrays in total). All expression data for the measured genes can be downloaded from the Gene Expression Omnibus (GEO) database.

2.2 Identification of Differentially Expressed Genes

R software (version 4.3.1; https://www.r-project.org/) and Bioconductor software package (http://www.bioconductor.org/) are commonly used tools for normalization and analysis of raw data. We utilize the "limma" package to analyze and filter differentially expressed genes (DEGs) data, where a fold change (|log2FC|) greater than or equal to 1 and an adjusted P-value less than 0.05 are considered statistically significant. After conducting significance analysis on expression levels, we further process the data using the "pheatmap" and "ggplot2" packages to generate heatmaps and volcano plots illustrating the DEGs.

2.3 Weighted Gene Co-expression Network Analysis (WGCNA)

WGCNA is a bioinformatics technique used to construct gene co-expression networks and identify network modules closely linked with clinical traits [5]. Utilizing the "WGCNA" R package, we constructed a gene co-expression network for SLE. Initially, we generated a sample clustering tree to evaluate potential outliers. Subsequently, the adjacency matrix (AM) was transformed into a topological overlap measure (TOM) matrix. We then identified hierarchical clustering genes using the dynamic tree cutting algorithm and grouped similar genes into the same module based on a TOM-based dissimilarity measure. Following this, we calculated the module membership (MM) and gene significance (GS), and analyzed the correlation between sample modules and clinical traits using the Pearson correlation test. Conclusively, the module genes most highly correlated with clinical traits were selected for further analysis.

2.4 Screening of candidate key genes and GO/KEGG analysis

We utilized the "VennDiagram" package in R software to identify the intersection of the most relevant genes and differentially expressed genes in the WGCNA module. These intersecting genes are considered candidate key genes related to the pathogenesis of SLE. Following this, we conducted gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the R package "ClusterProfiler". This allowed us to evaluate the functions and pathways of these intersecting genes, thereby facilitating a deeper understanding of the potential mechanisms underlying the progression of SLE.

2.5 Protein-Protein Interaction (PPI) Network Analysis and Identification of Core Genes

In PPI network analysis, we utilized STRING (http://string-db.org) to retrieve interacting genes and set a threshold of 0.4. Additionally, we employed Cytoscape software (version 3.7.1;http://www.cytoscape.org/) for the visualization of the PPI network. Lastly, to identify core genes, we utilized the CytoHubba tool and ranked important genes in the PPI network using the Degree algorithm.

2.6 Construction of the nomogram model

We employed the "rms" package in R software to construct a nomogram model for predicting the risk of systemic lupus erythematosus [8]. Subsequently, we used the "ROC" package to generate receiver operating characteristic (ROC) curves to validate the diagnostic efficacy of the candidate core genes. The accuracy was evaluated using the area under the ROC curve (AUC), with values within the range of 0.9 ≤ AUC < 1 considered indicative of high accuracy.

2.7 Immunocyte Analysis

In order to study the functionality of immune cells in systemic lupus erythematosus (SLE), we evaluated the immune cell infiltration levels of 130 immune cells in the SLE group using the CIBERSORT analysis method [9]. The "corrplot" and "pheatmap" packages in R software were used to perform correlation analysis and visualization of these 20 types of infiltrating immune cells. Moreover, the "vioplot" package in R software was used to generate violin plots, visually demonstrating the differences in immune cell infiltration between the two groups of samples.

2.8. Correlation analysis for identifying candidate core genes and infiltrating immune cells

We conducted a correlation analysis between the candidate core genes and the levels of infiltrating immune cells using the "reshape2", "ggpubr", and "ggExtra" packages in R software. This analysis produced bubble plots and scatter plots for visual representation of the findings.

2.9. Mendel randomization analysis

To conduct MR analysis, we utilized the GWAS ID of core genes and systemic lupus erythematosus obtained from publicly available genome-wide association studies (GWAS). We employed various methods such as IVW, MR Egger, Weighted median, and Weighted mode to determine the causal relationship between candidate core genes and systemic lupus erythematosus [10]. The primary method used for MR analysis was IVW, with inverse variance weighting and MR-Egger regression employed for heterogeneity detection. Horizontal pleiotropy was evaluated by examining the intercept term in MR-Egger regression. The analysis considered the intercept term to be close to zero (< 0.1) with a P-value greater than 0.05 as an indication of no evidence of horizontal pleiotropy. This interpretation supports the reliability of the Mendelian randomization analysis results [11]. Additionally, we conducted sensitivity analysis using leave-one-out method to investigate the impact of individual SNPs on the results. All of the aforementioned analysis methods were performed using the TwoSampleMR package.

3.1 DEG screening

We initially acquired the SLE dataset (GSE61635) from the GEO database and identified the DEGs within this dataset. In comparison to the control group, we identified 687 DEGs (comprising 477 upregulated and 210 downregulated) in the SLE group (Fig. 1A, 1B, and Supplementary Table S1).

Figure 1 The differential gene expression between the SLE group and the control group. (A) The heatmap represents the differential expression analysis from GSE61635. High gene expression is denoted by the color red, while low expression is indicated by the color blue. (B) The volcano plot depicts the differential expression analysis from GSE61635. Upregulated genes are represented by the color red, downregulated genes by the color green, and non-differentially expressed genes are shown in gray.

3.2 Construction of the WGCNA network and identification of modules associated with SLE

In order to determine the potential association between gene modules and SLE, we conducted a WGCNA analysis on all candidate genes from GSE61635 (Fig. 2A, 2B). We use a power of β = 10 (scale-free R2 = 0.8) as a soft threshold to ensure a scale-free network. During this process, we identified 14 distinct modules (Fig. 2C). After analysis, we obtained positive correlation coefficients. Finally, we screened the magenta module in the GSE61635 dataset ( Fig. 2D, Supplementary Table S2) and selected key genes within the magenta module for further analysis based on the criteria of MM > 0.8 and GS > 0.5 (Supplementary Table S3).

Figure 2 WGCNA was used to identify gene modules associated with SLE in the GEO dataset. (A) The clustering of genes based on their expression profiles. Hierarchical clustering grouped genes with similar expression patterns into distinct co-expression modules, represented by different colors in the dendrogram. (B) Hierarchical clustering of all genes in the GSE61635 dataset based on the topological overlap matrix (1-TOM), where each branch represents a gene and co-expression modules are indicated by different colors. (C) Module-trait heatmap showing the correlation between clustered gene modules in the GSE61635 dataset and SLE, with corresponding correlation coefficients and p-values for each module. (D) Scatter plot of the magenta module (highlighted in green color) in the GSE61635 dataset, showing the strongest positive correlation with SLE.

3.3 GO/KEGG analysis

We identified 144 overlapping genes between the most relevant module genes and differentially expressed genes from the WGCNA analysis. These overlapping genes may serve as candidate core genes playing important roles in the pathogenesis of SLE (Fig. 3A). To explore the potential functions of these 144 genes, we performed GO and KEGG enrichment analyses (Fig. 3B, 3C). The GO analysis revealed enrichment for biological processes related to hematopoietic stem cell proliferation, hematopoietic precursor cell differentiation, Golgi-to-lysosome transport, transcriptional regulation, apoptosis, microtubule dynamics, DNA repair, and recombination. Enriched cellular components included spindle and microtubule structures, DNA polymerase complexes, and transferase complexes. Molecular functions included various transferase activities, DNA binding and repair activities, and microtubule binding. The KEGG analysis showed enrichment for pathways related to porphyrin metabolism, lifespan regulation, endoplasmic reticulum protein processing, AMPK signaling, and amino acid metabolism.

Figure 3 Filtering candidate core genes and conducting verification. (A) A Venn diagram shows 144 overlapping candidate key genes. (B) GO enrichment analysis of candidate hub genes. BP represents biological processes. CC represents cellular components. MF represents molecular functions. (C) KEGG pathway analysis of candidate hub genes. The redder the color, the more obvious the enrichment.

3.4 PPI network analysis of candidate core genes

To identify key hub genes among the 144 overlapping genes from the WGCNA and differential expression analyses, we constructed a protein-protein interaction (PPI) network using the STRING online database (Fig. 4A). The top 10 ranked genes by connectivity degree were visualized using Cytoscape software (Fig. 4B). These hub genes were SLC4A1, EPB42, FECH, GYPB, ALAS2, AHSP, GATA1, KLF1, SNCA, and DMTN. We selected the top 5 ranked hub genes (SLC4A1, EPB42, FECH, GYPB, and ALAS2) as potential key drivers of SLE pathogenesis for further functional characterization.

Figure 4 Construction of PPI network. (A) PPI network graph of overlapping core genes. (B) Top 10 core genes of the interactome network identified by degree algorithm. The deeper the color red, the higher the score.

3.5 Construction of a column-line chart model for risk prediction of systemic lupus erythematosus (SLE).

To evaluate the diagnostic potential of the 5 hub genes, we developed a nomogram model to predict risk of SLE (Fig. 5A-B). The nomogram showed good predictive performance for SLE. We also assessed the diagnostic accuracy of the individual genes by plotting ROC curves and calculating the AUC values (Fig. 5C). The AUCs were 0.968 for SLC4A1, 0.974 for EPB42, 0.930 for FECH, 0.931 for GYPB, and 0.907 for ALAS2, indicating high diagnostic accuracy. This nomogram model utilizing the 5 hub genes successfully distinguished SLE cases from controls.

Figure 5 Using the nomogram model to predict the risk of systemic lupus erythematosus. (A) Column chart of the core genes in the nomogram model. (B) Calibration curve of the nomogram model. (C) Evaluating the diagnostic effectiveness of the nomogram model and key hub genes using ROC curves.

3.6 Evaluation of immune cell infiltration in systemic lupus erythematosus

To assess the immune microenvironment in SLE, we utilized CIBERSORT to estimate the relative abundance of 20 immune cell types in each sample (Fig. 6A). Correlation analysis revealed complex relationships between the immune cell types (Fig. 6B). Positive correlations existed between multiple cell types including memory B cells, plasma cells, Tfh cells, neutrophils, M1/M2 macrophages, eosinophils, resting/activated memory CD4 + T cells, naive B cells, resting NK cells, and activated dendritic cells. Negative correlations were found between neutrophils and monocytes, dendritic cells, B cells, NK cells, gamma delta T cells, CD8 + T cells, and activated CD4 + T cells. Additional negative correlations occurred between M1 macrophages and monocytes, between naive and CD8 + T cells, and between resting and CD8 + memory T cells. Comparing SLE and control samples, significant differences emerged in naive B cells, plasma cells, resting/activated memory CD4 + T cells, gamma delta T cells, M0 macrophages, dendritic cells, and neutrophils (Fig. 6C).

Figure 6 Systemic lupus erythematosus immune relation. (A) The relative distribution of 20 immune cells in systemic lupus erythematosus samples. (B) Correlation relationship of 20 immune cells. (C) Differences of immune cells between the control group and the experimental group.

3.7 The relationship between candidate core genes and infiltrating immune cells

We examined correlations between the 5 candidate hub genes and immune cell abundances to uncover potential relationships. SLC4A1 showed negative correlation with activated dendritic cells (R=-0.24, P = 0.015) and positive correlation with eosinophils (R = 0.32, P = 0.0013) (Fig. 7). Similarly, EPB42 was negatively correlated with activated dendritic cells (R=-0.25, P = 0.013) and positively correlated with eosinophils (R = 0.28, P = 0.0045) (Fig. 8). FECH positively correlated with eosinophils (R = 0.37, P = 0.00019), negatively correlated with neutrophils (R=-0.41, P = 3.7e-05), and positively correlated with activated memory CD4 + T cells (R = 0.25, P = 0.012), CD8 + T cells (R = 0.38, P = 8.5e-05), and γδT cells (R = 0.21, P = 0.033) (Fig. 9). GYPB showed positive correlation with plasma cells (R = 0.24, P = 0.015) (Fig. 10). Finally, ALAS2 negatively correlated with activated dendritic cells (R=-0.24, P = 0.015) and positively correlated with eosinophils (R = 0.32, P = 0.0013) (Fig. 11).

Figure7 The correlation between candidate core genes and infiltrating immune cells. (A) The correlation between SLC4A1 and infiltrating immune cells. (B) The correlation between SLC4A1 and Dendritic cells activated. (C) The correlation between SLC4A1 and Eosinophils.

Figure 8 The correlation between candidate core genes and infiltrating immune cells. (A) The correlation between EPB42 and infiltrating immune cells. (B) The correlation between EPB42 and Dendritic cells activated. (C) The correlation between EPB42 and Eosinophils.

Figure 9 The correlation between candidate core genes and infiltrating immune cells.(A) The correlation between FECH and infiltrating immune cells.(B) The correlation between FECH and Eosinophils.(C) The correlation between FECH and Neutrophils. (D) The correlation between FECH and T cells CD4 memory activated.(E) The correlation between FECH and T cells CD8. (F) The correlation between FECH and T cells gamma delta.

Figure 10 The correlation between candidate core genes and infiltrating immune cells. (A) The correlation between GYPB and infiltrating immune cells. (B) The correlation between GYPB and plasma cells

Figure 11 The correlation between candidate core genes and infiltrating immune cells. (A) The correlation between ALAS2 and infiltrating immune cells. (B) The correlation between ALAS2 and eosinophils.

3.8 MR analysis.

We performed further analysis on the hub gene GYPB due to the lack of GWAS data for the other genes. Table S4 shows the SNP features of GYPB associated with SLE, with no evidence of weak instrumental variables. The causal relationship between GYPB variants and SLE risk was assessed through MR analysis (Fig. 12A-B). The IVW method demonstrated a negative correlation between GYPB and SLE risk (OR 0.620, 95% CI 0.406–0.948, p = 0.027). The weighted median method confirmed this protective effect (OR 0.583, 95% CI 0.343–0.995, p = 0.048). While there was no evidence of pleiotropy based on the Egger intercept test (p = 0.706). The leave-one-out analysis retained consistent results after removing each SNP individually (Fig. 12C), indicating the causal effects are not driven by any single variant. In summary, multiple MR methods converged to suggest GYPB is causally associated with reduced SLE risk.

Figure 12 Mendelian randomization study results. (A) Scatter plot shows the causal relationship between GYPB and the risk of systemic lupus erythematosus. (B) Forest plot shows the causal relationship between each SNP and the risk of systemic lupus erythematosus. (C) Visualization graph of the causal effect of GYPB on the risk of systemic lupus erythematosus when one SNP is omitted.

SLE is a chronic autoimmune systemic inflammatory disease characterized by high heterogeneity, which may involve multiple organ systems [12]. Currently, the pathogenesis of SLE is not fully understood, but it is known to involve immune dysfunction, including abnormalities in the number and function of immune cells [13]. The main treatment options for SLE include non-specific therapies such as glucocorticoids or immunosuppressants. However, conventional treatments have limited effectiveness and poor prognosis. Therefore, exploring the molecular features and mechanisms underlying the development of SLE is crucial for effective prevention, diagnosis, and treatment of the disease.

The purpose of diagnostic biomarkers is to detect patients with pathological changes. Currently, researchers are working on validating biomarkers associated with systemic lupus erythematosus (SLE). In a recent study, MX2 was found to be an immune-related biomarker for SLE, which can be used for predicting the diagnosis and disease activity of SLE. MX2 activates the NOD-like receptor signaling pathway, promotes neutrophil infiltration, and may worsen the condition of SLE [14]. Another study identified 28 hub genes, such as RPS7, RPL19, RPS17, and RPS19, playing important roles in the occurrence, development, and prognosis of SLE [15]. Lv Juan et al. found common expression interactions among SLE-associated gene CD40LG extracted from CTD and seven differentially expressed long non-coding RNAs (HCG27, LINC02555, LINC02210, DHRS4-AS1, MIR600HG, DANCR, and LINC01278). The research observed downregulation of CD40LG and HLA-DOA expression, upregulation of FCGR3A and HIST1H2BE expression in the validation dataset GES46907. Additionally, downregulation of CD40LG was also confirmed by qRT-PCR [16]. Xingwang Zhao et al. found through comprehensive bioinformatics analysis that IFI27 is positively correlated with immune function, increased monocyte infiltration, and decreased resting NK cell infiltration, which may be related to the pathogenesis of SLE [17].Researcher Zhihang Jiang et al. utilized bioinformatics methods and integrated multiple data sets to discover 10 potential diagnostic biomarkers for systemic lupus erythematosus (SLE). These biomarkers are IFI44, IFI44L, EIF2AK2, IFIT3, IFITM3, ZBP1, TRIM22, PRIC285. Among them, IFI44 was identified as the most promising biomarker using five machine learning algorithms. To validate the diagnostic performance of IFI44 in specific cohorts, the researchers conducted qRT-PCR and ROC curve analysis. In addition, results from immune cell infiltration showed a significantly higher proportion of central memory CD8 + T cells in SLE patients, which correlated positively with the selected biomarkers [18].In order to identify new biomarkers as diagnostic markers or therapeutic targets for systemic lupus erythematosus (SLE) and achieve personalized and differentiated precision medicine, this study employed bioinformatics analysis methods. We downloaded gene expression profiling data of whole blood samples from SLE patients and healthy controls from the GEO database, and used the limma package to identify differentially expressed genes between SLE and normal population [15].In order to identify key genes and pathways in the pathogenesis of SLE, we performed WGCNA analysis and identified 14 co-expression modules. Among them, we selected the magenta module that was most relevant to SLE and identified it as the most important module associated with SLE. Next, we compared the genes extracted from this module with differentially expressed genes and identified 144 overlapping genes. These overlapping genes were considered candidate key genes related to the pathogenesis of SLE.We further conducted functional analysis and protein-protein interaction network analysis for these overlapping genes. Finally, we screened out 10 candidate core genes, which are SLC4A1, EPB42, FECH, GYPB, ALAS2, AHSP, GATA1, KLF1, SNCA, and DMTN. These genes play important roles in systemic lupus erythematosus and have the potential to become therapeutic targets or diagnostic markers. Our research results demonstrate excellent performance of the predictive model for SLE. We ultimately confirmed SLC4A1, EPB42, FECH, GYPB, and ALAS2 as core genes and calculated their ROC curves to evaluate their diagnostic accuracy. ROC analysis demonstrated high accuracy and sensitivity of each biomarker in the diagnosis of systemic lupus erythematosus, and the combination of these five biomarkers also showed promising diagnostic efficacy.

SLC4A1 is one of the commonly differentially expressed genes and belongs to the family of anion exchangers, primarily expressed in the red blood cell membrane [19]. SLC4A1 can stimulate the functional activity of innate immune cells and participate in immune response, T cell activation, regulation of toll-like receptor binding pathways, and granulocyte activation [20].According to reports, overexpression of SLC4A1 can cause certain diseases, such as hereditary spherocytosis caused by the instability of red blood cell membranes and hereditary distal renal tubular acidosis [21, 22]. Kawamoto et al. found that colorectal cancer cells stimulate the production of antibodies binding to red blood cell membranes through the overexpression of band 3, leading to immune-related anemia [23, 24].EPB42, also known as erythrocyte membrane protein band 4.2, is a protein that binds to ATP and regulates the interaction between protein 3 and strong protein. According to reports, EPB42 plays a role in modulating the shape and mechanical properties of red blood cells, and is closely associated with hereditary spherocytosis and autosomal recessive hemolytic anemia [25, 26]. Ferrochelatase (FECH) is the final enzyme in heme biosynthesis and plays an important role in choroidal neovascularization, retinal neovascularization, and congenital erythropoietic porphyria [27–29]. GYPB is the main intrinsic membrane protein of red blood cells [30], which specifically binds to native myeloperoxidase (MPO) and causes various changes in the biophysical properties of cells. These changes include the phenomenon of "lipid freezing", alteration of transmembrane potential, dynamic changes in red blood cell morphology and size, increased sensitivity to acidosis and osmotic hemolysis, and decreased cell deformability [31]. The 5-aminolevulinic acid synthase (ALAS) is encoded by the ALAS2 gene, and it catalyzes the synthesis of 5-aminolevulinic acid by reacting succinyl-CoA with glycine, which is the first step in hemoglobin synthesis [32]. Increasing hemoglobin synthesis can promote hematopoietic cell differentiation [33]. Upregulation of ALAS2 expression enhances heme production, hemoglobinization, and red blood cell formation [34]. Recent studies have shown that the ALAS2 gene is primarily associated with disorders such as X-linked sideroblastic anemia and primary myelofibrosis. The above-mentioned genes (SLC4A1, EPB42, FECH, GYPB, ALAS2) may be involved in the occurrence and development of systemic lupus erythematosus, although research in this field is still lacking.

Research has shown that immune cell infiltration plays an important role in the development of SLE [35]. Therefore, it is of profound significance to search for specific diagnostic markers and analyze the pattern of immune cell infiltration in SLE for improving the prognosis of SLE patients. To further investigate the role of immune cell infiltration in SLE, we used CIBERSORT to evaluate the expression levels of 20 immune cells. The results showed a close correlation between immune cell infiltration and the development of SLE. At the same time, we also observed a connection between the infiltration of certain immune cells and the expression of biomarker genes. The study by Huang Xinfang et al. mentioned that plasmacytoid dendritic cells (pDCs) are immune cells that primarily produce type I interferon (IFN-I) and are widely associated with SLE pathogenesis and activity [36]. IFN-I production by pDCs leads to immune activation and inflammation in SLE. The study by Chen Pingmin et al. found that CD8 + T cells play a role in peripheral tolerance, and reversing their depleted phenotype may alleviate systemic autoimmune and infection risks in SLE patients [37].

In this study, we found that five candidate biomarkers were significantly correlated with infiltration of multiple immune cells, including Eosinophils, Dendritic cells activated, T cell CD8, Memory activated T cells, Gamma delta T cells, Neutrophils, and Plasma cells. These results suggest that immune cell infiltration plays an important role in SLE pathogenesis and will be the focus of further research into diagnostic and prognostic biomarkers.

Glycophorin B (GYPB) is a glycoprotein that is highly expressed on the surface of red blood cells and serves as a receptor for Plasmodium, the pathogen causing malaria, making it a key factor in Plasmodium invasion [38]. Red blood cells (RBCs) must maintain optimal deformability to efficiently deliver oxygen throughout the body. RBC deformability depends on the fluidity and structural integrity of the cell membrane. Myeloperoxidase (MPO) is an enzyme secreted during inflammation that catalyzes reactions producing reactive oxidants and diffusible radicals. MPO was found to specifically bind major RBC membrane proteins, including GYPB. This MPO binding altered biophysical properties of RBCs, including decreased membrane lipid fluidity, shifted transmembrane potential, and dynamic changes in cell morphology and size. Furthermore, MPO binding heightened susceptibility to acidic and osmotic hemolysis and substantially reduced RBC deformability. Impaired RBC deformability is associated with several cardiovascular diseases including atherosclerosis, ischemic heart disease, sepsis, and diabetes. The marked reduction in RBC deformability induced by MPO binding provides a potential mechanistic link between chronic inflammation and increased risk of cardiovascular diseases. MPO-mediated rigidity in RBCs could hinder their passage through microvasculature blood vessels, resulting in tissue ischemia. [39–41]. Among the five core genes identified, GYPB emerged as potentially protective against SLE risk. While the mechanisms are still unclear, our findings suggest GYPB levels may be relevant to SLE pathogenesis. Further experimental studies are needed to validate the association between GYPB and SLE and uncover the potential molecular interactions involved.

Our research holds several advantages. Primarily, our MR study design provides evidence for a causal relationship between five core biomarkers and SLE, thereby circumventing common bias issues inherent in observational epidemiological studies. Secondly, our dataset comprises core biomarker and outcome data from European populations, which aids in minimizing the potential bias caused by differences in genetic backgrounds. Lastly, to ensure the stability of our results, we carried out an MR-Egger regression test, which found no evidence of directional pleiotropy.

However, our research is not without limitations. Although we leveraged bioinformatics to analyze the relationship between core genes and potential functions associated with SLE, further biological experiments are required to validate the specific mechanisms of the selected core genes. Additionally, while MR is a powerful tool for inferring causal relationships, the results still necessitate further validation through experimental research.

In this study, we leveraged weighted gene co-expression network analysis to identify candidate genes and pathways involved in the pathogenesis of SLE. WGCNA revealed a module of co-expressed genes that was significantly associated with SLE. Overlapping this module with differentially expressed genes highlighted 144 candidate hub genes for further characterization. Protein-protein interaction network analysis identified 5 top-ranked hub genes (SLC4A1, EPB42, FECH, GYPB, and ALAS2) that may play key roles in SLE development. These genes showed potential as diagnostic biomarkers based on ROC curve analysis. Additional correlation analyses uncovered relationships between the hub genes and infiltrating immune cell populations, including dendritic cells, eosinophils, and CD4 + T cell subsets. Finally, MR analysis provided evidence for a causal protective effect of GYPB variants against SLE risk. In summary, integration of WGCNA with multi-omics datasets enabled us to elucidate candidate mechanisms and biomarkers for SLE pathogenesis. The identified genes and pathways warrant further functional validation but may ultimately lead to improved diagnostic and therapeutic strategies for this complex autoimmune disorder.

Ethical statement：This article was obtained from a public database and has no ethical concerns.

Conflict of Interest ：The authors declare that they have no conflict of interest.

Funding Declaration：There was no Funding.

Author Contribution declaration: YHY designed the study. YHY and AMZ performed the data analysis and wrote the manuscript. HYJ supervised this study and was responsible for the visualization of the results. All authors read and approved the final manuscript.

Li H-D, You Y-K, Shao B-Y, Wu W-F, Wang Y-F, Guo J-B, et al. Roles and crosstalks of macrophages in diabetic nephropathy. Front Immunol. 2022;13:1015142.
Tsokos GC, Lo MS, Costa Reis P, Sullivan KE. New insights into the immunopathogenesis of systemic lupus erythematosus. Nat Rev Rheumatol. 2016;12:716–30.
Durcan L, O’Dwyer T, Petri M. Management strategies and future directions for systemic lupus erythematosus in adults. Lancet. 2019;393:2332–43.
Kiriakidou M, Ching CL. Systemic Lupus Erythematosus. Ann Intern Med. 2020;172:ITC81–96.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Yuan S, Mason AM, Carter P, Vithayathil M, Kar S, Burgess S, et al. Selenium and cancer risk: Wide-angled Mendelian randomization analysis. Int J Cancer. 2022;150:1134–40.
Wang Z, Monteiro CD, Jagodnik KM, Fernandez NF, Gundersen GW, Rouillard AD, et al. Extraction and analysis of signatures from the Gene Expression Omnibus by the crowd. Nat Commun. 2016;7:12846.
Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol. 2008;26:1364–70.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.
Bowden J, Davey Smith G, Haycock PC, Burgess S. Consistent Estimation in Mendelian Randomization with Some Invalid Instruments Using a Weighted Median Estimator. Genet Epidemiol. 2016;40:304–14.
Verbanck M, Chen C-Y, Neale B, Do R. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet. 2018;50:693–8.
Yan R, Jiang H, Gu S, Feng N, Zhang N, Lv L, et al. Fecal Metabolites Were Altered, Identified as Biomarkers and Correlated With Disease Activity in Patients With Systemic Lupus Erythematosus in a GC-MS-Based Metabolomics Study. Front Immunol. 2020;11:2138.
Zhang X, Ouyang X, Xu Z, Chen J, Huang Q, Liu Y, et al. CD8 + CD103 + iTregs Inhibit Chronic Graft-versus-Host Disease with Lupus Nephritis by the Increased Expression of CD39. Mol Ther. 2019;27:1963–73.
Meng X-W, Cheng Z-L, Lu Z-Y, Tan Y-N, Jia X-Y, Zhang M. MX2: Identification and systematic mechanistic analysis of a novel immune-related biomarker for systemic lupus erythematosus. Front Immunol. 2022;13:978851.
Hu H, He C. Identification of Diagnostic Gene Markers and Immune Infiltration in Systemic Lupus. Comput Math Methods Med. 2022;2022:3386999.
Lv J, Chen L, Wang X, Gao Q, Zhao L. Immune-relevant genes of systemic lupus erythematosus by transcriptome profiling analysis. Cytokine. 2022;158:155975.
Zhao X, Zhang L, Wang J, Zhang M, Song Z, Ni B, et al. Identification of key biomarkers and immune infiltration in systemic lupus erythematosus by integrated bioinformatics analysis. J Transl Med. 2021;19:35.
Jiang Z, Shao M, Dai X, Pan Z, Liu D. Identification of Diagnostic Biomarkers in Systemic Lupus Erythematosus Based on Bioinformatics Analysis and Machine Learning. Front Genet. 2022;13:865559.
Arakawa T, Kobayashi-Yurugi T, Alguel Y, Iwanari H, Hatae H, Iwata M, et al. Crystal structure of the anion exchanger domain of human erythrocyte band 3. Science. 2015;350:680–4.
Zhang D, Zhang N, Wang Y, Zhang Q, Wang J, Yao J. Analysis of differentially expressed genes in individuals with noninfectious uveitis based on data in the gene expression omnibus database. Med (Baltim). 2022;101:e31082.
Lux SE, John KM, Kopito RR, Lodish HF. Cloning and characterization of band 3, the human erythrocyte anion-exchange protein (AE1). Proc Natl Acad Sci U S A. 1989;86:9089–93.
Mohebbi N, Wagner CA. Pathophysiology, diagnosis and treatment of inherited distal renal tubular acidosis. J Nephrol. 2018;31:511–22.
Kawamoto S, Kamesaki T, Masutani R, Kitao A, Hatanaka K, Imakita M, et al. Ectopic expression of band 3 anion transport protein in colorectal cancer revealed in an autoimmune hemolytic anemia patient. Hum Pathol. 2019;83:193–8.
Kitao A, Kawamoto S, Kurata K, Hayakawa I, Yamasaki T, Matsuoka H, et al. Band 3 ectopic expression in colorectal cancer induces an increase in erythrocyte membrane-bound IgG and may cause immune-related anemia. Int J Hematol. 2020;111:657–66.
Chonat S, Risinger M, Sakthivel H, Niss O, Rothman JA, Hsieh L, et al. The Spectrum of SPTA1-Associated Hereditary Spherocytosis. Front Physiol. 2019;10:815.
Qin L, Nie Y, Zhang H, Chen L, Zhang D, Lin Y, et al. Identification of new mutations in patients with hereditary spherocytosis by next-generation sequencing. J Hum Genet. 2020;65:427–34.
Long Z-B, Wang Y-W, Yang C, Liu G, Du Y-L, Nie G-J, et al. Identification of FECH gene multiple variations in two Chinese patients with erythropoietic protoporphyria and a review. J Zhejiang Univ Sci B. 2016;17:813–20.
Basavarajappa HD, Sulaiman RS, Qi X, Shetty T, Sheik Pran Babu S, Sishtla KL, et al. Ferrochelatase is a therapeutic target for ocular neovascularization. EMBO Mol Med. 2017;9:786–801.
Pran Babu SPS, White D, Corson TW. Ferrochelatase regulates retinal neovascularization. FASEB J. 2020;34:12419–35.
Yang S-Y, Zeng L-Y, Li C, Yan H. Correlation between an ABO Blood Group and Primary Femoral Head Necrosis: A Case-Control Study. Orthop Surg. 2020;12:450–6.
Algady W, Louzada S, Carpenter D, Brajer P, Färnert A, Rooth I, et al. The Malaria-Protective Human Glycophorin Structural Variant DUP4 Shows Somatic Mosaicism and Association with Hemoglobin Levels. Am J Hum Genet. 2018;103:769–76.
Kafina MD, Paw BH. Intracellular iron and heme trafficking and metabolism in developing erythroblasts. Metallomics. 2017;9:1193–203.
Rio S, Gastou M, Karboul N, Derman R, Suriyun T, Manceau H, et al. Regulation of globin-heme balance in Diamond-Blackfan anemia by HSP70/GATA1. Blood. 2019;133:1358–70.
La P, Oved JH, Ghiaccio V, Rivella S. Mitochondria Biogenesis Modulates Iron-Sulfur Cluster Synthesis to Increase Cellular Iron Uptake. DNA Cell Biol. 2020;39:756–65.
Moulton VR, Tsokos GC. T cell signaling abnormalities contribute to aberrant immune cell function and autoimmunity. J Clin Invest. 2015;125:2220–7.
Huang X, Dorta-Estremera S, Yao Y, Shen N, Cao W. Predominant Role of Plasmacytoid Dendritic Cells in Stimulating Systemic Autoimmunity. Front Immunol. 2015;6:526.
Chen P-M, Tsokos GC. T Cell Abnormalities in the Pathogenesis of Systemic Lupus Erythematosus: an Update. Curr Rheumatol Rep. 2021;23:12.
Dankwa S, Chaand M, Kanjee U, Jiang RHY, Nobre LV, Goldberg JM, et al. Genetic Evidence for Erythrocyte Receptor Glycophorin B Expression Levels Defining a Dominant Plasmodium falciparum Invasion Pathway into Human Erythrocytes. Infect Immun. 2017;85:e00074–17.
Gorudko IV, Sokolov AV, Shamova EV, Grigorieva DV, Mironova EV, Kudryavtsev IV, et al. Binding of human myeloperoxidase to red blood cells: Molecular targets and biophysical consequences at the plasma membrane level. Arch Biochem Biophys. 2016;591:87–97.
Li Z, Hu S, Cheng K. Platelets and their biomimetics for regenerative medicine and cancer therapies. J Mater Chem B. 2018;6:7354–65.
Su T, Huang K, Ma H, Liang H, Dinh P-U, Chen J, et al. Platelet-Inspired Nanocells for Targeted Heart Repair After Ischemia/Reperfusion Injury. Adv Funct Mater. 2019;29:1803567.

No competing interests reported.

Supplementarytables.xlsx

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of Key Genes in Systemic Lupus Erythematosus through integrated bioinformatics

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials and Methods

2.1 Data source

2.2 Identification of Differentially Expressed Genes

2.3 Weighted Gene Co-expression Network Analysis (WGCNA)

2.4 Screening of candidate key genes and GO/KEGG analysis

2.5 Protein-Protein Interaction (PPI) Network Analysis and Identification of Core Genes

2.6 Construction of the nomogram model

2.7 Immunocyte Analysis

2.8. Correlation analysis for identifying candidate core genes and infiltrating immune cells

2.9. Mendel randomization analysis

3. Results

3.1 DEG screening

3.2 Construction of the WGCNA network and identification of modules associated with SLE

3.3 GO/KEGG analysis

3.4 PPI network analysis of candidate core genes

3.5 Construction of a column-line chart model for risk prediction of systemic lupus erythematosus (SLE).

3.6 Evaluation of immune cell infiltration in systemic lupus erythematosus

3.7 The relationship between candidate core genes and infiltrating immune cells

3.8 MR analysis.

4. Discussion

5. Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1