- Open Access
Integrated analysis of competing endogenous RNA networks in peripheral blood mononuclear cells of systemic lupus erythematosus
Journal of Translational Medicine volume 19, Article number: 362 (2021)
Systemic lupus erythematosus (SLE) is an autoimmune disease with a complicated pathogenesis, and its aetiology has not been clearly unveiled. The lack of effective diagnosis and treatment methods makes it necessary to explore the molecular mechanism of SLE. We aimed to identify some critical signalling pathways and key competing endogenous RNAs (ceRNAs) underlying the molecular mechanism of SLE and to map out the systematic signalling networks by integrating the data on different kinds of RNAs.
Peripheral blood mononuclear cells (PBMCs) were collected from both SLE patients and healthy subjects, RNA was extracted from the PBMCs, and RNA libraries including ribosomal RNA-depleted strand-specific libraries and small RNA libraries were built for deep RNA sequencing (RNA-seq). RNA-seq yielded differential expression profiles of lncRNAs/circRNAs/miRNAs/mRNAs related to SLE. The DAVID database (v. 6.8) was employed for Gene Ontology (GO) and KEGG pathway analysis. ceRNA networks (circRNA/lncRNA-miRNA-mRNA) were constructed and visualized using Cytoscape software (v. 3.5.0). The TargetScan and miRanda databases were used to predict target relationships in ceRNA networks. qRT-PCR was used to verify our data.
Differential expression of ceRNAs related to SLE was detected in SLE patients’ PBMCs: 644 mRNAs (384 upregulated, 260 downregulated), 326 miRNAs (223 upregulated, 103 downregulated), 221 lncRNAs (79 upregulated, 142 downregulated), and 31 circRNAs (21 upregulated, 10 downregulated). We drew ceRNA signalling networks made up of the differentially expressed mRNAs/miRNAs/lncRNAs/circRNAs mentioned above, and the hub genes included IRF5, IFNAR2, TLR7, IRAK4, STAT1, STAT2, C2, and Tyk2. These hub genes were involved in ceRNA signalling pathways, such as the IL-17 signalling pathway and type I interferon signalling pathway.
We explored the differential expression profiles of various kinds of ceRNAs and integrated signalling networks constructed by ceRNAs. Our findings offer new insights into the pathogenesis of SLE and hint at therapeutic strategies.
Systemic lupus erythematosus (SLE) is a complex systemic autoimmune disease that is characterized by overproduction of autoantibodies against self-antigens, complement consumption, and chronic inflammation . According to statistics updated in September 2016, the prevalence of SLE varies between different areas in the world, and the prevalence could be as high as 241/100,000 people (95% CI 130,352) in specific areas, such as North America . In addition, SLE is a disease with a clear sex tendency that mainly affects women of childbearing age .
Multiorgan pathologies caused by immune complex deposition in SLE result in diverse clinical manifestations, such as arthritis and skin disease . Disease heterogeneity and diverse manifestation characteristics make early diagnosis and treatment difficult . The study of aetiology in SLE will improve our understanding of this disease and thus help with the early diagnosis and treatment of the disease.
Multiple factors, including endocrine factors, genetic factors, infections and environmental factors, contribute to the development of SLE, but the exact aetiology of this disease has not yet been clearly unveiled [6, 7].
Many studies have shown that various kinds of RNAs, including coding and noncoding RNAs, contain abnormal expression in SLE, and their distinct expression profiles are correlated with the distinct pathophysiological states in SLE patients. These coding and noncoding RNAs include messenger RNAs (mRNAs), long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs). According to earlier studies, lncRNAs and circRNAs regulate the expression of targeted genes (mRNAs) by competing for the binding of mRNAs to miRNAs; thus, all of these RNAs were defined as competing endogenous RNAs (ceRNAs). Different kinds of ceRNAs can be distinguished by their unique characteristics, such as their length. miRNAs are single-stranded RNAs of 18–22 nucleotides  and lncRNAs longer than 200 nucleotides , while circRNAs contain closed loops and are much more stable than other kinds of RNAs .
Signalling networks made of ceRNAs
With the development of bioinformatics and omics technologies, an increasing number of ceRNA networks have been identified, which will enrich our understanding of the pathogenesis of SLE. As summarized below, the ceRNA networks include lncRNA–miRNA–mRNA and circRNA–miRNA–mRNA.
For lncRNA–miRNA–mRNA interactions, Ye et al. found a ceRNA network consisting of 155 validated mRNAs, 15 miRNAs and 7 lincRNAs, and the network had a robust or weak “IFN signature”, which was closely related to the pathogenesis of SLE .
For circRNA–miRNA–mRNA networks, Tian et al. identified a ceRNA network consisting of mmu_circRNA_34428, miR-338-3p, miR-670-3p, miR-3066-5p, miR-210-5p and their corresponding mRNA targets in an NZB/W F1 lupus nephritis mouse disease model . Zhang et al. identified a ceRNA network made up of hsa_-circ_0012919, miR-125a-3p, CD70 and CD11a in SLE CD4+ T cells . Another study showed that circIBTK in a ceRNA network could inactivate protein kinase B (AKT) by binding to miR-29b in SLE, which impedes the development of SLE . Zhang et al. found that hsa_-circ_0049224 and hsa_-circ_0012919 affect the expression of DMNT1, CD70 and CD11a by sponging miR-125a-3p in SLE CD4+ T cells .
Studies also showed that different circRNAs play a different role in specific canonical signalling pathways in SLE. For example, hsa_circ_100226 targeted p65 by regulating miR-138 in the NF-κB inflammation pathway , while upregulated circRNA_002453 in the plasma of lupus nephritis (LN) was associated with the severity of renal damage .
From this evidence, researchers have recognized some pathways that are closely related to the pathogenesis of SLE, such as JAK-STAT (STAT1, STAT4) , oxidation-related superoxide dismutase (SOD), NADPH , apoptosis (Fas, Bcl2, TNF and IFN)  and complement and coagulation cascades.
ceRNA and ceRNA networks unveiled by former studies help us better understand the molecular pathogenesis of SLE. However, most ceRNA studies in SLE have focussed on a single ceRNA network, such as circRNA/lncRNA–miRNA–mRNAs. Discussion of both circRNA–miRNA–mRNA networks and lncRNA–miRNA–mRNA networks in the same study is rare, not mention to the crosstalk between the two different kinds of ceRNA networks.
To construct ceRNA regulatory networks of SLE, PBMCs were collected from both SLE patients and healthy subjects in this study. RNA was extracted from PBMCs, and then RNA libraries, including ribosomal RNA-depleted strand-specific libraries and small RNA libraries, were built for deep RNA sequencing (RNA-seq). The RNA-seq data between SLE patients and healthy subjects were compared to identify the differential expression profiles of RNAs (lncRNAs/circRNAs/miRNAs/mRNAs), and these RNAs were submitted to integrative analysis. The GO and KEGG pathway databases were utilized for functional and pathway analysis to identify some key RNAs and signalling pathways that play critical roles in the pathogenesis of SLE. Our findings will offer new, deeper insight into the molecular mechanisms of SLE pathogenesis and provide empirical support to therapeutic target investigations for SLE.
Blood samples collection
Blood samples were collected from 40 healthy subjects (H_B) and 40 SLE patients (SL_B) in Shenzhen People’s Hospital (Guangdong, China) from June 2018 to December 2019. All SLE patients and healthy subjects were of Han Chinese ancestry. Inclusion criteria: SLE patients diagnosed with SLE and healthy subjects without SLE or other immune diseases who were not treated with immunosuppressants. The SLE diagnostic criteria in this study were according to the ACR (the 1997 revised criteria of American College Rheumatology) . The SLE classification criteria and the characteristics of all patients are presented in Additional file 7: Table S5. The study was approved by the ethics committee of Shenzhen People's Hospital, and the experiments were carried out under the guidance of the Declaration of Helsinki of 1995. All participants signed informed consent forms.
Ethylene diamine tetraacetic acid (EDTA)-anticoagulated venous blood samples (5 ml for each donor) were obtained from SLE and healthy patients, and peripheral blood mononuclear cells (PBMCs) were prepared as described in a previous study . Briefly, blood samples were centrifuged at 3200 rpm for 10 min at room temperature. Plasma from blood samples was removed to an RNase-free tube for cryopreservation. PBMCs were obtained from remaining anticoagulated whole blood by using the Ficoll density gradient separation method (Sigma, USA). PBMCs were stored at −80 °C until further treatment.
Construction and sequencing of small-RNA library
Total RNA was extracted from PBMCs with TRIzol reagent (Invitrogen, CA, USA) or the Total RNA Purification Kit (LC Sciences, Houston, USA). The RNA amount and purity of each sample were quantified using a NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). RNA integrity was assessed by a Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0 (RIN > 8.4 in our experiment, data not shown) and was confirmed by electrophoresis with denaturing agarose gel. TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA) were used for small RNA library construction, and then the miRNA library was sequenced by an Illumina HiSeq 2500 with 50-base single-end reads following the vendor's recommended protocol. For more information, please refer to our Additional file 10: Methods-ceRNA analysis protocol).
Ribosomal RNA-depleted strand-specific library construction and sequencing
After total RNA extraction from PBMCs by TRIzol reagent, construction of two ribosomal RNA-depleted strand-specific libraries was conducted as previously reported [21, 22]. Briefly, the quantity and quality of RNA samples were analysed as described above. Then, RNA was purified by removing ribosomal RNA from total RNA using the Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, USA). After purification, the poly(A)- and poly(A) + RNA fractions were randomly fragmented into small pieces. Finally, the cDNA libraries were constructed by reverse transcription of the sheared products using the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). The average insert size for the final cDNA library was 300 ± 50 bp. Then, cDNA libraries were sequenced on an Illumina NovaSeq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) with 2 × 150 bp paired-end sequencing (PE150). To find more detailed information on the library construction, please refer to our Additional file 10: Methods-ceRNA analysis protocol).
Differential expression analysis and identification of RNAs
Different methods were utilized for expression analysis and identification of different kinds of RNAs. For miRNA, an in-house program of ACGT101-miR (LC Sciences, Houston, TX, US) or Cutadapt  was used to remove sequencing adaptors, low-quality bases and undetermined bases from the raw sequencing datasets. For other kinds of RNAs (mRNAs, circRNAs, lncRNAs), fastp was utilized (Additional file 10: Methods-ceRNA analysis protocol).
Then, the quality of expressed sequence tags was estimated using FastQC, and the valid reads were mapped to the human reference genome (GRCh38) by Bowtie 2  and TopHat2 . The differential expression levels of mRNAs/miRNAs/lncRNAs were calculated with StringTie  and Fisher’s exact test (F-test). circRNA expression levels from different samples or groups were calculated by SRPBM. Expression data with |log2 (fold change) |> 1 and p-value < 0.05 were defined as statistically significant. The value setting of |log2 (fold change) |> 1 for the threshold of fold change for RNAs can generate enough RNA candidates. Detailed information is provided in the Additional file 10: Methods-ceRNA analysis protocol.
Functional enrichment analysis
By volcano plot filtering and GO (http://www.geneontology.org) and KEGG analyses (http://www.genome.jp/kegg), the biological functions and pathways of the candidate target genes involved in ceRNA networks (lncRNAs/circRNAs-miRNAs-mRNAs) were predicted. Then, the ceRNA network-related target genes were put into the Database for Annotation, Visualization and Integrated Discovery (DAVID, v. 6.8) for GO and KEGG pathway analysis. Each GO and KEGG pathway term with a p value < 0.05 and |log2(fold change)|> 1 was defined as significant.
ceRNA network construction
TargetScan (v. 5.0) and miRanda (v. 3.3a)  were used to predict the potential targeting relationship of miRNA-mRNA (3′UTR) pairs and miRNA-lncRNA/circRNA pairs. UTR sequence data were acquired from Ensembl (http://ensemblgenomes.org). The regulatory networks of interest in our study included lncRNA-miRNA, miRNA-mRNA, circRNA-miRNA, lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA in SLE PBMCs. We set the context score cut-off as ≥ 50 in TargetScan and max energy to < −10 in miRanda, and a p value < 0.05 indicated significant binding possibilities between ceRNAs. Then, the overlapping predictions between the two programs were considered effective target pairs. All the ceRNA networks were drawn by Cytoscape (v. 3.5.0) .
Quantitative real-time polymerase chain reaction (qRT-PCR) was used to verify the identified sequencing data of circRNAs/lncRNAs/mRNAs/miRNAs from PBMCs of both H_B and SL_B.
To obtain miRNAs, total RNA was extracted from the PBMCs described above and subjected to first-strand cDNA synthesis using the TRUEscript 1st Strand cDNA SYNTHESIS Kit (Aidlab, Beijing, China). The specific primers used in the qRT-PCR are listed in Additional file 8: Table S6. qRT-PCR was performed using a SYBR PrimeScript miRNA RT-PCR Kit (TianGen Biotech, Beijing, China). Each qRT-PCR mix was referenced from a previous study, with minor modifications: 5 μl 2 × SYBR® Green Supermix, 0.5 μl specific forward primer, 0.5 μl specific reverse primer, 1 μl cDNA template and 3 μl ddH2O. The qRT-PCR programme was as follows: 95 °C for 3 min followed by 39 cycles of 95 °C for 10 s and 60 °C for 30 s, with a melt curve analysis. The experiment was performed in a real-time PCR system (analytikjena-qTOWER2.2, Germany). U6 was used as an internal control gene for miRNAs. Three independent biological replicates were run for each gene. The qRT-PCR data were analysed by the comparative 2-ΔΔCt method , and the average of three independent experiments for each gene was calculated as the relative expression level by GraphPad Prism (version 8.0.1).
To obtain lncRNAs and mRNAs, we used a PrimeScript™ RT Reagent Kit (Perfect Real Time, Takara, Japan) to prepare cDNA. The qRT-PCR was conducted as described above. GAPDH expression was used as the endogenous control for normalization of sample data.
To obtain circRNA , extracted total RNA was reverse-transcribed to cDNA with random primers using a RevertAid™ First-Strand cDNA Synthesis Kit (Fermentas, Waltham, MA). Real-time quantitative PCR was performed as described above. Divergent primers were designed for circRNAs (Additional file 8: Table S6). We used GAPDH as an internal control and calculated the relative gene expression with the 2-ΔΔCt comparative threshold cycle method.
The expression profile of mRNAs/miRNAs/lncRNAs/circRNAs in SLE
To explore the potential functions of mRNAs/miRNAs/lncRNAs/circRNAs in SLE, we analysed the expression profiles of these RNAs. From the transcriptional data, 644 transcripts were differentially expressed; 384 genes were upregulated, while 260 were downregulated (Fig. 1A, Additional file 3: Table S1). Further analysis showed that among all the differentially expressed transcripts, 99 were significant (60 genes upregulated, 39 downregulated) (Fig. 1B). The details of the differential expression of mRNAs are shown in a heatmap (Fig. 1C).
For miRNAs, the data are presented in a bar graph that summarizes the data for the differential expression profile of miRNAs in SL_B and H_B under the filter threshold (p < 0.01) (Fig. 1D, Additional file 4: Table S2). Detailed miRNA information is presented in the heatmap (Fig. 1F). When compared to H_B, SL_B contained 223 upregulated miRNAs and 103 downregulated miRNAs (p < 0.01). These data are also presented in volcano plots (Fig. 1E).
For lncRNAs, 221 differentially expressed lncRNAs were identified in our study. Among the differentially expressed lncRNAs, 79 lncRNAs were upregulated, 142 downregulated (Fig. 1G, Additional file 5: Table S3). The results are shown in volcano plots (Fig. 1H), and some top differentially expressed lncRNAs are shown in the heatmap (Fig. 1I).
circRNA differential expression between H_B and SL_B is shown in the volcano plots, bar chart and supplemental file (Fig. 2A–B, Additional file 6: Table S4), and some of the significant ones are presented in the heatmap (Fig. 2C). According to the data, there were 21 upregulated circRNAs and 10 downregulated circRNAs in SL_B. Compared to the H_B, there were more circRNAs in SL_B generated from exons of the genome (Fig. 2D). In contrast, more circRNAs in H_B (35.27%) than in SL_B (30.8%) came from introns.
GO enrichment and KEGG pathway enrichment analysis for the targeted genes in the ceRNA network
Gene Ontology (GO) was conducted for hub gene functional analysis (Fig. 3B). Hierarchical clustering of hub genes according to biological process (BP), molecular function (MF) and cellular component (CC) categories of the GO nomenclature, with screening criteria of p value < 0.05 and |log2 (fold changes) |≥ 1, was conducted. Twenty-five BP, 10 MF and 15 CC GO terms were found. BP terms included apoptotic process, immune system process, negative regulation of apoptotic process, cellular response to DNA damage stimulus, cell differentiation and others. CC included mitochondrion and other GO terms. MF contained ATP binding.
To reveal pathways involving the differentially expressed genes targeted by miRNAs in SLE, KEGG pathway analysis was performed (Fig. 3A). KEGG pathway analysis showed that some pathways participated in the development of SLE, such as cell apoptosis, inflammation pathways (Toll-like receptor signalling pathway, NF-kappa B signalling pathway and NOD-like receptor signalling pathway), and Th17 cell differentiation. Our results also indicated that the genes targeted by miRNAs function in SLE partially through apoptosis-, cell differentiation- and inflammation-related pathways.
Regulation networks of lncRNA-miRNA, circRNA-miRNA and miRNA-mRNA
When we examined the details of the ceRNA network, we found that the lncRNA/circRNA/miRNA/mRNA network could be divided into three subsets of interaction networks. The subsets of the network were distinguished by the correlation between two different types of RNAs, such as miRNA-mRNA, miRNA-lncRNA and miRNA-circRNA networks (Fig. 4A–C).
Regulatory network for lncRNAs/miRNAs/mRNAs and circRNAs/miRNAs/mRNAs
The interaction network summarized for lncRNA-miRNA-mRNA (Fig. 5A) showed that 6 of 13 of the top differentially expressed mRNAs were downregulated, and the 6 downregulated mRNAs correlated with 105 upregulated miRNAs and 141 downregulated lncRNAs in this regulated network (Fig. 5B).
At the same time, 7 mRNAs were upregulated, and these 7 upregulated mRNAs correlated with 57 downregulated miRNAs and 79 upregulated lncRNAs in this regulatory network (Fig. 5C).
The interaction network summarized for ceRNAs (circRNA-miRNA-mRNA) (Fig. 5D) contained 5 of the top differentially expressed and downregulated mRNAs, 18 upregulated miRNAs and 10 downregulated circRNAs in this regulated network (Fig. 5E).
At the same time, 7 mRNAs were upregulated, which correlated with 40 downregulated miRNAs and 14 upregulated circRNAs in this regulatory network (Fig. 5F).
ceRNA networks analysis by integrative method
The interaction network was summarized for ceRNAs (lncRNA/circRNA-miRNA-mRNA) (Fig. 6A–B). In the ceRNA network, there were 13 top differentially expressed mRNAs (DE-mRNAs), 31 differentially expressed circRNAs, 240 differentially expressed miRNAs and 221 differentially expressed lincRNAs.
A specific regulatory network including four types of ceRNAs (lncRNA/circRNA/miRNA/mRNA) was identified in our study. This regulatory network was made up of 1 circRNA, 5 mRNAs, 7 lncRNAs, and 8 miRNAs. Several hub genes identified in this ceRNA network have been reported previously, such as Tyk2 , TLR7 , and IRF5 .
Data validation by RT-PCR
The validation of RNA sequencing data was conducted by RT-PCR. A panel of 4 genes (mRNAs), 2 circRNAs, 2 lncRNAs and 2 miRNAs were selected from among the most differentially expressed genes. Our technical validation data (Additional file 1: Figure S2) showed that circRNA2013 was downregulated in SLE patient PBMCs, while in our RNA sequencing data, it was upregulated 5.35-fold (FC), which was contrary to our validated RT-PCR data. The expression level of circRNA3016 was lower in SLE patients than in control subjects in both sequencing and RT-PCR data; nevertheless, the data were significant only in the sequencing data and not in RT-PCR validation result. lncRNA-LIN02009 showed significantly higher levels in the sequencing data from SLE patients’ PBMCs, which was consistent with our RT-PCR data. Another lncRNA, SLC25A5-AS1, shared a similar expression trend in both sequencing and RT-PCR data: lower expression in SLE patient PBMCs than in controls, but the difference in the RT-PCR data was not significant. For microRNA RNA sequencing data validation, both sha-miR-125a_R + 2 and pc-3p-9143_61 were downregulated in SLE patients, but neither one was significant by RT-PCR. Among messenger RNAs, the C1QA gene was upregulated in PBMCs of SLE patients, while TYK2 was downregulated in both their RT-PCR and sequencing data. The expression level of STAT2 (ENST00000557235) presented a similar trend between the RT-PCR and sequencing data, while the expression of C2 (ENST00000299367) in our RT-PCR results was opposite to the sequencing result. In the sequencing data, C2 was 4.99-fold higher in SLE subjects than in healthy controls, while it presented a downregulated trend in our RT-PCR validation data.
SLE is an autoimmune disease with complicated pathogenesis and aetiology. The existing treatment methods, including immunosuppressive agents and new therapeutic interventions, have increased the survival rate of SLE patients. Nevertheless, both the main molecular mechanisms of immune dysfunction and the key driver of pathogenesis in SLE remain unclear. Therefore, there are still some challenges in both the early diagnosis and individualized treatment for this disease .
Abnormal epigenetic regulation has been reported previously and potentially plays an important role in the pathogenesis of SLE, and ceRNAs are important forms of epigenetic modifications . In our study, we used RNA-seq to detect differential expression levels of lncRNAs, circRNAs, miRNAs and mRNAs in both SLE patients and healthy individuals. We discovered some hub genes and disease-related signalling pathways, such as IRF5, IFNAR2, TLR7, IRAK4, STAT1, STAT2, C2, and Tyk2 and their pathways, and signalling networks, such as circRNA2692-miRNA125b-stat2, circRNA2692-miRNA125b-C2, circRNA2013-PC-3P-9143_61-IFNAR2, and circRNA2013-hsa-mir-6741-nfkb1.
By KEGG and GO enrichment analysis and integrated methods, some regulatory pathways containing hub genes and related ceRNA signalling networks are presented in Table 1.
ceRNAs identified in previous studies
lncRNAs such as linc0949 and linc0597 , GAS5, lnc7074, lnc0640, lnc5150, NEAT-1, MALAT1 and lnc-DC [36,37,38] have a potential role in the pathogenesis of SLE. Among them, GAS5, NEAT-1, and MALAT-1 were verified by our study (Additional file 5: Table S3).
Some circRNAs identified in previous studies have potential functions in the pathogenesis of SLE, such as hsa_circ_0045272 , hsa_circ_0049224, hsa_circ_0049220 , hsa_circRNA_001308, and hsa_circRNA_407176 . In our study, we also unveiled several new circRNAs in SLE, such as circRNA 2092, circRNA 2013, and circRNA 3016 (Additional file 6: Table S4). Some new identical circRNAs in this study have not been reported in SLE studies, but they were demonstrated to function in other biological processes. For example, circRNA 2092 plays a significant role in the development of secondary hair follicles .
Some miRNAs related to SLE were also recognized in both previous studies and our current study (Additional file 4: Table S2): miR-126, miR-21, miR-146a, miR-451a, miR-148a, miR-873, miR-1246, miR-25, miR-302, miR-93, miR-9b, miR-22, miR-151a-3p and miR-766-3p .
mRNAs acknowledged by previous reports to be associated with the pathogenesis of SLE include stat4, Itgax, Lyn, Tnfaip3 Cdkn1b, Intergenic, Ptpn22, Trex1, Hla-a, Ifi44, Ifi44l, Ifit1, Ly6e and Irf7 . Interestingly, these genes were well represented in the current study (Additional file 3: Table S1).
The pathways involved in SLE
By using the KEGG and GO databases, 8562 KEGG pathways and GO terms were discovered in this study. We also summarize 16 frequently mentioned KEGG pathways and GO terms from previous studies (Fig. 7, Table 2). Some of the KEGG pathways and GO terms identified in both previous studies and our current study are shown in Table 1. At the same time, common elements of GO terms and KEGG pathways between the earlier studies and this study are also summarized (Table 1, Additional file 1: Figure S1).
As discussed by previous studies, there are many pathways that play critical roles in the pathogenesis of SLE, such as IL-17 signalling , the Jak-STAT pathway , type I interferon signalling , apoptosis , the NF-kappa B pathway , the TNF signalling pathway , and the Toll-like receptor signalling pathway .
The KEGG pathways identified by our current study included cell apoptosis, inflammation pathways (Toll-like receptor signalling pathway, NF-kappa B signalling pathway and NOD-like receptor signalling pathway), Th17 cell differentiation, apoptosis, Th1 and Th2 cell differentiation, TNF signalling pathway, Jak-STAT signalling pathway, natural killer cell mediated cytotoxicity, antigen processing and presentation, T cell receptor signalling pathway, B cell receptor signalling pathway, etc.
Some SLE-related GO terms identified by here include apoptotic process, immune system process, negative regulation of apoptotic process, DNA damage stimulus, cell differentiation, innate immune response, T cell homeostasis, regulation of B cell differentiation, etc.
Summary of overlapping pathways between earlier studies and ours
ceRNA networks modify or control the expression of key genes in many diseases. For instance, ceRNA network analysis has been mostly used in cancer research . With ceRNA studies extending to autoimmune diseases, many important molecules and key signalling networks have been unveiled, such as the Jak-STAT signalling pathway , complement and coagulation cascades , and apoptosis. The pathways above are summarized in Additional file 9: Table S7.
Several significant ceRNA-related pathways were unveiled by our current study, such as osteoclast differentiation, legionellosis, antigen processing and presentation, natural killer cell-mediated cytotoxicity, the IL-17 signalling pathway, Th17 cell differentiation, and acute myeloid leukaemia.
Some common pathways identified by our study and earlier ones were identified. For example, the IL-17 signalling pathway  and type I interferon signalling pathway  were identified in both previous reports and our study, which strongly suggests that these two pathways might play key roles in the pathogenesis of SLE. On the other hand, some pathways recognized by our study have not been well studied before, and these new pathways might provide new insight for aetiology studies of SLE. The ceRNAs identified by our study also add candidates for potential diagnostic biomarker screening.
As more research has been conducted on ceRNAs, a number of genes (molecules) and regulatory pathways that potentially play important roles in the development of SLE have been found. Earlier studies in ceRNA networks were mostly focussed on single regulatory pathways, such as circRNA–miRNA–mRNAs or lncRNA–miRNA–mRNAs. There were no reports discussing the possible relationship between circRNA–miRNA–mRNA sets and lncRNA–miRNA–mRNA sets. Our study conducted an integrated analysis of the ceRNA network for lncRNAs/circRNAs–miRNAs–mRNAs, and some important ceRNA networks were identified. For example, miRNAs (PC-3p-9143_61, hsa-miR-548ay-3p,) lncRNAs (LINC02009) and circRNA 2013 were found in the same regulatory network and correlated with the hub gene TYK2.
Therefore, our ceRNA network study of PBMCs from SLE patients recognized some new ceRNAs and signalling pathways, which will promote the identification of potential diagnostic biomarkers and ceRNA-related regulatory pathways in SLE. Our findings will help us better understand the pathogenesis of SLE, which will aid in research on therapeutic targets.
Limitations of this article
First, methodological flaws: the sample size of the SLE subjects and healthy controls was relatively small. Larger samples are needed to improve the statistical validity of future studies. Another potential shortcoming is the lack of detailed patient selection criteria and descriptions of their clinical characteristics. SLE is a complex autoimmune disease with a highly varied clinical presentation, and the pathogenic mechanisms involved in SLE are varied. Therefore, more precise criteria, such as activity markers (C3, C4, dsDNA), should be adapted in future studies rather than merely pooling all the patients together for analysis. Finally, the validation methods should be expanded. The hub ceRNAs (genes) were identified by RNA-seq and verified by RT-PCR. Various kinds of biotechnologies for testing PBMCs are needed for further verification of biomarkers, such as ELISA and western blotting.
In conclusion, this integrated analysis of ceRNA networks identified several hub ceRNAs (genes) and signalling pathways associated with the pathogenesis of SLE. By combining some important signalling pathways and genes identified by previous studies and our study, we have identified some hub genes, KEGG pathways and GO terms that may play critical roles in SLE development. The hub genes and ceRNAs recognized in our study included mRNAs (IRF5, IFNAR2, TLR7, STAT1, C2, Tyk2) and some circRNAs, miRNAs and lncRNAs. The pathways identified here included the Jak-STAT signalling pathway, apoptosis, NF-kappa B signalling pathway, TNF signalling pathway, Toll-like receptor signalling pathway, IL-17 signalling pathway, and type I interferon signalling pathway. Notably, the IL-17 signalling pathway and type I interferon signalling pathway were two of the most significant pathways shown in our study (Fig. 7).
Availability of data and materials
The raw data for small RNA sequencing and ribosomal RNA-depleted sequencing in our study are available from [Gene Expression Omnibus: GSE146410], but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available.
Systemic Lupus Erythematous
Peripheral blood mononuclear cells
Long non-coding RNAs
Quantitative real-time polymerase chain reaction
Molecular functions (MF)
Kyoto encyclopaedia of genes and genomes
Database for Annotation, Visualization, and Integrated Discovery
Activate protein kinase B
DNA methyltransferase 1
Nicotinamide adenine dinucleotide phosphate
Systemic Lupus Erythematosus Disease Activity Index
Erythrocyte sedimentation rate
Ethylene Diamine Tetraacetic Acid
Singh RP, Waldron RT, Hahn BH. Genes, tolerance and systemic autoimmunity. Autoimmun Rev. 2012;11(9):664–9.
Rees F, Doherty M, Grainge MJ, Lanyon P, Zhang W. The worldwide incidence and prevalence of systemic lupus erythematosus: a systematic review of epidemiological studies. Rheumatology (Oxford). 2017;56(11):1945–61.
Zhou Z, Sun B, Huang S, Zhao L. Roles of circular RNAs in immune regulation and autoimmune diseases. Cell Death Dis. 2019;10(7):503.
Fava A, Petri M. Systemic lupus erythematosus: diagnosis and clinical management. J Autoimmun. 2019;96:1–13.
Smith PP, Gordon C. Systemic lupus erythematosus: clinical presentations. Autoimmun Rev. 2011;10(1):43–5.
Tsokos GC. Systemic lupus erythematosus. N Engl J Med. 2011;365(22):2110–21.
Zou YF, Feng CC, Zhu JM, Tao JH, Chen GM, Ye QL, et al. Prevalence of systemic lupus erythematosus and risk factors in rural areas of Anhui Province. Rheumatol Int. 2014;34(3):347–56.
Eulalio A, Huntzinger E, Izaurralde E. Getting to the root of miRNA-mediated gene silencing. Cell. 2008;132(1):9–14.
Wu GC, Li J, Leng RX, Li XP, Li XM, Wang DG, et al. Identification of long non-coding RNAs GAS5, linc0597 and lnc-DC in plasma as novel biomarkers for systemic lupus erythematosus. Oncotarget. 2017;8(14):23650–63.
Zhang M-Y, Wang J-B, Zhu Z-W, Li L-J, Liu R-S, Yang X-K, et al. Differentially expressed circular RNAs in systemic lupus erythematosus and their clinical significance. Biomed Pharmacother. 2018;107:1720–7.
Ye H, Wang X, Wang L, Chu X, Hu X, Sun L, et al. Full high-throughput sequencing analysis of differences in expression profiles of long noncoding RNAs and their mechanisms of action in systemic lupus erythematosus. Arthritis Res Ther. 2019;21(1):1–17.
Tian S, Liu X, Fan Q, Ma J, Yao L, Li Y. Microarray expression and functional analysis of circular RNAs in the glomeruli of NZB/W F1 mice with lupus nephritis. Exp Ther Med. 2019;18(4):2813–24.
Zhang C, Wang X, Chen Y, Wu Z, Zhang C, Shi W. The down-regulation of hsa_circ_0012919, the sponge for miR-125a-3p, contributes to DNA methylation of CD11a and CD70 in CD4(+) T cells of systemic lupus erythematous. Clin Sci. 2018;132(21):2285–98.
Xin W, Zhang C, Wu Z, Yue C, Shi W. CircIBTK inhibits DNA demethylation and activation of AKT signaling pathway via miR-29b in peripheral blood mononuclear cells in systemic lupus erythematosus. Arthritis Res Ther. 2018;20(1):1–10.
Wei ZJ, Liu J, Qin J. miR-138 suppressed the progression of osteoarthritis mainly through targeting p65. Eur Rev Med Pharmacol Sci. 2017;21(9):2177–84.
Qingqing O, Qin H, Zhenlan J, Jinjun Z, Guo-Ping S, Min Y. Using plasma circRNA_002453 as a novel biomarker in the diagnosis of lupus nephritis. Mol Immunol. 2018;101:531–8.
Fang TJ, Lin YZ, Liu CC, Lin CH, Yen JH. Methylation and gene expression of histone deacetylases 6 in systemic lupus erythematosus. Int J Rheum Dis. 2016;19(10):968–73.
Shah D, Mahajan N, Sah S, Nath SK, Paudyal B. Oxidative stress and its biomarkers in systemic lupus erythematosus. J Biomed Sci. 2014;21(1):1–13.
Munoz L, van Bavel C, Franz S, Berden J, Herrmann M, Vlag DVJ. Apoptosis in the pathogenesis of systemic lupus erythematosus. Lupus. 2008;17(5):371–5.
Hochberg M. Updating the American College of Rheumatology revised criteria for the classification of systemic lupus erythematosus. Arthritis Rheum. 1997;40(9):1725.
Yuan Y, Jiaoming L, Xiang W, Yanhui L, Shu J, Maling G, et al. Analyzing the interactions of mRNAs, miRNAs, lncRNAs and circRNAs to predict competing endogenous RNA networks in glioblastoma. J Neuro Oncol. 2018;137(3):493–502.
Liu S, Xie X, Lei H, Zou B, Xie L. Identification of key circRNAs/lncRNAs/miRNAs/mRNAs and pathways in preeclampsia using bioinformatics analysis. Med Sci Monitor. 2019;25:1679–93.
Kechin A, Boyarskikh U, Kel A, Filipenko M. cutPrimers: a new tool for accurate cutting of primers from reads of targeted next generation sequencing. J Comput Biol. 2017;24(11):1138–43.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):1–13.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Manyam G, Ivan C, Calin GA, Coombes KR. targetHub: a programmable interface for miRNA-gene interactions. Bioinformatics (Oxford, England). 2013;29(20):2657–8.
Liu S, Wang D, Liu Y. Extracellular RNA in systemic lupus erythematosus. ExRNA. 2019;1(1):1–7.
Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3(6):1101–8.
Cunninghame Graham DS, Morris DL, Bhangale TR, Criswell LA, Syvänen AC, Rönnblom L, et al. Association of NCF2, IKZF1, IRF8, IFIH1, and TYK2 with systemic lupus erythematosus. PLoS Genetics. 2011;7(10):e1002341.
Santiago-Raber ML, Dunand-Sauthier I, Wu T, Li QZ, Uematsu S, Akira S, et al. Critical role of TLR7 in the acceleration of systemic lupus erythematosus in TLR9-deficient mice. J Autoimmun. 2010;34(4):339–48.
Tada Y, Kondo S, Aoki S, Koarada S, Inoue H, Suematsu R, et al. Interferon regulatory factor 5 is critical for the development of lupus in MRL/lpr mice. Arthritis Rheumatol. 2011;63(3):738–48.
Lerang K, Gilboe IM, Steinar Thelle D, Gran JT. Mortality and years of potential life loss in systemic lupus erythematosus: a population-based cohort study. Lupus. 2014;23(14):1546–52.
Gérard C, Gonze D, Lemaigre F, Novák B. A model for the epigenetic switch linking inflammation to cell transformation: deterministic and stochastic approaches. PLoS Comput Biol. 2014;10(1):e1003455.
Wu Y, Zhang F, Ma J, Zhang X, Wu L, Qu B, et al. Association of large intergenic noncoding RNA expression with disease activity and organ damage in systemic lupus erythematosus. Arthritis Res Ther. 2015;17(1):1–11.
Yang H, Liang N, Wang M, Fei Y, Sun J, Li Z, et al. Long noncoding RNA MALAT-1 is a novel inflammatory regulator in human systemic lupus erythematosus. Oncotarget. 2017;8(44):77400–6.
Wu GC, Hu Y, Guan SY, Ye DQ, Pan HF. Differential plasma expression profiles of long non-coding RNAs reveal potential biomarkers for systemic lupus erythematosus. Biomolecules. 2019;9(6):206.
Zhao CN, Mao YM, Liu LN, Li XM, Wang DG, Pan HF. Emerging role of lncRNAs in systemic lupus erythematosus. Biomed Pharmacother. 2018;106:584–92.
Li LJ, Zhu ZW, Zhao W, Tao SS, Li BZ, Xu SZ, et al. Circular RNA expression profile and potential function of hsa_circ_0045272 in systemic lupus erythematosus. Immunology. 2018;155(1):137–149.
Zhang C, Huang J, Chen Y, Shi W. Low expression and clinical value of hsa_circ_0049224 and has_circ_0049220 in systemic lupus erythematous patients. Med Sci Monitor. 2018;24:1930–5.
Yin R, Wang Y, Wang Z, Zhu Y, Bai W. Discovery and molecular analysis of conserved circRNAs from cashmere goat reveal their integrated regulatory network and potential roles in secondary hair follicle. Electronic J Biotechnol. 2019;41:37–47.
Husakova M. MicroRNAs in the key events of systemic lupus erythematosus pathogenesis. Biomed Papers Med Faculty Univ Palacky Olomouc Czechoslovakia. 2016;160(3):327–42.
Song W, Tang D, Chen D, Zheng F, Huang S, Xu Y, et al. Advances in applying of multi-omics approaches in the research of systemic lupus erythematosus. Int Rev Immunol. 2020;39(4):163–73.
Chen S-Y, Liu M-F, Kuo P-Y, Wang C-R. Upregulated expression of STAT3/IL-17 in patients with systemic lupus erythematosus. Clinical Rheumatology. 2019;38(5):1361–1366.
Salmon M, Gordon C. The role of apoptosis in systemic lupus erythematosus. Rheumatology. 1999;38(12):1177–1183.
Brightbill HD, Suto E, Blaquiere N, Ramamoorthi N, Sujatha-Bhaskar S, Gogol EB, et al. NF-κB inducing kinase is a therapeutic target for systemic lupus erythematosus. Nat Commun. 2018;9(1):1–14.
Zhu LJ, Landoltmarticorena C, Li T, Yang X, Yu XQ, Gladman DD, et al. Altered expression of TNF-alpha signaling pathway proteins in systemic lupus erythematosus. J Rheumatol. 2010;37(8):1658–66.
Celhar T, Fairhurst AM. Toll-like receptors in systemic lupus erythematosus: potential for personalized treatment. Front Pharmacol. 2014;5:265.
Li Y, Ma B, Yin Z, Liu P, Liu J, Li J. Competing endogenous RNA network and prognostic nomograms for hepatocellular carcinoma patients who underwent R0 resection. J Cell Physiol. 2019;234(11):20342–53.
Liang Y, Xie SB, Wu CH, Hu Y, Ye D-Q. Coagulation cascade and complement system in systemic lupus erythematosus. Oncotarget. 2017;9(19):14862–81.
Mary K, Crow, Mikhail, Olferiev, Kyriakos A, Kirou. Type I interferons in autoimmune disease. Ann Rev Pathol. 2019;14:369–393.
This work was supported by the National Natural Science Foundation for Young Scientists of China (Grant No32000641), the Key Research and Development Program of Guangdong Province (No. 2019B020229001), the Science and Technology Plan of Shenzhen (No. JCYJ20200109144218597 and JCYJ20190807145815129), Sanming Project of Medicine in Shenzhen (No. SYJY201704 and No. SYJY201705), Shenzhen Key Medical Discipline Construction Fund (No. SZXK011).
Ethics approval and consent to participate
The study was approved by the ethics committee of Shenzhen People's Hospital (the committee’s reference number: LL-KY-2019589), and the experiments were undertaken and carried out by the guidelines of the Declaration of Helsinki of 1995 (as revised in Edinburgh 2000). Besides, all participants signed informed consent forms.
Consent for publication
All the authors have read and approved the paper and declare no potential conflicts of interest in the paper. All the authors agree to publish this paper.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Overlapping common signalling pathways are shown in a Venn diagram. Overlapping of KEGG pathways and GO terms was analysed by comparing both the KEGG and GO terms from previous studies, our current study and some of the significant terms in our current study.
Validation of RNA sequencing data by RT-PCR. Four genes, 2 circRNAs, 2 lncRNAs and 2 miRNAs were analysed by qRT-PCR, and their relative expression levels were normalized to the housekeeping gene GAPDH or U6. Data are represented as mean ± SEM. *p < 0.05 and **p < 0 01 indicate that the differential expression of genes was significant, while ns indicates no significance (p ≥ 0.05).
SL_B vs H_B_transcripts-mRNAs
SL_B vs H_B differential expression of miRNAs
SL_B vs H_B_lncRNA_differential_expression
SL_B vs H_B_circRNA_differential_expression
Clinical data of participants.
The primers using for PCR amplification.
Common elements of GO terms and KEGG pathways
Protocols for RNA extraction, library construction, sequencing and Bioinformatics analysis
About this article
Cite this article
Song, W., Qiu, J., Yin, L. et al. Integrated analysis of competing endogenous RNA networks in peripheral blood mononuclear cells of systemic lupus erythematosus. J Transl Med 19, 362 (2021). https://doi.org/10.1186/s12967-021-03033-8
- Integrative analysis
- RNA sequencing