Dysregulation of Signaling Pathways Due to Differentially Expressed Genes From the B-Cell Transcriptomes of Systemic Lupus Erythematosus Patients - A Bioinformatics Approach.

Systemic lupus erythematosus (SLE) is an autoimmune inflammatory disorder that is clinically complex and has increased production of autoantibodies. Via emerging technologies, researchers have identified genetic variants, expression profiling of genes, animal models, and epigenetic findings that have paved the way for a better understanding of the molecular and genetic mechanisms of SLE. Our current study aimed to illustrate the essential genes and molecular pathways that are potentially involved in the pathogenesis of SLE. This study incorporates the gene expression profiling data of the microarray dataset GSE30153 from the Gene Expression Omnibus (GEO) database, and differentially expressed genes (DEGs) between the B-cell transcriptomes of SLE patients and healthy controls were screened using the GEO2R web tool. The identified DEGs were subjected to STRING analysis and Cytoscape to explore the protein-protein interaction (PPI) networks between them. The MCODE (Molecular Complex Detection) plugin of Cytoscape was used to screen the cluster subnetworks that are highly interlinked between the DEGs. Subsequently, the clustered DEGs were subjected to functional annotation with ClueGO/CluePedia to identify the significant pathways that were enriched. For integrative analysis, we used GeneGo MetacoreTM, a Cortellis Solution software, to exhibit the Gene Ontology (GO) and enriched pathways between the datasets. Our study identified 4 upregulated and 13 downregulated genes. Analysis of GO and functional enrichment using ClueGO revealed the pathways that were statistically significant, including pathways involving T-cell costimulation, lymphocyte costimulation, negative regulation of vascular permeability, and B-cell receptor signaling. The DEGs were mainly enriched in metabolic networks such as the phosphatidylinositol-3,4,5-triphosphate pathway and the carnitine pathway. Additionally, potentially enriched pathways, such as the signaling pathways induced by oxidative stress and reactive oxygen species (ROS), chemotaxis and lysophosphatidic acid signaling induced via G protein-coupled receptors (GPCRs), and the androgen receptor activation pathway, were identified from the DEGs that were mainly associated with the immune system. Four genes (EGR1, CD38, CAV1, and AKT1) were identified to be strongly associated with SLE. Our integrative analysis using a multitude of bioinformatics tools might promote an understanding of the dysregulated pathways that are associated with SLE development and progression. The four DEGs in SLE patients might shed light on the pathogenesis of SLE and might serve as potential biomarkers in early diagnosis and as therapeutic targets for SLE.

Introduction

Systemic lupus erythematosus (SLE), also known as lupus, is a rare systemic autoimmune disease that mostly affects middle-aged women, mainly of Asian, African, American, and Hispanic origin (Costa-Reis and Sullivan, 2013; Cui et al., 2013; Gurevitz et al., 2013). SLE affects an estimated 5 million people across the world, with an incidence of 1–10 per 100,000 person-years (Pons-Estel et al., 2010). SLE is characterized by a wide range of different autoantibodies, deposition of immune complexes, and immune system infiltration and inflammation within damaged organs. SLE autoantibodies invade the patient’s kidneys, heart, skin, joints, and brain, leading to various typical clinical symptoms. The most common clinical symptoms of lupus are rash, arthritis, and fatigue. Severe complications of SLE lead to nephritis, anemia, neurological symptoms, and thrombocytopenia, eventually leading to severe morbidity and mortality.

SLE is characterized by its clinical heterogeneity, with a wide range of clinical manifestations reflecting its complex etiopathogenesis (Tan et al., 1982). The clinical heterogeneity of SLE highlights the contribution of genetic and environmental factors to the susceptibility to the disease (Prokunina and Alarcon-Riquelme, 2004; Harley et al., 2009; Yang and Lau, 2015; Dang et al., 2016; Wang et al., 2017). To date, the reason for phenotypic variation in SLE is unknown. Understanding the molecular mechanisms behind the pathogenesis of SLE phenotypes could help in developing more efficient therapeutic approaches and preventive strategies.

With the extensive use of gene detection methods, high-throughput sequencing and extensive microarray data profiling studies on SLE have been conducted, and several differentially expressed genes (DEGs) and cellular pathways in SLE have been identified (Borrebaeck et al., 2014; Zhu et al., 2015). Nevertheless, until now, no particular gene has been recognized to act as a potential marker for the diagnosis of SLE. In addition, a large amount of data obtained from microarray technology and high-throughput sequencing have not been fully used. Ducreux et al. (2016) collected blood samples from SLE patients and healthy volunteers to identify differentially expressed genes (Ducreux et al., 2016). However, the interactions among differentially expressed genes and key genes involved in the signaling pathways of SLE remain to be elucidated. In addition, previous studies of genetic factors primarily focused on single genes; nevertheless, interactions among multiple genes may result in the multisystem invasion characteristics observed in SLE (Smith et al., 2017). Remarkably, studies have shown that disease−associated gene expression networks have a potential role in the immune response, which highlights their mechanism and therapeutic value for SLE (Deng and Tsao, 2010; Bentham et al., 2015).

Integrating and reanalyzing the data using bioinformatics methods may help in identifying gene regulatory pathways, essential genes, and their associated networks in SLE disease, which can provide new and valuable ideas for understanding the molecular mechanisms and identifying reliable diagnostic and therapeutic targets of SLE. Therefore, in this study, we first conducted a comprehensive collection of genes associated with SLE from the GEO dataset with ID GSE30153. Then, we performed a bioinformatics analysis of these genes with the MCODE (Molecular Complex Detection), GeneGo, and ClueGO tools. To further explore the pathogenesis of SLE in a more specific manner, functions and pathways identified by the modules were used to indicate the biological processes and biochemical pathways related to the immune system. Finally, the genes potentially associated with arthritis, pleurisy, and myocarditis, which are the common complications of SLE, were compared with SLE-related genes to identify common genes that participated in the development of SLE. To interpret the biological relevance of these changes in gene expression, we analyzed the microarray data via an integrated bioinformatic analysis expanding on traditional microarray analysis methods, namely, Gene Ontology (GO) and pathway analysis, thereby allowing the construction of interaction networks that might identify novel prognostic markers and therapeutic targets.

Materials and Methods

Acquisition of Array Data and Processing

Gene expression profiling data from microarray array analysis of the GSE30153 dataset were downloaded from the NCBI GEO database (Gene Expression Omnibus database)[1]. The database accommodates gene expression datasets from a variety of experiments, such as DNA-seq, ChIPs, RNA-seq, microarray, and high-throughput hybridization array (Edgar et al., 2002; Barrett et al., 2013). GSE30153 contains 26 samples, including 17 patients with SLE and 9 healthy controls of human sorted B-cells obtained by using the platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array (Garaud et al., 2011). The downloaded gene expression profiling data are freely available in the public database, and there were no human or animal experiments conducted by any of the authors in this study.

Preprocessing of Data and DEG Identification

Using the robust multiarray standard model, the initial information from the dataset was subjected to quantile normalization, background correction, and log transition (Irizarry et al., 2003). Preprocessing included changing to gene symbols from probe IDs using the Gene ID converter from Entrez (Alibés et al., 2007). The statistical online tool GEO2R uses the R/Bioconductor, and limma package v3.26.8 was used to screen the raw gene expression data (Smyth, 2005; Barrett et al., 2013; Ritchie et al., 2015). We performed a Benjamini–Hochberg test (to determine the false discovery rate) and T-tests to compute the false discovery rate (FDR) and p-values to identify the DEGs between SLE patients and healthy control human sorted B-cells (Benjamini and Hochberg, 1995; Aubert et al., 2004). We set the primary criteria of | log (2 fold change) | > 1 and p < 0.05 to obtain significant DEGs from the dataset, whereas cutoffs of log2FC ≥ 1 and log2FC ≤ −1 were used to denote upregulated and downregulated DEGs, respectively. For high-throughput sequencing, a logarithm to base 2 is widely used and in the initial scaling, the doubling is equivalent to a log2FC of 1 (Love et al., 2014). A volcano plot was constructed using a web-based tool[2]. The resulting DEGs were used for further analysis.

Constructing PPI Networks

To assess the relationships between the DEGs from the GSE30153 dataset, we constructed a protein–protein interaction (PPI) network by using Search Tool for the Retrieval of Interacting Genes (STRING v11.0)[3] (Szklarczyk et al., 2017, 2019). The cutoff criterion was set to a high confident interaction score of ≥0.7 to eliminate inconsistent PPIs from the dataset. We then incorporated the results from the STRING database into Cytoscape software (v3.7.2)[4] to envisage the PPIs within the statistically relevant DEGs (Shannon et al., 2003). The MCODE plugin from Cytoscape was utilized to identify the interconnected regions or clusters from the PPI network. The cluster finding parameters were adopted, such as a degree cutoff of 2, a node score cutoff of 0.2, a kappa score (K-core) of 5, and a max depth of 100, which limits the cluster size for coexpressing networks (Bader and Hogue, 2003). The top clusters from MCODE were subjected to ClueGO v2.5.5/CluePedia v1.5.5 analysis to obtain comprehensive GO and pathway results from the PPI network. ClueGO combines GO and pathway analyses from KEGG and BioCarta and provides a fundamentally structured GO or pathway network from the PPI network (Bindea et al., 2009).

Metacore GeneGo Analysis of DEGs

Metacore, a Cortellis Solution software (Clarivate Analytics, London, United Kingdom)[5], was used to perform curated pathway enrichment analysis and GO analysis. GeneGo facilitates the rapid assessment of metabolic pathways, protein biological networks, and pathway maps from high-throughput experimental data (MetaCoreLogin | Clarivate Analytics). Based on a significance threshold of p < 0.05, a pictorial representation of the molecular interactions of DEGs from the study groups is generated. Determination of a hypergeometric p-value enables the estimation of the chance that an intersection between DEGs and ontological elements is random. An FDR < 0.05 was used as a criterion to calculate if statistically significant DEGs constituted a processor pathway.

Results

Identification of DEGs From the Dataset

Our study contained the gene expression profiles of the GSE30153 dataset from the GEO database, which were submitted by Garaud et al. (2011) based on analysis with the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) (Garaud et al., 2011). The dataset encompasses 26 samples, including 17 patients with SLE and 9 healthy controls (Table 1). By utilizing the GEO2R online tool, we obtained the differentially expressed genes (DEGs) from the GSE30153 dataset by comparing the SLE samples with control samples. By calculating p-values and | log2FC | values, the top 250 DEGs were identified. A volcano plot was constructed using the Rstudio web server ShinyVolcanoPlot to identify DEGs by comparing the SLE and control groups from the dataset. The volcano plot in Figure 1 depicts all the DEGs with a log2FC against the – log10 (p-value) between the two groups. With cutoffs of p < 0.05 and log2FC ≥ 1.0 or ≤−1, we found 4 and 13 genes that were upregulated and downregulated, respectively, between the two groups (Table 2). The genes that were differentially expressed between the two groups are shown in Supplementary Table S1.

TABLE 1

The primary characteristics of 26 studies in GSE30153 procured from the Gene Omnibus Expression database.

Group	Accession	Title	Organism	Disease state	Tissue	Cell type
Patient	GSM746726	Patient 1	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746727	Patient 2	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746728	Patient 3	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746729	Patient 4	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746730	Patient 5	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746731	Patient 6	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746732	Patient 7	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746733	Patient 8	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746734	Patient 9	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746735	Patient 10	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746736	Patient 11	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746737	Patient 12	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746738	Patient 13	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746739	Patient 14	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746740	Patient 15	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746741	Patient 16	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
	GSM746742	Patient 17	Homo sapiens	Systemic lupus erythematosus (SLE)	Blood	Human sorted B cell
Control	GSM746743	Control 1	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746744	Control 2	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746745	Control 3	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746746	Control 4	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746747	Control 5	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746748	Control 6	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746750	Control 8	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746751	Control 9	Homo sapiens	Control	Blood	Human sorted B cell
	GSM746752	Control 10	Homo sapiens	Control	Blood	Human sorted B cell

FIGURE 1

Pictorial representation of volcano plot for differentially expressed genes (DEGs) in systemic lupus erythematosus (SLE) compared to controls from the GSE30153 dataset. The X-axis represents Log2FC, large magnitude fold changes; Y-axis represents −log10 of a p-value, high statistical significance. Each black dot represents one gene. Black dots above red and beside blue line (left-sided and right-sided) are log2FC ≥ 1 and p-value <0.05, representing SLE related DEGs.

TABLE 2

Significantly upregulated and downregulated DEGs between two groups from GSE30153 dataset are tabulated.

GENE SYMBOL	log2FC	p-value
Upregulating Genes
EGR1	1.22	0.00074
DSE	1.125	0.00291
CD1C	1.068	0.00053
GPM6A	1.052	0.00097
GPM6A*	1.043	0.002981
Downregulating Genes
RRM2	–2.406	0.0027527
RRM2*	–2.152	0.0030096
TYMS	–1.923	0.0032032
CD38	–1.702	0.0031747
CAV1	–1.516	0.0048324
MIR7110	–1.4	0.0027212
ELL2	–1.354	0.0035256
SLC44A1	–1.219	0.0035298
SAR1B	–1.176	0.0049953
MAN1A1	–1.111	0.004425
CHAC2	–1.071	0.004354
ERAP1	–1.047	0.005321
ARF4	–1.044	0.0058274
PDIA4	–1.014	0.0041988

Screening of Module and Construction of Interlinking PPI Network

To assess the protein–protein connections among the DEGs, we used the STRING tool to compute the protein interactions and plotted them using Cytoscape v3.7.2. Figure 2 depicts the PPI network with 103 nodes and 201 edges. The DEGs are represented as nodes, and the edges are interactions between the DEGs. A combined node score of >0.4 was considered to be significant. MCODE plugin v1.5.1 from Cytoscape was utilized to identify the densely interlinked regions within the protein network. As a result, we obtained the top two significant clusters from the DEGs protein network with MCODE scores of 5.043 and 3.625. A graphical representation of these clusters is shown in Figures 3A,B. The subnetworks, scores, number of nodes and edges, and node IDs are tabulated in Table 3.

FIGURE 2

The network demonstrates the protein–protein interaction between the DEGs identified from GSE30153 using Cytoscape. The nodes represented as ellipse (robin’s blue) and edges as lines (gray).

FIGURE 3

The MCODE (Molecular Complex Detection) plugin from Cytoscape analyzed the top two clusters derived from the network of interactions between protein and protein. (A) Cluster 1; (B) Cluster 2. The MCODE cluster score > 3. The nodes represented as ellipse (green) and edges as lines (gray).

TABLE 3

The interconnected regions are clustered from the GSE30153 dataset using MCODE plugin in Cytoscape.

Cluster	Score (density × No. of nodes)	Nodes	Edges	Node IDs
1	5.043	45	116	OBBP2, IL6ST, CD1C, ABHD15, PTPRJ, ITGAX, UNC5B, TLR10, CD38, GAS6, NCF4, MAPK9, DDAH2, PTPN6, GAB1, ARHGAP24, CUL3, PROX1, CYTH4, E2F6, TNFSF9, VEGFA, TYMS, IL7, RRM2, PRKCB, MPEG1, MARCKS, SLC2A5, ARHGAP35, BMP6, TCF3, AKT1, EIF4EBP2, GNG11, CAV1, FYN, EGR1, SIGLEC10, CD24, CHEK1, E2F7, CD84, CDK6, SRGN
2	3.625	15	29	RRM2, CKAP2, MGLL, TCF3, SUB1, EGR1, POLA2, RPA1, CHEK1, E2F7, CASC5, DP2, E2F6, CDK6, TYMS

Enrichment Analysis by ClueGO

The top two subnetworks from MCODE were used as an input for analyzing the functional enrichment of PPI subnetworks using the ClueGO/CluePedia plugin from Cytoscape. In a biologically clustered subnetwork, ClueGO helps to visualize the biological terms of broad gene clusters. The subnetwork enrichment analyses of MCODE cluster 1 and cluster 2 are depicted in Figures 4A,B. For functional enrichment analysis, we set the statistical options based on a two-sided hypergeometric test with a Benjamini–Hochberg correction, p ≤ 0.05, and kappa scores ≥ 0.4 as criteria. The DEGs from cluster 1 were shown to be enriched mostly in T-cell costimulation (GO: 0031295), lymphocyte costimulation (GO: 0031294), negative regulation of vascular permeability (GO: 0043116), the metaphase/anaphase transition of the mitotic cell cycle (GO: 0007091), regulation of the transcription involved in the G1/S transition of the mitotic cell cycle (GO: 0000083), negative regulation of signal transduction in the absence of ligand (GO: 1901099), and KEGG pathways such as hematopoietic cell lineage (KEGG: 04640), B-cell receptor signaling pathway (KEGG: 04662), ErbB signaling pathway (KEGG: 04012), and AGE-RAGE (advanced glycation end products and receptor for AGE) signaling pathway in diabetic complications (Figure 4A). The DEGs from cluster 2 were mainly enriched in the regulation of the transcription involved in the G1/S transition of the mitotic cell cycle (GO: 0000083), the negative regulation of the G0 and G1 transitions (GO: 0070317), and the p53 signaling pathway (KEGG: 04115) (Figure 4B). The pathways that were activated in the enrichment analysis were highly related to B-cell pathophysiology, resulting in events associated with the immune system, vasculopathy, and kidney.

FIGURE 4

Visualization of Gene Ontology (GO) enrichment profiles from DEGs using Cytoscape software based on network analysis of ClueGO/CluePedia inferred from MCODE cluster 1 (A) and cluster 2 (B). The plugin provides a combined enrichment analysis of clusters, including the GO biological process, molecular function, and pathway from KEGG. The GO term/pathway network connectivity defined by edges and functional clusters on genes shared between terms (kappa score ≥ 0.4) and displaying pathways only with p≤ 0.05. The size of the node indicates the p-value. The color code of nodes represents the functional group that they belong to. The most important functional terms specify the pathway names within each class are indicated in bold colored characters. (A) The network enrichment analysis of cluster 1. Each node constitutes a precise term for cluster 1; (B) The network enrichment analysis of cluster 2. Each node constitutes a precise term for cluster 2.

MetacoreTM GeneGo for Enrichment Analysis of DEGs

Further functional enrichment analysis was carried out using MetacoreTM GeneGo software from Clarivate Analytics to comprehensively dissect the pathways associated with the DEGs. Using the functional ontology feature in GeneGo, the IDs of potential genes that were involved in the target pathways were identified. Based on hypergeometric p-values, the probability that the intersection of a gene set and associated ontological objects was random was evaluated. A decreased p-value indicated that the entity would be more significant to the DEGs, suggesting a better score. The functional enrichment analysis of the DEGs defined the top 10 metabolic networks, and canonical pathway maps are depicted in Figures 5A,B. For each classification, the significant statistical data rely on a low p-value. The pathway maps with the lowest p-value are shown in Figures 6A–C. These are the top-scoring signaling pathways based on the gene enrichment distribution, which emphasizes that the DEGs from human sorted B-cells are triggered via oxidative stress and ROS-induced cellular signaling (Figure 6A), chemotaxis and lysophosphatidic acid signaling via GPCRs (Figure 6B), and androgen receptor activation and downstream signaling in prostate cancer (Figure 6C). The well-distinguished proteins and complexes of proteins are shown as specific symbols[6]; all experimental data are displayed and have corresponding thermometer-like symbols on all the maps. The upregulated genes are indicated by a red thermometer facing upwards.

FIGURE 5

The top 10 metabolic networks and pathway maps were annotated using GeneGo enrichment analysis for the genes that are differentially expressed from SLE patients vs. healthy controls, respectively. (A) The content of these metabolic networks was annotated and defined by GeneGo Cortellis Solution software. Each process represents a pre-set network of protein interactions characteristic for the process, and sorting was performed for the metabolic networks that are statistically significant. (B) The pathway maps (canonical) of GeneGo display a series of signals and metabolic charts that cover human in a structured manner. The significant expression of a gene/protein represented in histogram height.

FIGURE 6

The enrichment analysis from GeneGo showed three regulated pathways with the highest score that are triggered in the SLE human sorted B-cells. (A) Oxidative stress ROS induced cellular signaling. (B) Chemotaxis lysophosphatidic acid signaling via G protein-coupled receptors (GPCRs). (C) Androgen receptor activation and downstream signaling in prostate cancer. The image depicts the protein and protein complexes that are well characterized as a specific symbol; laboratory data from all reports are correlated and shown on the maps as thermometer-like indicators. The red or blue color upward/downward thermometers indicate gene transcripts with upregulation/downregulation, respectively. The proteins connected by arrows demonstrate the stimulating and inhibitory effect of the protein. Further details are given at https://portal.genego.com/help/MC_legend.pdf.

Pathway Map Interaction Results From Clarivate

From the MetacoreTM results, we extracted the key genes from the enriched pathways that were differentially expressed, such as EGR1, CD38, CAV1, and AKT1. The differential expression of these genes was involved in the activation or inhibition of specific protein complexes in the enriched pathway maps (Figure 6 and Table 4). Early growth response 1 (EGR1) is a transcription factor that interacts with the IGF-2, APEX, SRD5A1, CD44, and EGFR genes and activates them through transcriptional regulation. The cyclic adenosine diphosphate (ADP) ribose hydrolase CD38 is an enzyme involved in the activation of the genes Semaphorin 4D, CD19, and c-Cbl through physical interactions. Caveolin 1 (CAV1) is a binding protein shown to interact with the ErbB2, MDR1, HTR2A, and androgen receptor genes, and inhibition or activation is followed by specific binding to its corresponding proteins. RAC-alpha serine/threonine-protein kinase (AKT1) is a protein kinase that interacts with the FKHR, mTOR, Bcl-10, FOXO3A, HNF3-beta, and GSK3 beta genes via phosphorylation, resulting in inhibition or activation. These genes were differentially expressed between sorted B-cells from controls and sorted B-cells from SLE patients and result in transcriptional regulation and inhibition of genes/proteins within the top-scored pathway maps.

TABLE 4

The interaction reports of key genes from pathway maps by Clarivate Analytics.

Network object “from”	Object type	Network object “to”	Object type	Effect	Mechanism	Link info	Input IDs	Signal	P-value	PMID
EGR1	Transcription factor	IGF-2, APEX, SRD5A1, CD44, EGFR	Receptor ligand, generic enzyme, generic enzyme, generic receptor, receptor with enzyme activity	Activation	Transcription regulation	EGR1 increases IGF II expression, EGR1 binds to gene APEX promoter and activates APEX expression, Egr-1 trans-activates the 5alpha-R1 promoter via the Egr-1-binding site at position −60/−54, Putative EGR1 binding site is found in gene CD44 promoter, EGR1 binds to gene EGFR promoter and activates EGFR expression.	EGR1	1	0.0032807	8584025; 9925986; 10606246; 11336542; 16043101; 19276347; 29092905; 29170465; 15788231; 15936112; 17194527; 18215136; 8628295; 9300687; 12670907; 15155664; 15923644; 19195913; 20357818; 25673149; 1417865; 11830539; 16750517; 17230532; 19032775; 20190820; 23763269
CD38	Generic enzyme	SEMA4D, CD19, c-Cbl	Generic receptor, generic binding protein, generic enzyme	Activation	Unspecified, Binding	CD31-induced activation of CD38 up-regulates Semaphorin 4D cell-surface expression in B cells, CD19/CD81 complex interacts with CD38 but this interaction is not required to induce proliferation in mouse B-lymphocytes, Fluorescence resource energy transfer and coimmunoprecipitation showed that c-Cbl and CD38 bind each other.	CD38	1	0.0031747	15613544; 17327405; 20570673; 22564057; 8695807; 18974118; 19635790
CAV1	Generic binding protein	ErbB2, MDR1, HTR2A, Androgen receptor	Receptor with enzyme activity, transporter, GPCR, transcription factor	Unspecified, Inhibition, activation	Binding	HER2 physically interacts with caveolin-1, Caveolin-1 interacts with p-gp, Down-regulation of caveolin-1 by siRNA reduced the interaction between p-gp and caveolin-1, followed by a decrease in [3H]-Taxol and [3H]-Vinblastine accumulation in RBE4 cells, Caveolin-1 physically interacts with HTR2A and increases its activity, Highly conserved 9 amino acid motif in the ligand binding domains (E domains) was identified in human/mouse ER alpha and ER beta, progesterone receptors A and B, and the androgen receptor. The localization sequence mediated palmitoylation of each SR, which facilitated caveolin-1 association, subsequent membrane localization, and steroid signaling.	CAV1	1	0.0048324	9374534; 9685399; 11697880; 22389470; 14622130; 15239129; 15498565; 17326770; 18485890; 19099191; 22389470; 25788263; 15190056; 8703009; 11278309; 17535799; 17940184; 18786521; 19931639; 22771325; 24375805
AKT1, AKT (PKB)	Protein kinase	FKHR, mTOR, Bcl-10, FOXO3A, HNF3-beta, GSK3 beta	Transcription factor, protein kinase, generic binding protein	Inhibition, activation	Phosphorylation	AKT1 phosphorylates FKHR1 and decreases its activity, Increased AKT phosphorylation regulates different metabolic pathways in liver, including increases in protein synthesis through activation of mTOR/p70 (S6kinase), AKT1 phosphorylates Bcl-10 and increases its activity, AKT1 phosphorylates FOXO3A and decreases its activity, AKT1 phosphorylates HNF3-beta and decreases its activity, AKT (PKB) inhibits GSK3 alpha by phosphorylation at Ser-9.	AKT1	1	0.0010146	10102273; 10358014; 10358075; 10377430; 11030146; 12393870; 16076959; 16099987; 16230533; 16603397; 17186497; 18388859; 18391970; 18420577; 18687691; 18786403; 19703413; 20940043; 21106439; 21157483; 21238503; 21407213; 21440577; 21708191; 21779512; 26053093; 27966458; 30413788; 10567225; 10910062; 11357143; 11438723; 12767043; 14970221; 15208671; 15549092; 16818631; 17660512; 18505677; 18566586; 18566587; 18678273; 21097843; 21177249; 21302298; 21343617; 22084251; 22595285; 23686889; 23872070; 26958938; 29221131; 16280327; 10102273; 12130673; 12767043; 17570479; 17577629; 17957242; 17960591; 18391970; 18687691; 19703413; 20223831; 20399660; 21106439; 21157483; 21440577; 21621563; 21708191; 21775285; 21779512; 24518891; 27966458; 14500912; 11584303; 11701324; 12124352; 12750378; 12808085; 14966899; 14985354; 15016802; 15297258

Discussion

DNA microarrays and next generation sequencing (NGS) approaches are high-throughput technologies that have resulted in the emergence of new biomedical discoveries. Data from microarray and gene expression profiles have enabled a deeper understanding of the intrinsic molecular pathways of complex mechanisms of biological systems and their responses (Russo et al., 2003; Babu, 2004; Perez-Diez et al., 2013; Kumar et al., 2019). It is therefore highly relevant to examine the peripheral B-cell transcriptomes of SLE patients and healthy controls to determine genes that are differentially regulated and their target pathways. Our current study extracted DEGs from 17 SLE patients and 9 healthy controls from the GEO database (GSE30153) (Garaud et al., 2011). The top 250 DEGs were identified, including 4 upregulated and 13 downregulated genes from the groups through bioinformatics strategies (Table 2 and Supplementary Table S1). These identified DEGs were subjected to ClueGO and GeneGo MetacoreTM analysis for GO and pathway annotation, and constructed the interacting networks of PPI and used for cluster analysis. In the network, the nodes were considered proteins, and the edges were their interactions. Using network topology features, the PPI network can be analyzed to distinguish the core proteins that are involved in the pathways (Barabási and Oltvai, 2004; Ideker and Sharan, 2008; Keskin et al., 2016; Kumar et al., 2019). The identified DEGs from the present study were analyzed with STRING to exploit the complex interactions between the DEGs via text mining, evidence from experiments, and repositories (Figure 2). We performed module screening of the PPI networks using the MCODE plugin from Cytoscape. As a result, we obtained significant clusters that are densely interlinked regions in the PPI network (Figures 3A,B). Screening of these clusters from the network might help to identify the essential genes that are involved in the pathogenesis and progression of SLE. The obtained clusters mostly contained protein complexes or proteins present in the pathways in the PPI network, and cluster visualization is essential for comprehending the properties of the network functionally and systematically (Krogan et al., 2006; Rahman et al., 2013).

Furthermore, to identify the functional enrichment of these subnetworks from MCODE, we implemented the ClueGO plugin for analysis. This revealed that the DEGs were enriched in most essential pathways, which are highly associated with the immune system. The GO and KEGG enrichment analyses of the DEGs from cluster 1 showed that they were mostly enriched in T-cell costimulation, lymphocyte costimulation, negative regulation of vascular permeability, the metaphase/anaphase transition of the mitotic cell cycle, regulation of the transcription involved in the G1/S transition of mitotic cell cycle, the hematopoietic cell lineage, the B-cell receptor signaling pathway, the ErbB signaling pathway, the AGE-RAGE signaling pathway in diabetic complications, and pancreatic cancer. Interestingly, the costimulation of T-cell and lymphocyte receptors is recognized to be important in SLE pathogenesis by enabling communicating with B-cells for the production of autoantibodies (Shlomchik et al., 2001; Mak and Kow, 2014). In SLE, negative regulation of vascular permeability may be induced by different mechanisms; the dysregulated genes from the cluster 1 subnetwork might lead to endothelial cell damage and vasculopathy (Favero et al., 2014; Lee et al., 2019). The differential cell signaling results in the recruitment of various proteins and inappropriate activation of B-cells (Zhou et al., 2009; Comte et al., 2015). Oxidative stress is common in inflammatory disorders and results in the increased production of reactive carbonyl groups that are partially converted to AGEs, and the DEGs in the AGE-RAGE signaling pathway might also be involved in the accumulation of AGEs in SLE patients and lead to diabetic complications (de Leeuw et al., 2007; Li et al., 2007; Kurien and Scofield, 2008; Nienhuis et al., 2008). Interestingly, our enrichment analysis found that the identified differential expression of the genes (AKT1, VEGFA, CDK6, and MAPK9) that were involved in the risk of developing pancreatic cancer in SLE patients was due to chronic inflammation, suggesting that these genes might be involved in the pathogenesis of SLE. Our findings are therefore consistent with the roles of genes that are differentially expressed in SLE-causing pathways (Figure 4A). The enrichment analysis of the cluster 2 subnetwork showed that the DEGs were mostly enriched in the regulation of transcription involved in the G1/S transition of the mitotic cell cycle, the negative regulation of the G0 and G1 transitions, and the p53 signaling pathway. It has been reported that the proliferation of T-cells is followed by lowered levels of cyclin-dependent kinase (CDK) inhibitors, and alterations in the expression of CDKs in the G0/G1 phase were seen in the lymphocytes of SLE patients (Yamauchi and Bloom, 1997; Tang et al., 2009). The DEGs involved in the cluster 2 subnetwork might negatively regulate these pathways. Alterations in cyclin-CDK complex behavior and cyclin-dependent kinase inhibitors (CDKIs) have been reported to alter the proliferation of T-cells, oxidative stress, and immune responses (Santiago-Raber et al., 2001; Tang et al., 2009). p53 signaling is essential for various cellular mechanisms, and defects in this signaling pathway are associated with SLE development. Considerably elevated levels of p53 protein are found in SLE patients with active inflammatory disorders (Miret et al., 2003; Veeranki and Choubey, 2010). Apoptosis dysregulation appears to be another cause of SLE pathogenesis because the possible sources of autoantigens are cell debris from apoptosis in SLE, and excessive cellular senescence of the immune cells, especially T-cells, was reported in SLE patients with peripheral blood mononuclear cells (PBMCs) and skin lesions (Colonna et al., 2014; Sáenz-Corral et al., 2015). Thus, our identified DEGs (RRM2, APC, CHEK1, E2F6, TYMS, E2F7, and CDK6) from the cluster 2 subnetwork are highly related to and consistent with the members of the signaling pathways associated with the immune system, apoptosis, the cell cycle, and vasculopathy.

To clearly define the interactions between the proteins and signaling pathways examined from the interpretation of STRING, Cytoscape, MCODE, and ClueGO analyses, we implemented the GeneGo Metacore software, which incorporates extensive data on metabolic signaling pathways and their regulatory mechanisms and contains accurately complied networks of biological systems. By utilizing the GeneGo Metacore software, we obtained a detailed description of the DEGs that participate in SLE pathogenesis based on the determined p-values. Among the top 10 metabolic networks, the phosphatidylinositol-3,4,5-triphosphate pathway, O-hexadecanoyl-(L)-carnitine pathway, 1,2-didocosapentaenoyl-sn-glycerol 3-phosphate pathway, and 1-linoleoyl-glycerol-3-phosphate pathway were profoundly enriched and significant in the SLE DEGs (Figure 5A). The increased activity of phosphatidylinositol-3,4,5-triphosphate stimulates essential cell signaling pathways such as the pathways involved in cell division, survival, and the rapid increase in T-lymphocytes in SLE (Comte et al., 2015). PI3K (phosphatidylinositol 3-kinase) is a protein kinase that phosphorylates phosphatidylinositol 4,5-phosphate to regulate the signaling of T-lymphocytes; an increased amount of PI3K was also observed in an animal model of lupus (Liu et al., 1998; Grolleau et al., 2000; Niculescu et al., 2003; Joseph et al., 2014). Modification of the carnitine signaling pathway results in various organ failures by producing effective responses to pathogens (Famularo and De Simone, 1995; Famularo et al., 2004). Thus, the DEGs involved in the O-hexadecanoyl-(L)-carnitine pathway might lead to increased immune responses. In addition, the top three pathways associated with the DEGs of sorted B-cells from SLE patients were mostly enriched in oxidative stress- and ROS-induced cellular signaling (Figure 6A), chemotaxis and lysophosphatidic acid signaling via GPCRs (Figure 6B), and androgen receptor activation and downstream signaling in prostate cancer (Figure 6C). Recent findings have shown that oxidative stress and ROS induce molecular alterations that have adverse effects in SLE patients (Choi et al., 2016; Tsokos et al., 2016; Lightfoot et al., 2017). Elevated oxidative stress in SLE patients leads to the accumulation of higher amounts of oxidative lipoproteins, which are harmful in zebrafish models and cause additional oxidative damage to the system (Chung et al., 2007; Park et al., 2016; Lightfoot et al., 2017). Interestingly, our study identified the EGR1 gene as downregulated in the SLE patients in comparison to controls, and it also plays a role in ROS signaling. This clearly indicates that EGR1 might be required to maintain the oxidative stress and ROS signaling pathways.

Moreover, the DEGs involved in the oxidative stress signaling pathway might contribute to peripheral neuropathy, damage to blood vessels, and cardiovascular events, which are the prominent clinical conditions found in SLE patients. Chemotaxis and lysophosphatidic acid (LPA) signaling are essential pathways in autoimmune inflammatory disorders, and GPCRs are responsible for regulating immune cells via LPA receptors (Yang et al., 2005; Skoura and Hla, 2009). G2A gene knockout resulted in the hyperresponsiveness of T-cells to T-cell receptor stimulation, manifesting as an increased proliferation of T-cells, which may promote inflammatory phenotypes in G2A-deficient mice (Le et al., 2001). Various studies have suggested that LPA plays a vital role in atherosclerosis progression and development by promoting neutrophil and monocyte adherence and enhancing inflammatory events (Siess et al., 1999; Smyth et al., 2008; Skoura and Hla, 2009). The androgen receptor (AR) is a transcription factor that is activated by a ligand and is essential for cells targeted by the androgen response (Robeva et al., 2013; Gubbels Bupp and Jorgensen, 2018). AR also regulates immune function in SLE via transcriptional regulation of various genes. Our study identified the transcription factor AR, which positively regulates the c-Myc, SCAP, prosaposin, and KLF5 genes, which are responsible for inflammatory responses, and promotes tumor growth factors and cytokine signaling when activated (Figure 6C). The enriched terms from ClueGO modules and the GeneGo-identified terms correlated well in this study and validate the significance of the findings from the pathway maps. The combined results from these two enrichment analyses suggest that B-cells from SLE patients and B-cells from healthy controls undergo differential gene expression associated with positive regulation of kidney development, the hematopoietic cell lineage, positive regulation of vasoconstriction, T-cell costimulation, and regulation of the transcription involved in the G1/S transition of the mitotic cell cycle.

Furthermore, the interaction results from the GeneGo analysis provided the essential genes (EGR1, CD38, CAV1, and AKT1) from the pathway maps constructed from the DEGs. Among them, EGR1 (early growth response 1) is a transcription factor shown to interact with the IGF-2 (insulin-like growth factor 2), APEX (apurinic/apyrimidinic endodeoxyribonuclease 1), SRD5A1 (steroid 5 alpha-reductase 1), CD44 (cell surface glycoprotein CD44), and EGFR (epidermal growth factor receptor) genes and transcriptionally regulate them by activating or promoting their expression in sorted B-cells from patients with SLE (Liu et al., 1995; Recio and Merlino, 2003; Lee et al., 2005; Pines et al., 2005; Blanchard et al., 2007; Rui et al., 2008; Cullen et al., 2010; Sauer et al., 2010). The cyclic ADP ribose hydrolase (CD38) is also known as cluster of differentiation 38 protein, can be found on several immune cells, and activates SEMA4D (semaphorin-4D or cluster of differentiation 100), CD19 (B-lymphocyte antigen CD19 or cluster of differentiation 19), and c-Cbl (Casitas B-lineage lymphoma proto-oncogene) (Deaglio et al., 2005, 2007; Shen and Yen, 2008; Vences-Catalán et al., 2012). These interactions with CD38 result in the activation of B-lymphocytes and increase immune responses in SLE patients. The protein caveolin 1 (CAV1) has been shown to interact with the ErbB2 (Erb-B2 receptor tyrosine kinase 2), MDR1 (multidrug resistance protein 1), HTR2A (5-hydroxytryptamine receptor 2A), and AR (androgen receptor) genes. Several studies have suggested that caveolin 1 physically interacts with HER2, p-gp, HTR2A, and AR and activates/inhibits them by binding to their specific caveolin-binding motif (Couet et al., 1997; Lu et al., 2001; Razani and Lisanti, 2001; Bhatnagar et al., 2004; Bennett et al., 2010, 2014; Yu et al., 2012). AKT1 is a protein kinase that interacts with FKHR (Forkhead box protein O1), mTOR (mechanistic target of rapamycin), Bcl-10 (B-cell lymphoma/leukemia 10), FOXO3A (Forkhead box O3), HNF3-beta (hepatocyte nuclear factor 3-beta), and GSK3 beta (glycogen synthase kinase three beta). The protein kinase AKT1 inhibits FKHR via phosphorylation and decreases its activity (Biggs et al., 1999; Rena et al., 1999; Tang et al., 1999; Hay, 2011), whereas it increases its activity via phosphorylation of mTOR (Navé et al., 1999; Sekuliæ et al., 2000; Ikenoue et al., 2008; Thirumal Kumar et al., 2019). Additionally, AKT1 phosphorylates Bcl-10 at the specific residues Ser231 and Ser218, increasing its activity (Yeh et al., 2006), while it inhibits the action of FOXO3A via phosphorylation and decreases its activity, increasing the survival of cells (Brunet et al., 1999; Linding et al., 2007; Calnan and Brunet, 2008; Li et al., 2008; Tzivion et al., 2011). AKT1 decreases HNF3-beta activity by phosphorylating it at Thr156 (Wolfrum et al., 2003), whereas phosphorylation of GSK3-beta by AKT1 occurs at Ser9 to inhibit its activity (Brazil and Hemmings, 2001; Salas et al., 2004; Kuemmerle, 2005; Shin et al., 2006; Markou et al., 2008). This suggests the vital genes we identified from the DEGs of patients with SLE play essential roles in the development and progression of SLE via different signaling pathways to increase autoimmune responses.

In addition to the interaction analysis, we carried out interrelation analysis for the essential genes to determine the relationships between the genes, which implicitly or explicitly interacted with each other. Interestingly, the identified genes indirectly communicated with each other via molecular signaling pathways related to mTOR signaling, apoptosis, PI3K-Akt signaling, the hematopoietic cell lineage, positive regulation of vasoconstriction, signaling by receptor tyrosine kinases, AGE-RAGE signaling, and lymphocyte and T-cell costimulation (Figure 7). EGR1 and AKT1 are directly involved in oxidative stress via ROS and AGE-RAGE signaling, whereas CAV1 is directly involved in tyrosine kinase receptor signaling and lymphocyte and T-cell costimulation. CD38 is directly associated with the hematopoietic cell lineage and positive regulation of vasoconstriction. Overall, the dysregulation of the indicated pathways in SLE patients is a result of differential gene expression. The essential genes are differentially expressed between cells from patients with SLE and cells from healthy controls and are present in important signaling cascades, which could be a crucial factor for SLE development.

FIGURE 7

The interrelation analysis of genes EGR1, CD38, CAV1, and AKT1 that strongly associated to SLE. Each gene involved in different pathways via interacting to each other. Inbuilt color code was provided to all the genes based on the STRING tool from Cytoscape.

Conclusion

Taken together, the results of our comprehensive bioinformatics analysis showed that the DEGs identified between sorted B-cells from patients with SLE sorted B-cells from controls could play a significant role in the growth, progression, and development of SLE. This study identified 4 upregulated and 13 downregulated genes, including essential genes (EGR1, CD38, CAV1, and AKT1), from the pathway enrichment analysis. Indeed, the identified pathways from the enrichment analysis were strongly related to the immune system, vasculopathy, cardiovascular functions, and inflammatory responses, which are processes that can lead to the development of SLE. The broad understanding of SLE pathophysiology from this study will allow us to identify and develop therapies targeting SLE and contribute to personalized treatment strategies. Collectively, the study findings could aid in enhancing our understanding of the fundamental molecular processes of SLE and provide possible strategies for early diagnosis in SLE; in addition, combinatorial therapeutic strategies using oxidative stress and ROS cellular signaling and lysophosphatidic acid signaling via GPCRs might have symbiotic effects on the molecular events in SLE.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30153.

Author Contributions

SU, DT, HZ, and CG were involved in the design of the study and the acquisition, analysis, and interpretation of the data. SU, DT, CG, SY, NY, MS, and HZ were involved in the interpretation of the data and drafting the manuscript. CG, RS, and HZ supervised the entire study and were involved in study design, the acquisition, analysis, and interpretation of the data, and drafting the manuscript. The manuscript was reviewed and approved by all the authors.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.