Systemic lupus erythematosus (SLE) is an autoimmune disease, characterized by the
production of autoantibodies and the involvement of multi-systems. In order to explore
its molecular mechanism, we, using bioinformatics analysis and wet lab experiments,
identified two key genes in SLE patients.
Briefly, we analyzed differentially expressed genes (DEGs) and differentially methylated
genes (DMGs) from the GEO database (GSE50772 and GSE59250). The key genes were screened
by constructing a weighted gene co-expression network analysis (WGCNA) network. We
further performed immune infiltration analysis, functional enrichment analysis, and
motif transcriptional regulation analysis. Wet lab verification, including quantitative
real-time polymerase chain reaction (qPCR) and western blotting, was performed to
further confirm the expression of the identified key genes in SLE patients. The overall
workflow is shown in Figure S1A and the detailed research methods are in Supplementary
Based on the screening criteria (adjusted P < 0.05), 6,817 DEGs in GSE50772 were screened
(Fig. 1A). The top 10 upregulated or downregulated genes were clustered by visual
inspection of the heatmap (Fig. 1B). GO analysis indicated the enriched pathways of
these DEGs in each category (Fig. S1B). KEGG pathway enrichment showed that they were
mainly involved in herpes simplex virus 1 infection (Fig. S1B). Using the R package
“ChAMP”, 238 differential methylated probes (DMPs) were screened out from GSE59250
(Fig. 1C). By comparing the DMGs and the DEGs, we found four hypermethylated downregulated
genes and 35 hypomethylated upregulated genes (Fig. S1C). To determine which specific
genes are highly associated with SLE, we used the expression profile of GSE50772 to
construct a WGCNA network and explore the disease modules related to SLE. Samples
were first clustered to detect outliners and none was excluded (Fig. S1D). Soft threshold
power β = 12 (scale-free R
2 = 0.9) was selected (Fig. S1E). In total, 19,746 genes were classified into 18 distinct
modules based on the TOM matrix (Fig. 1D). Through further analysis between modules
and traits, we found that the light-green module had the highest correlation with
SLE phenotype (cor = −0.69, P = 1e−12; Fig. 1E). The correlation between this module
and the gene expression profiles was calculated. The gene significance score and module
membership in the light-green module was also highly correlated (cor = 0.6, P < 1e−200;
Fig. 1F). Therefore, this module was used for a further selection of key driver genes.
By intersecting the WGCNA light-green module genes with DEGs-DMGs (Fig. S1F), we obtained
two key genes OSBPL3 and TNRC6C (Fig. 1G). By analyzing the correlation between core
genes and immune infiltration in the dataset, we further explored their potential
molecular mechanisms in disease progression with regard to the immune component. We
quantified the immune infiltration level of each patient using the CIBERSORT algorithm
and found that both OSBPL3 and TNRC6C were significantly positively correlated with
CD4 memory T cells resting and negatively correlated with neutrophils and monocytes
(Fig. 1H). OSBPL3 was also positively correlated with γδT cells and NK cells resting,
while negatively correlated with activated dendritic cells (Fig. 1H). TNRC6C was positively
correlated with naïve CD4 T cells and naïve B cells, while negatively correlated with
memory B cells and plasma cells (Fig. 1H). The specific signal pathways involving
the core genes in the disease progression were investigated too. GSVA data demonstrated
that upregulated OSBPL3 was mainly enriched in G2M checkpoint and UV response DNA
pathways. Upregulated TNRC6C was mainly enriched in the G2M checkpoint and bile acid
metabolism pathways. Downregulated OSBPL3 and TNRC6C were both mainly enriched in
TNFα signaling via NF-κB and interferon- α/γ response (Fig. 1I). We also analyzed
GO and KEGG enrichment by GSEA (Fig. 1J, K, respectively). GO enrichment indicated
that the tumor necrosis factor superfamily (TNFSF) cytokine production pathway was
enriched in low expression OSBPL3, while regulation of myeloid leukocyte mediated
immunity pathway was enriched in low expression TNRC6C (Fig. 1J). In the transcription
factors (TFs) analysis, we found these core genes were regulated by multiple TFs sharing
common mechanisms. These TFs were enriched by cumulative recovery curve, motif-TF
annotation, and selection of key genes (Fig. 1L, M). The motif with the highest normalized
enrichment score (NES: 8.79) was cisbp_M2425, and OSBPL3 was enriched in it (Fig. 1L).
The whole set of significant motifs (NES > 3) and corresponding TFs that were enriched
for either one or both of the two core genes are listed in Supplementary File 2. Lastly,
we performed qPCR to quantify the relative transcription levels of hub genes in SLE
and healthy controls (HCs). As shown in Figure 1N, OSBPL3 and TNRC6C transcription
levels were significantly decreased in SLE patients compared with HCs. To further
validate protein translation levels, we performed the western blotting of the PBMC
from SLE patients and HCs (each n = 6). We observed that the expression levels of
OSBPL3 and TNRC6C were greatly reduced in SLE patients (Fig. 1O).
Analysis of the study results. (
) Analysis of dataset GSE50772. (A) The volcano plot illustrates the DEGs between
control and SLE patients, including 3250 upregulated and 3567 downregulated genes.
(B) The heatmap clustering the expression of the top 10 upregulated or downregulated
DEGs among each sample. (
) The heatmap showing the DMPs between control and SLE patients, including 150 hypermethylated
probes (61 genes) and 88 hypomethylated probes (84 genes). (D
F) The WGCNA network analysis of 81 samples in dataset GSE50772. (D) Dendrogram of
co-expression modules (n = 18) identified by WGCNA, including black (n = 899), blue
(n = 1,858), cyan (n = 245), dark green (n = 177), dark orange (n = 126), dark red
(n = 307), green-yellow (n = 285), grey (n = 8,787), grey 60 (n = 477), light cyan
(n = 237), light green (n = 2,461), light yellow (n = 219), magenta (n = 564), orange
(n = 353), pink (n = 501), red (n = 702), tan (n = 280), and turquoise (n = 1,268)
modules. Different color branches of the cluster tree represent different modules,
and the color bands below the dendrogram show the cluster membership according to
different clustering methods. (E) Module–trait correlation heatmap. The color scale
on the right shows correlations from −1 (red) to 1 (blue). (F) A scatter plot showing
the correlation analysis between the gene significance score and module membership
in the light green module. (G) Two key genes identified by intersecting the light
green module with DEGs-hypermethylated genes. Relative expression of them in SLE and
control samples in the dataset was shown. (H–K) Functional analysis of the core genes.
(H) The correlation between OSBPL3 and TNRC6C and immune cell infiltration. (I) The
pathway enrichment analysis in high- or low-expression of OSBPL3 and TNRC6C in SLE
patients. (J) The enrichment of GO in high- or low-expression of OSBPL3 and TNRC6C
by GSEA analysis. (K) The enrichment of KEGG pathways in high- or low-expression of
OSBPL3 and TNRC6C by GSEA analysis. (L, M) Motif enrichment analysis of regulatory
regions in OSBPL3 and TNRC6C in SLE patients. (L) Analysis of global mean and SD estimation
as well as the recovery curves of motifs with high NES. The selection of significant
motifs is done based on NES, which is calculated for each motif based on the AUC distribution
of all the motifs for the gene set [(x-mean)/SD]. The red line is the mean of the
recovery curves of each motif. The green line is the mean value with standard deviation.
The blue line is the recovery curve of the indicated motif. (M) AUC histogram of hub
genes. (N, O) Relative transcription and translation levels of two hub genes. (N)
Relative transcription levels of OSBPL3 and TNRC6C in SLE patients (n = 22) or healthy
controls (n = 23) (both P = 0.0118). (O) Protein levels of OSBPL3 and TNRC6C in SLE
patients (n = 6) and healthy controls (n = 6). AUC, area under the curve; DC, dendritic
cell; GO, gene ontology; GSVA, gene set variation analysis; H, healthy control; HExp,
high expression; KEGG, Kyoto Encyclopedia of Genes and Genomes; LExp, low expression;
NK, natural killer; NS, no significant; NSE, normalized enrichment score; S, SLE patient;
SD, standard deviation; SLE, systemic lupus erythematosus. ∗P < 0.05. ∗∗P < 0.01.
To date, distinctive transcriptional signatures have been identified in SLE, showing
the importance of how key genes play a role in its pathogenesis. OSBPL3 belongs to
the oxysterol-binding protein family and functions in cell adhesion, the actin cytoskeleton,
cellular lipid metabolism, vesicle transport, and cell signaling.
TNRC6C belongs to the trinucleotide repeat containing the 6/GW182 family and plays
important roles in RNA-mediated gene silencing by micro-RNAs.
In our analysis, they were downregulated in SLE patients, possibly due to hypermethylation
in their promoter region, which is consistent with the regulatory mechanism reported
in other diseases.
We also showed that downregulated OSBPL3 was associated with interleukin (IL)-8,
and TNFSF cytokine production and signaling, while downregulated TNRC6C was associated
with TNFα signaling via NF-κB. IL-8 is an important proinflammatory factor that induces
neutrophil activation, chemotaxis, and neutrophil extracellular trap (NET) formation.
NETs are one of the sources of self-antigens that destroy immune tolerance.
TNFSF ligand-receptors signaling pathways have been demonstrated to be active in
inflammatory and autoimmune disease.
Taken together, our findings imply that decreased OSBPL3 and TNRC6C might be involved
in the pro-inflammatory process in SLE. By analyzing immune cell infiltration, we
found both OSBPL3 and TNRC6C were inversely correlated with monocytes and neutrophils.
OSBPL3 was also inversely associated with activated dendritic cells, while TNRC6C
was inversely associated with plasma cells and memory B cells. Dendritic cells are
crucial in regulating peripheral tolerance to self-antigens in SLE and could induce
B cell differentiation, which ultimately leads to the production of pathogenic autoantibodies.
The interaction between the two hub genes and immune cells in the disease pathogenesis
should be highlighted in further research. In the TFs analysis, we identified TF regulatory
network associated with OSBPL3 and TNRC6C from publicly available SLE datasets. We
found that cisbp_M2425 was the most significantly enriched motif in OSBPL3, suggesting
a potential binding site of TFs.
In conclusion, based on bioinformatic data analysis, two novel hypermethylated, lowly
expressed hub genes (OSBPL3 and TNRC6C) were first identified and confirmed in SLE
patients. They may be involved in the pathophysiologic mechanism of SLE through aberrant
cell adhesion, induction of pro-inflammatory cytokines, and interaction with imbalanced
immune cells. Our study provides a new methodology to study the molecular mechanisms
of SLE from both epigenome and transcriptome. Further studies are warranted to validate
the hypothesis as well as the predicted TFs by regulatory network analysis.
This study was approved by the Ethical Committee of Sichuan Provincial People's Hospital.
All subjects signed informed consent.
Conflict of interests
The authors declare that they have no competing interests.
This research was supported by the
National Natural Science Foundation of China
(No. 81802504, 81872207), Sichuan Science and Technology Bureau, China (No. 2019YFS0439,
2020JDJQ0067, 2020JDRC0118, 2021YJ0564, and 2022YFH0005), the Science and Technology
Innovation Project of Chengdu, China (No. 2021-YF05-00225-SN), and a Sichuan Medical
Association grant (China) (No. Q19037).
Consent for publication
All authors signed a consent for publication form.
Availability of data and materials
The datasets analyzed during the current study are available in the GEO repository,