Transcriptional Landscape of Enhancer RNAs in Peripheral Blood Mononuclear Cells from Patients with Systemic Lupus Erythematosus

Introduction

Systemic lupus erythematosus (SLE) is a multisystem autoimmune disease characterized by loss of immune tolerance, immune cell activation, and pathogenic autoantibody production.¹ The pathology of SLE, which is a highly clinically heterogeneous disease, remains largely unclear, but involves the dysregulation of multiple immune cell lineages and genes. Epigenetic alterations in peripheral blood immune cells, including T cells and B cells, have been implicated in the pathogenesis of SLE.² Epigenomic disease-related gene regulation in peripheral blood mononuclear cells (PBMCs) has been found to contribute to abnormal inflammatory cytokine production during the course of the disease.³ In particular, non-coding (nc)RNAs are involved in immune cell responses and their potential role in SLE has gained increased attention.⁴ We previously reported that microRNAs,⁵ circular RNAs,⁶ and long ncRNAs⁷ in PBMCs play significant roles in SLE pathogenesis and serve as potential biomarkers.

Recently, it has been reported that enhancers and super-enhancers (SEs) participate in transcriptional control and play roles in maintaining cell identity and homeostasis in all cell types, and may be involved in increased autoimmune disease risk and disease development through spatiotemporal regulation of gene transcription.^8–10 Enhancer (e)RNAs, which are a type of ncRNA, are transcribed from enhancers, which are cis-regulatory elements in the genome.¹¹ Some studies have examined eRNAs from monocytes of SLE patients, and discussed the roles of SEs and eRNAs in the pathogenesis of autoimmune diseases, including SLE.^7,12 Intriguingly, eRNAs have been identified as one of the most plausible candidates in the regulation of functional interactions between enhancers and promoters.¹³ Growing evidence suggests that at least some eRNAs play critical roles in enhancer functions by recruiting transcription factors (TFs) to specific enhancers, promoting the establishment of enhancer–promoter looping.¹³ Based on the intensity of TF binding, a distinct subset of enhancers have been termed SEs,¹⁴ which are more powerful in regulating gene expression than regular enhancers. Various studies have demonstrated the critical roles of eRNAs in target gene transcription.^10,15,16 A recent study showed that eRNAs exert key roles in cancer pathogenesis and are potential therapeutic targets worthy of further investigation.¹¹ However, the relationships between eRNAs, associated SEs, and TFs and SLE pathogenesis in PBMCs have not been explored. Given the prominent functions of eRNAs and SEs in both physiological and pathological conditions, we hypothesized that eRNAs may contribute to SLE and investigated the hypothesis based on epigenomic profiling.

In this study, we used RNA-sequencing (RNA-seq) data from 20 PBMC samples of 10 SLE patients and 10 healthy controls (HCs) to characterize the expression landscape of eRNAs and associated SEs. The putative functions of the eRNAs identified were investigated using the Functional Annotation of the Mammalian Genome (FANTOM) project,¹⁷ Roadmap Epigenomics project,¹⁸ and DENdb¹⁹ databases, which comprise annotations of numerous regulatory elements, including enhancers. We identified eRNA expression associated with the inflammatory environment in SLE. Additionally, we validated the differentially expressed eRNAs in SLE patients compared with HCs and we discussed the putative functions of the eRNAs in SLE.

Methods

Subjects and Study Design

RNA-seq data from PBMCs from 10 female SLE patients stratified by their disease activity characteristics and 10 female HCs, generated in our previous study,⁶ were used to profile and analyze eRNA expression in this study. Based on the SLE disease activity index (SLEDAI), patients were classified into a stable disease group (SLEDAI < 5, SSG) and an active disease group (SLEDAI ≥ 5, SAG). For validation, 31 SLE patients and 26 age- and sex-matched HCs without arthralgia, heart failure, renal failure, or autoimmune disease, and free from other inflammatory conditions, were recruited in this study. The research protocol was approved by the Medical Ethical Committees of the First Affiliated Hospital of Wenzhou Medical University. All subjects participating in this study provided written informed consent.

PBMC and RNA Isolation and Next-Generation Sequencing

Full protocol details were reported in our previous study.⁶ In brief, sequencing libraries were prepared using a NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, USA) according to the manufacturer’s protocol. Ribosomal (r)RNA was removed from the total RNA using Ribo-Zero™ rRNA Removal kit (Human/Mouse/Rat) (Illumina, San Diego, CA, USA). The libraries were sequenced on an Illumina HiSeq instrument (Illumina), generating 150-bp pair-end reads.

Identification and Quantification of eRNAs

Clean RNA-seq reads were mapped using the Hisat2 software (http://ccb.jhu.edu/software/hisat2).²⁰ Samples with a mapping rate of >80% and a proportion of multiple mapped reads of <20% were considered to be reliable and usable for subsequent analyses. The StringTie software was used to assemble the reads into transcripts and to quantify them based on genome alignment.²¹ Determining the eRNA transcription region, ie, the enhancer region, is essential for eRNA identification. We identified eRNAs by integrating various previously reported strategies for eRNA identification.^9,11,22 The entire process comprised five major steps. (1) Collection of enhancer annotation data. In this study, three databases that use multi-omics methods to infer enhancer regions, ie, FANTOM 5, Roadmap Epigenomics, and DENdb, were used to determine enhancer regions and eRNA transcription regions.^11,17,19 A dataset of putative enhancer regions in human PBMCs was constructed based on the human hg19 reference genome by cross-comparisons in the three databases. (2) Filtering of the enhancer annotation data. Putative enhancer regions were compared with annotated genes in the hg19 reference genome, and enhancer regions overlapping with known genes (including upstream 1-kb and downstream 1-kb regions) were filtered out to prevent noise generated by eRNAs from enhancers too close to promoters. The remaining enhancer regions (with >99% sequence similarity to the original sequence) were defined as hg19 genome eRNA transcription regions¹¹ and were used for subsequent analyses. (3) Transcript assembly and merging. Transcript assemblies from each sample were merged using the Cuffmerge software. (4) eRNA transcript filtering. The Cuffcompare software was used to compare and analyze the assembled transcripts and eRNA transcription regions. Transcripts overlapping with eRNA transcription regions and falling in the antisense strand of the eRNA transcription region were obtained. Next, transcripts that crossed into the coding region and its 1-kb flanking regions due to assembled transcripts being too long were excluded. Thus, candidate eRNA transcripts were obtained. (5) Single-exon screening. Finally, transcripts containing exons were screened out as eRNAs rarely undergo splicing. Based on the eRNA transcript data, corresponding enhancer positions were annotated. A human PBMC eRNA.gtf file containing locus data, enhancer annotation data, and exon data of each eRNA transcript was generated (Supplementary File, Table S1) and was used for subsequent eRNA quantification and differential expression analysis. eRNA expression levels in human PBMCs were calculated as the mean expression value as fragments per kilobase of transcript sequence per millions base pairs sequenced reads per million (FPKM). The Cuffdiff software was used to analyze differential expression, using threshold values of P < 0.05 and fold change >2.²³

Identification of TFs and SEs

Enhancers regulate target gene transcription by binding TFs, transcriptional initiation cofactors, RNA polymerase II, and other proteins involved. Recent studies on eRNAs have shown that TFs are also involved in the regulation of eRNA transcription. Through correlation analysis of the transcription status of eRNAs and related TFs, the correlations between TFs and eRNAs, enhancers, and specific cell types can be investigated.^11,24 Human TF genes were downloaded from the AnimalTFDB.^25,26 Pearson correlation coefficients of single eRNAs and TFs were calculated to obtain a set of TFs co-expressed with eRNAs. Based on the table of critical values for Pearson correlation and the degrees of freedom determined by the number of samples²⁷ (https://conjointly.com/kb/correlation-statistic/), a 0.05 cutoff level for significance for Pearson correlation was selected. Key TFs were selected based on the proportion of eRNAs co-expressed with TFs to all eRNAs. Based on the correlations between eRNAs and TFs, TFs showing a significant correlation with >10% of eRNAs were identified as key TFs. The FIMO software in the MEME suite (http://meme-suite.org/doc/fimo.html?man_type=web) was used to obtain human TF motif data from the Jaspar database (http://jaspar.genereg.net/), and TF-binding sites/motifs were predicted based on the enhancer loci data corresponding to the eRNAs obtained from the analysis.^28,29

SEs are cis-regulatory elements with determinative regulatory effects on cell fate. The ROSE algorithm, a tool for SE identification developed by the Young Lab at MIT (http://younglab.wi.mit.edu/super_enhancer_code.html), uses the quantitative.bam file from chromatin immunoprecipitation (ChIP)-seq data to determine the peak size of binding between TFs or cofactors and a 12.5-kb locus. Through ranking and background calculations, putative SE regions exceeding the background cutoff value are determined.^14,30 A recent study combined eRNA expression levels and the ROSE algorithm to predict SEs.³¹ Therefore, these methods were combined for de novo SE prediction. The SEs identified were annotated using the SEdb and dbSuper databases.

Bioinformatics Analysis of eRNA Targets

eRNAs participate in the regulation of enhancer activation of adjacent target genes and the transcription of nearby target genes in a cis-regulatory manner through their transcriptional behavior. In addition, they can regulate distant target genes in trans by binding with trans-acting elements.⁹ eRNA expression levels are positively correlated with the activity of their enhancers and thus, correlations between eRNA expression levels and gene expression levels can be used to predict target genes of enhancers.¹⁵ Therefore, we predicted target genes of eRNAs based on co-location and co-expression.

As for co-location analysis, enhancers can regulate the transcription and expression of upstream and downstream genes from a distance, and the effective distance based on enhancer activity is generally <1 Mb. Therefore, we set the co-location threshold to 1 Mb up- and downstream of the eRNA, and genes within 1 Mb upstream and downstream were searched based on the location of the eRNA. For eRNA co-location analysis, the most recent version of the BEDTools software was used, with the parameters -D b -a eRNA.bed -b coding-gene.bed. -D b indicates that the distance is based on coding genes; a negative number indicates that the eRNA is downstream of the gene, and a positive number indicates that it is upstream.

Co-expression analysis is based on the principle that if an eRNA has a regulatory relationship with a target gene, their expression levels in different samples and groups should be highly correlated. We used the Python SciPy package to determine correlations between eRNAs and genes by calculating Pearson correlation coefficients, based on which target genes were predicted. The table of critical values for Pearson correlation was used for confirmation, with degrees of freedom = number of samples – 2.²⁷

Next, we performed functional enrichment analysis of the target genes of differential eRNAs to predict the biological functions of the eRNAs. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome functional enrichment analyses were performed on co-located and co-expressed target genes of up- and downregulated eRNAs using the over-representation analysis method in the clusterProfiler package. Additionally, the putative functions of target genes and master TFs of candidate eRNAs were analyzed using the KEGG Orthology Based Annotation System (KOBAS) web software version 3.0 (http://kobas.cbi.pku.edu.cn/kobas3) based on pathway enrichment analysis and DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/tools.jsp), respectively. P_adj < 0.05 or P < 0.05 was regarded significant enrichment.

Quantitative Reverse Transcription PCR (RT-qPCR) Validation of the eRNA Expression Data

Detailed protocols are provided in our previous report.⁶ RT-qPCRs were run on a qTOWER 2.2 Real-Time PCR system (Analytik Jena, Jena, Germany) using reaction mixtures containing 1 μL of diluted cDNA and 5 μL of 2× SYBR^® Green Supermix (DBI, Ludwigshafen, Germany). The thermal cycling conditions were as follows: 95 °C for 3 min followed by 39 cycles of 95 °C for 10s and 60 °C for 30s. Relative eRNA levels in PBMCs were calculated using the following equation: amount of target = 2^−ΔCt, where delta Ct = Ct _eRNAs – Ct _GAPDH. The sequences of the eRNA-specific primers used are provided in Supplementary File, Table S2. Each biological replicate was analyzed in triplicate.

Statistical Analysis

Statistical analysis was conducted using the SPSS 22.0 software (SPSS, Inc., Chicago, IL, USA). Means of two groups were compared using the independent-samples t-test or Mann–Whitney U-test, as relevant. A two-sided P-value < 0.05 was considered statistically significant. The Wilcoxon test was used for differential eRNA expression analysis. Receiver operating characteristic (ROC) curve analysis of PBMC eRNA expression was conducted using MedCalc^® Statistical Software version 19.6.1 (MedCalc Software Ltd., Ostend, Belgium; https://www.medcalc.org; 2020).

Results

eRNA Profiling of PBMCs from SLE Patients and HCs

A flow chart of the transcriptome-wide eRNA profiling analysis is shown in Figure 1. In total, 32,308 eRNA transcripts, comprising 16,594 (51.36%) sense strands and 15,714 (48.64%) antisense strands, were identified in the 20 samples (10 SLE patients and 10 HCs) by eRNA identification using the FANTOM, Roadmap Epigenomics, and DENdb databases (Supplementary File, Table S1). The SLE cohort comprised 10 patients, 5 of whom had stable disease (SSG) and five of whom had active disease (SAG). Analysis of the distribution of eRNA transcripts in the HC, SSG, and SAG groups revealed that eRNAs were transcribed from all chromosomes (Figure 2A and B and Supplementary File, Figure S1A). The median number of eRNAs identified was 14,022 in HC PBMCs, 13,065 in SSG PBMCs, and 11,706 in SAG PBMCs, with no significant difference between the HC group and the SSG and SAG groups, nor between the SSG and SAG groups (Figure 2C). Analysis of the eRNA length distribution revealed that the majority of eRNAs in the HC, SSG, and SAG groups were 200–800 nt (Supplementary File, Figure S1B). The mean eRNA transcript expression levels (in FPKM) in HCs were lower than those in the two SLE groups, especially in the SAG group (P < 0.05) (Figure 2B and D). Mean eRNA transcript expression levels did not significantly differ between the SSG and SAG group (Figure 2D). Next, we classified all eRNAs identified into three groups based on their expression: 5873 ubiquitous eRNAs were expressed in >80% (16/20) of the samples, 22,233 median-ubiquitous eRNAs were expressed in at least two samples, and 2908 special eRNAs were expressed in only one sample (Figure 2E, Supplementary File, Table S3). Special eRNAs accounted for only 9% of the eRNAs detected in PBMCs from the SLE and HC groups, and for 91% of the eRNAs in at least two samples.

Figure 1 Flow chart of eRNA identification.

Figure 2 Expression landscape of eRNAs in PBMCs from healthy controls (HCs) and patients with systemic lupus erythematosus (SLE). (A) Distribution of eRNAs in the chromosomes of HCs, active SLE patients (SAG group), and stable SLE patients (SSG group). The outer circle is a chromosome positional map, and the numbers outside the circle indicate the chromosome numbers. Different colors are used to distinguish different chromosomes. Tick marks indicate chromosome coordinates. Unit are in millions. The inner circle indicates the circular RNA density of each coordinate interval in the chromosomes. Data are shown as a scatter plot, with each point representing a coordinate interval in the chromosome, and values decrease from the outside to the inside. (B) Mean eRNA expression in chromosomes of the HC, SAG, and SSG groups. The lower and upper hinges correspond to the 25th and 75th percentiles. The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge. The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). Outliers are plotted individually. (C) Numbers of eRNAs in the HC, SAG, and SSG groups. Each dot represents an individual sample. (D) Total mean eRNA expression in the HC, SSG, and SAG groups. Each dot represents an individual sample. *P < 0.05; ns, no statistical significance. (E) Expression profile of eRNA in human PBMCs. Blue, yellow, and green bars indicate ubiquitous, semi-ubiquitous, and unique eRNAs, respectively. Purple, red, and blue-green bars denote the HC, SAG, and SSG groups, respectively. The red color gradient represents the eRNA expression level (in FPKM). Abbreviation: ns, no statistical significance.

SE Profiling of PBMCs from SLE Patients and HCs

SEs, an emerging and distinct class of enhancers, have broader and more profound effects in pathological and physiological conditions than regular enhancers,³² and they drive master genes crucial for cell identity. Therefore, we explored the landscape of SEs in PBMCs from SLE patients and HCs. In total, 7068 functional SEs (Supplementary File, Table S4) were detected in the total samples. In Supplementary File, Figure S2A, they are ranked according to their expression levels in SSG and SAG patients and HCs. We annotated the identified functional SEs using the SEdb and dbSuper databases (Supplementary File, Table S5). Intriguingly, each SE was generally associated with one to several eRNAs, whereas each eRNA was related to only one SE (Supplementary File, Table S5). For example, one SE was associated with four eRNAs, including TCONS_00378828, TCONS_00392560, TCONS_00378829, and TCONS_00378832, whereas TCONS_00034326 was related to only one SE (chr1:89238617–89240917). Analysis of the length distribution of the SEs revealed that the majority of SEs in the HC, SSG, and SAG groups were 2000–9000 nt, and the number of SEs was higher in HCs than in the two SLE groups (Supplementary File, Figure S2B). We classified the SEs identified into three groups: 16 ubiquitous SEs were expressed in >80% (16/20) of the samples, 4387 semi-ubiquitous SEs were expressed in at least two samples, and 2665 unique SEs were expressed in only one sample (Supplementary File, Figure S2C, Table S4). In contrast to the findings for eRNAs, ubiquitous SEs accounted for only 0.2% of the SEs detected in PBMCs from SLE patients and HCs, and for 38% of SEs in only one sample. Intriguingly, the number of SE transcripts was lower in both SLE groups than in HCs, especially in the SAG group (P < 0.05), whereas there was not significant difference between the SSG and SAG groups (Supplementary File, Figure S2D). The majority of SEs were associated with 1–10 eRNA(s), and the number of SEs associated with different eRNAs decreased with increasing disease activity (Supplementary File, Figure S2E).

Analysis of SLE-Related Master TFs of eRNAs

SEs consist of multiple enhancer-like elements that are densely occupied by master TFs.¹⁴ TFs play an important role in mediating eRNA biogenesis;⁹ however, SLE-related TFs of eRNAs have not been fully unraveled. We collected SLE-related TFs from the AnimalTFDB^25,26 and evaluated Pearson correlations between individual TF and eRNA expression levels in the SLE PBMC samples. Putative master regulators were defined as those significantly correlated with ≥10% of individual eRNAs. Based on computational predictions, we found 281 SLE-related putative master TFs, including several well-known regulators, such as JUNB, JUN, and STAT5 (Figure 3A, Supplementary File, Table S6). JUNB is lowly expressed in SLE and induces interleukin-6 secretion,³³ suggesting its potential importance in eRNA regulation in SLE.

Figure 3 Putative regulation of eRNA biogenesis in SLE PBMCs. (A) Putative regulators of eRNAs in SLE PBMCs. Each dot represents one transcription factor (TF). Red dots denote putative master TFs significantly correlated to ≥10% of individual eRNAs. The X-axis indicates the chromosomes and the y-axis the percentage of eRNAs correlated to each TF. (B) Gene Ontology (GO) analysis of putative SLE-related master TFs. Biological processes are indicated in red, cellular components in Orange, and molecular functions in black. (C) KEGG pathway analysis of putative SLE-related master TFs. The Y-axis indicates pathway names and the X-axis the richness factor.

Next, we performed GO analyses of the identified TFs. In the biological process category, TFs were mostly enriched in transcriptional processes, in the CC category, they were mostly linked to the nucleus and TF complex, and in the molecular function category, they were mostly associated with nucleic-acid-related molecular binding functions (Figure 3B, Supplementary File, Table S6). KEGG analyses revealed that the TFs, including CREB3, JUN, and JUNB, were significantly enriched in the cGMP-PKG and tumor necrosis factor (TNF) signaling pathways (Figure 3C, Supplementary File, Table S6). KEGG pathway analysis results obtained using DAVID are detailed in Supplementary File, Table S6. The master regulators were generally enriched in functions related to SLE-related pathways, which provides a potential explanation for the high mean expression levels of eRNAs observed in the SLE groups.

Enhancers activate their target genes by binding TFs, transcriptional initiation cofactors, and related proteins. TFs have a domain that binds to a specific DNA sequence, and the sequence on the corresponding enhancer is known as the TF-binding site or conserved domain (motif). By predicting the motifs in a TF binding an enhancer, the binding sites between specific TFs and the enhancer locus and their probabilities can be inferred, allowing prediction of their regulatory relationship in the cell. Putative enhancer motifs in SLE-related master TFs (q < 0.05) were predicted and are listed in Supplementary File, Table S7.

Differentially Expressed eRNAs in PBMCs from SLE Patients and HCs

To find key eRNA biomarkers of SLE, we determined eRNAs differentially expressed in the SSG vs HC, SAG vs HC, and SLE (SSG and SAG) vs HC groups. We found 66 up- and 45 downregulated eRNAs in SLE vs HCs (Figure 4A and B, Supplementary File, Figure S3A and B). There were 109 up- and 149 downregulated eRNAs in SSG vs HCs (Figure 4A and B, Supplementary File, Figure S2C and D), and 105 up- and 69 down-regulated eRNAs in SAG vs HCs (Figure 4A and B, Supplementary File, Figure S3E and F). Based on the intersections in a Venn diagram of the differentially expressed eRNAs from comparisons of the HC vs the SLE, SAG, and SSG groups, we were able to narrow down the upregulated candidate SLE-related eRNAs to 13 and the downregulated candidate SLE-related eRNAs to 23 (Figure 4A and B). Clustering heatmaps of these eRNAs are shown in Figure 4C and D. Importantly, 10 of the 13 upregulated candidate SLE-related eRNAs and 17 of the 23 downregulated candidate SLE-related eRNAs were associated with functional SEs, which suggests that these candidate eRNAs have key roles in SLE pathogenesis.

Figure 4 Differentially expressed eRNAs in PBMCs from HCs and SLE patients. (A and B) Venn diagram showing the overlap between significantly differentially expressed eRNAs in PBMCs of SLE patients (SAG and SSG) vs HCs, SAG vs HCs, and SSG vs HCs. (C and D) Hierarchical clustering of the 13 commonly upregulated and 23 commonly downregulated eRNAs. Red represents relatively high expression and green relatively low expression.

Pathway Analysis of the Targets of Candidate SLE-Related eRNAs

Based on the eRNA transcriptome data, we explored the functions of putative co-expressed and co-located target genes of the candidate SLE-related differentially expressed eRNAs using GO, Reactome, and KEGG pathway analyses. In GO analysis, co-expressed targets of the 13 upregulated eRNAs were significantly enriched in cellular response to type I interferon (IFN) and the IFN-related signaling pathway, whereas those of the 23 downregulated eRNAs were markedly enriched in cytosolic ribosome and nuclear-transcribed mRNA catabolic process (Supplementary File, Figure S4A). In Reactome pathway analysis, the co-expressed targets of the 13 upregulated eRNAs were involved in IFN-related regulation signaling, cytokine signaling, and DNA methylation, whereas those of the 23 downregulated eRNAs were associated with eukaryotic translation-related regulation signaling (Supplementary File, Figure S4B). In KEGG pathway analysis, the co-expressed targets of the 13 upregulated eRNAs mostly participated in the systemic lupus erythematosus pathway, cytokine-cytokine receptor interaction, and Toll-like receptor signaling pathway, whereas those of the 23 downregulated eRNAs were involved in RNA degradation and other metabolism-related regulatory processes (Supplementary File, Figure S4C). Interestingly, the co-expressed targets of candidate upregulated eRNAs were largely associated with SLE-related pathways, including IFN-related, systemic lupus erythematosus, and cytokine signaling.

According to GO analysis, the co-located targets of the 13 upregulated eRNAs were markedly enriched in pyrophosphatase activity and other binding-related processes (Supplementary File, Figure S4D). Reactome analysis revealed that the co-located targets of the candidate upregulated eRNAs were involved in gene expression and DNA methylation (Supplementary File, Figure S4E), and KEGG pathway results indicated that the co-located targets of the candidate upregulated eRNAs participated in the systemic lupus erythematosus pathway (Supplementary File, Figure S4F). In contrast, the co-located targets of the candidate downregulated eRNAs were not significantly enriched in any pathway. Based on the data, both co-expressed targets and co-located targets of upregulated eRNAs were likely to participate in SLE-related pathways. KEGG analysis of the integrated co-expressed and co-located targets of the candidate upregulated eRNAs revealed that they were commonly involved in systemic lupus erythematosus (P_adj < 0.05) as shown in Supplementary File, Figure S5. Most of these are involved in autoantigens, cytokine–cytokine receptor interaction, immune complexes, and tissue injury and damage. Together, these results suggested that the targets of upregulated rather than those of downregulated eRNAs play a critical role in regulating SLE-related signaling.

Validation of Candidate Upregulated eRNAs

To evaluate the clinical utility of the upregulated candidate SLE-related eRNAs, we measured their expression in a validation cohort consisting of six SLE patients and six HCs. Among the 13 candidate eRNAs, four eRNAs were excluded based on the following criteria: 1) eRNA was not detected in several HC samples; 2) eRNA > 10,000 bp. Additionally, RT-qPCR of TCONS_00053400 failed. Thus, only eight upregulated eRNAs (TCONS_00034326, TCONS_00416100, TCONS_000219381, TCONS_00315667, TCONS_00359113, TCONS_00378832, TCONS_00463820, TCONS_00378833) were evaluated in the first validation cohort (Figure 5A, Supplementary File, Figure S6A). Based on eRNA expression being upregulated in the majority of samples of the first validation cohort, TCONS_00034326, TCONS_00416100, and TCONS_000219381 were selected for validation in another, larger cohort, including 25 SLE patients and 20 HCs. Only TCONS_00034326 was significantly increased in expression in SLE patients in the second validation cohort (P < 0.05; Figure 5B, Supplementary File, Figure S6B). To confirm the potential utility of TCONS_00034326 as a diagnostic biomarker, we conducted ROC curve analysis. The AUC of TCONS_00034326 for distinguishing SLE patients from HCs was 0.691, with a diagnostic sensitivity of 0.833 and specificity of 0.55 (Figure 5C).

Figure 5 Expression of candidate SLE-related upregulated eRNAs. (A) Expression of three selected eRNAs in a validation cohort comprising six SLE patients and six HCs as indicated by RT-qPCR. The three eRNAs were selected for analysis in a larger validation cohort. Data are presented as 2−ΔCt relative to GAPDH expression (mean ± SEM). (B) Expression of TCONS_00034326 in a large validation cohort comprising 25 SLE patients and 20 HCs. One outlier in the SLE group was excluded. *P < 0.05. (C) ROC curve of TCONS_00034326 in the large validation cohort for SLE diagnosis.

Potential Functions of TCONS_00034326 in SLE

Functional analysis of TCONS_00034326 in SLE was conducted using KOBAS version 3.0. Pathway analysis indicated that the targets of TCONS_00034326 were significantly enriched in vital SLE-related pathways, including TNF signaling, immune signaling, metabolic pathways, and NF-kappa B signaling (Figure 6A, Supplementary File, Figure S7A–C, Table S8). TCONS_00034326 was strongly related to an SE (chr1:89238617–89240917), which is more powerful in regulating gene expression in SLE than the eRNA (Supplementary File, Table S4). To investigate the potential biological function of TCONS_00034326 further, we applied two methods to identify TCONS_00034326-related master TFs. On the one hand, considering that TCONS_00034326 was associated with the SE in chr1:89238617–89240917, we combined TCONS_00034326-related TFs with SE-associated TFs and found three TFs, including SP1, SPI1, and ZNF467, which activate the transcription of various genes (Figure 6B).^34–36 On the other hand, we obtained 26 TCONS_00034326-related TFs (|PCC| > 0.8). Based on a Venn diagram of TCONS_00034326-related TFs and SLE-related master TFs, we found 24 common candidate TFs (Figure 6B). Pathway analysis of these 24 TFs showed that they are involved in TNF signaling and IFN-α/β and -γ signaling (Figure 6C, Supplementary File, Table S9). Notably, three of the 24 (IRF1, CREB5 and JUN) TFs participate in the TNF signaling pathway, which is involved in mediating immunopathogenic injury in SLE.^37,38 Two of the 24 TFs (IRF1 and IRF2) are involved in signaling of not only IFN-α/β, but also IFN-γ, which are inducers and effectors in SLE pathogenesis.^39,40 Altogether, these findings indicate that TCONS_00034326 and its related master TFs may be of great significance in modulating SLE pathogenesis.

Figure 6 Putative regulation and regulatory network of eRNA biogenesis in SLE. (A). Putative regulatory pathways of TCONS_00034326 target genes in SLE. Three well-known SLE-related pathways (TNF signaling, immune system and metabolic pathways, and NF-κB signaling) are detailed in Supplementary File, Table S8 and Figure S7. (B) Overlap of TCONS_00034326-related TFs and SE-associated TFs (upper panel) and of TCONS_00034326-related TFs and SLE-related master TFs (lower panel). (C) Putative regulatory pathways of 24 TCONS_00034326-related master TFs in SLE. Three master TFs participate in TNF signaling, two master TFs are involved in interferon α/β and γ signaling (Supplementary File, Table S9).

Discussion

Recently, the importance of eRNAs and SEs in the regulation of gene transcription has been increasingly recognized based on their wide transcription in various conditions. Some eRNAs transcriptome studies in human cancers^11,41 and neuronal systems⁴² have highlighted the unique role of eRNAs in gene regulation in human tissues and organs and demonstrated that eRNAs have value as potential therapeutic targets. In this study, we utilized a novel computational framework to reveal the global expression landscapes of eRNAs in PBMC samples from SLE patients and HCs.

Zhang et al identified cancer-specific eRNAs by integrating multi-omics data from the large-scale ENCODE, FANTOM, and Roadmap Epigenomics projects.¹¹ Carullo et al combined “assay for transposase-accessible chromatin”-seq and RNA-seq, which enabled genome-wide enhancer identification and putative enhancer–gene pair prediction based on transcriptional output correlations.⁴² Similarly, by integrating enhancer annotation data from the FANTOM, Roadmap Epigenomics, and DENdb databases, we gained novel insights into the eRNA and SE expression landscapes of human PBMCs. We utilized these three project databases to identify eRNAs as comprehensively and accurately as possible, to ensure that our results are reliable. eRNAs were first detected in large-scale transcriptome profiling,⁴³ and RNA polymerase II ChIP-seq⁴⁴ studies. Subsequently, in a study of genome-wide eRNA synthesis, transcripts from enhancers were detected by RNA-seq using total RNAs, but not polyadenylated RNAs alone.⁴⁵ Polyadenylated RNA-seq detects accumulated RNA end products rather than the dynamic, transient transcriptional activity of a cell, and is not sufficiently robust to detect eRNAs, which generally lack poly(A) tails.³² Compared to other tools used to detect eRNAs, RNA-seq is a well-established and simple tool that allows high-throughput analysis; however, it detects RNA end products rather than the dynamic transcriptional activity, and most reads are from highly expressed coding genes.⁹ Therefore, we utilized genome-wide transcriptomics along with strand-specific library construction and multiple bioinformatics analyses to study the eRNA composition of the PBMC population after ribosomal RNA depletion to obtain complete eRNA profiles.

We observed that eRNAs were transcribed from all chromosomes in HCs and SLE patients, and more than 90% of eRNA genes were detected in at least two samples from SLE patients and HCs, whereas 9% of eRNAs were detected in only one sample, which is slightly different from findings in a cancer study.¹¹ Consequently, RT-qPCR did not detect the unique eRNAs in some of the samples, which may be due to the very low or lack of expression of these eRNAs. eRNAs are retained in the nucleus, and they are very unstable and therefore, can degrade within minutes.⁸ The total numbers of eRNAs did not significantly differ between the SLE and HC groups, but the mean eRNA expression was increased in the SAG group when compared with that in the HC group. The recently reported heart-specific eRNA HERNA1 promotes pathological cardiac metabolism and contractile dysfunction, and knockdown of HERNA1 protected mice from heart disease and even reversed cardiac pathology in mice.⁴⁶ Likewise, increased eRNA expression is likely to be correlated with SLE pathology.

SEs have been shown to recruit transcription complexes and generate eRNAs. Consistent with a previous report,¹⁰ we found that SEs were generally substantially longer than the majority of eRNAs. The number of SEs was negatively correlated with disease activity, which implicates that an imbalance in SEs promotes SLE pathogenesis. In contrast to the findings for eRNAs, unique SEs accounted for nearly 38%, whereas ubiquitous SEs represented only 0.2% of the SEs identified in the total SLE and HC PBMC samples. The large proportion of unique SEs may be due to disease heterogeneity and differences in disease activity among the patients, or the small sample size in this study. Accumulating evidence suggests that eRNAs are not only associated with SEs, but also play functional roles in SE biology.⁸ The fact that SE transcript numbers were negatively correlated with disease activity and eRNA expression levels were positively correlated with disease activity may suggest that part of the eRNAs play a key role in the progression of lupus. One study revealed that differential eRNA expression in patients with autoimmune disease was especially localized to SEs near known genetic variants for autoimmune disease risk.⁴⁷ Thus, SEs appear to play an important role in the biology of SLE. As regulators of the transcriptional machinery, eRNAs can recruit and regulate TFs and promote enhancer and gene transcription.⁴⁸ We identified various SLE-specific TFs as potential regulators of eRNA biogenesis.²⁴ Specific master TFs widely activate the gene expression program to define cell identity, and contribute to multiple diseases.^14,49 In accordance herewith, we observed that putative SLE-related master regulators of eRNAs generally displayed a significant enrichment in functions related to transcriptional misregulation, suggesting a potential mechanistic relation between eRNA biogenesis and genome instability. In line herewith, in a recent report,¹¹ putative master regulators of eRNAs in cancer were generally involved in modulating genomic instability. Interestingly, we found that many SLE-related TFs were involved in modulating the TNF signaling pathway, which mediates immunopathogenic injury in SLE, greatly expanding our knowledge of eRNA biogenesis in SLE.

By classifying the samples according to SLE disease activity, we identified 13 and 23 eRNAs that were commonly upregulated and downregulated, respectively, in active and stable disease. For comprehensive and integral pathway analysis of the candidate SLE-related differentially expressed eRNAs, we used the GO, Reactome, and KEGG pathway analysis tools. GO and Reactome analyses showed that co-expressed target genes of the upregulated eRNAs were largely involved in IFN signaling, which has key immune effector functions in SLE.^50,51 Intriguingly, KEGG analysis indicated that co-expressed and co-located target genes of the upregulated candidate SLE-related eRNAs participate in systemic lupus erythematosus pathways, including IFN-γ, which corroborates that eRNAs function in SLE pathogenesis. Importantly, ten of the 13 candidate upregulated eRNAs were related to SEs. Accumulating evidence suggests that eRNAs can provide insights into diagnosing, prognosticating, treating, and understanding the pathogenesis of diseases that cannot be achieved with current biomarkers.⁸ Considering the putative biological functions of the candidate eRNAs, we further explored their potential functionality as diagnostic biomarkers of SLE through validation analysis. TCONS_00034326 expression was markedly increased in SLE patients compared with HCs. We used ROC curve analysis to confirm the usability of TCONS_00034326 as a biomarker. We demonstrated the SLE-specific expression of certain eRNAs, suggesting that these may be powerful diagnostic and/or prognostic markers in SLE.

TCONS_00034326 was strongly related to an SE (chr1:89238617–89240917), which is more powerful in regulating gene expression in SLE than an eRNA (Supplementary File, Table S5). Using two different methods, we identified 26 TCONS_00034326-related SLE master TFs, including SP1, SPI1, ZNF467, IRF1, and IRF2. SP1 has been implicated in inflammation and genome instability, as well as in epigenetic silencing.⁵² Using the KOBAS software and DAVID pathway enrichment, we further explored the biological roles of the target genes and master TFs of TCONS_00034326. Pathway enrichment analysis indicated that both the target genes and master TFs of TCONS_00034326 were relevant to the regulation of SLE pathological pathways, including TNF signaling and IFN-α/β and γ signaling. TNF signaling is involved in various immunoregulatory and proinflammatory processes in SLE.⁵³ To some degree, TNF apparently has dual effects in checking autoimmunity and fostering inflammation.⁵⁴ In fact, anti-TNF-induced lupus is an adverse effect related to the use of TNF blockers in autoimmune arthritis, with manifestations ranging from mild cutaneous to serious fatal conditions. Type-I IFNs (IFN-α/β) are regarded key IFNs in the pathogenesis of SLE,⁵⁵ and mounting evidence suggests that IFN-γ also plays an essential role in SLE development and severity.⁴⁰ Our data indicate that TCONS_00034326 plays an important role in SLE pathogenesis by recruiting these master TFs. Further, these results suggest that appreciable expression of eRNAs are clinically valuable and inform and motivates future investigation of the molecular mechanisms of eRNA in the immunopathogenesis of SLE.

This study had some limitations. First, our sample size was relatively small. A larger population for validation is required to further explore the clinical value of eRNAs, including TCONS_00034326. Second, the mechanisms of the eRNAs identified in regulating SLE were not revealed. The SEs in our study were identified from analytical databases and may not represent the actual SEs in SLE, and further experimentation is needed to validate these results. We are currently developing a roadmap to evaluate the undetermined mechanisms of eRNA-mediated SLE-related regulation.

Conclusion

Our study portrayed transcriptional eRNA landscapes in PBMCs from SLE patients and HCs and identified several eRNAs that are relevant to SLE. We suggest these may serve as prospective diagnostic and therapeutic targets in SLE pathogenesis.