Identification of genes and pathways related with cerebral small vessel disease based on a long non-coding RNA-mediated, competitive endogenous RNA network.

This study aimed to investigate the cerebral small vessel disease (SVD), an intrinsic disorder of the brain's perforating arterioles, which could cause serious cognitive decreasing and trigger dementia. In this study, we present a multi-step computational approach to construct a functional SVD long non-coding RNA (lncRNA)-mediated ceRNA network (LMCN) by integrating genome-wide lncRNA and mRNA expression profiles, miRNA-target interactions, and functional analyses. We used hypergeometric test to evaluate enrichment significance of miRNAs and we obtained the LMCN, which contained 27 lncRNAs, 7,229 mRNA, and 28,871 lncRNAs-mRNA interrelationship pairs. What's more, co-expression analysis was utilized to constructe a competitive endogenous RNAs (ceRNAs) interaction network which comprised of 21 lncRNAs, 129 mRNAs and 141 interaction pairs. We determined that metastasis-associated lung adenocarcinoma transcript 1 and MIR155 host gene acted synergistically to regulate mRNAs in a network module of the competitive LMCN. Moreover, 7 genes were obtained through Gene Ontology enrichment. In addition, 29 mRNA enriched pathways significantly associated with lncRNAs was obtained via Fisher test. In conclusion, we identified 7 potential lncRNAs and 29 possible lncRNA-mediated pathways associated with SVD. Thus, ceRNAs could accelerate biomarker discovery and therapeutic development in SVD.

Introduction

Cerebral small vessel disease (SVD), an intrinsic disorder of the brain's perforating arterioles, which could induce vascular dementia, cognitive decline or disability (1,2). SVD refers to a set of pathological processes the two most frequent types of which were cerebral amyloid angiopathy and hypertension-associated small vessel diseases (3–5). The damages of SVD on brain are usually lacunar infarcts, microbleeds and white matter lesions. Common SVD can be detected via neuroimaging, whereas small vessels (including small arteries and venules) cannot be diagnosed by neuroimaging. In this condition, postmortem plays an important role in detecting microvascular pathology (6,7). Owing to the development of high-throughput sequencing technology, long non-coding RNAs (lncRNAs) have been discovered in a wide range of biological processes. The functions of lncRNAs in SVD are not well characterized, and the identification of lncRNA biomarkers is challenging. lncRNAs compete with competitive endogenous RNAs (ceRNAs) by acting as miRNA sponges to regulate expression level of other transcripts (8). Several ceRNA related databases have also been developed to facilitate interference of lncRNA function. Collectively, this data underscores the importance of lncRNA interactions with ceRNAs, and indicates that integration of expression profiles and network analysis could enable identification of risk lncRNAs and the underlying tumor pathology of SVD. Up to now, the pathogenesis of SVD is basically unknown. The pathogenesis of SVD is still not clear, so we need to find molecular markers to prevent.

It is reported that human gene transcripts which encode proteins account for approximatly 2% of all transcripts, and the remaining transcripts are non-coding RNAs (ncRNAs), which based on length could divide into two subclasses lncRNA and small ncRNA (9,10). lncRNAs, which often longer than 200 nucleotides, lack open reading frames compared with mRNAs and through chromatin structural alternation and transcriptional modulation regulate gene transcriptions (11–14). Recent studies have identified that lncRNAs played a pivotal role in tumorigenesis and tumor metastasis (15). Prostate cancer (PC) associated loci were mostly located at lncRNA regions. Jin et al reported that 52 loci were enriched in lncRNA regions and SNP rs3787016 which was located at lncRNA regions was risk-associated with PC (16). In human ovarian cancer SKOV3.ip 1 cell line, lncRNA 578 was downregulated while lncRNA 583 was upregulated, compared with parental SKOV3 cell line respectively (17).

ceRNAs modulate gene expression in trans and account for major parts of gene regulators (18). MicroRNA recognition elements (MREs), a common sequences of ceRNAs, which recognized by microRNAs through sequence complementarity and usually induce gene repression (19). Complicated interaction network between gene transcripts have MREs which regulate each other are called ceRNA crosstalk (20). Zhou et al constructed a breast cancer-related ceRNA network by integrating miRNA expression data and mRNA expression data (21). In human ovarian cancer, 10 lncRNAs which related to ceRNA network were determined as potential biomarkers (22).

According to the criteria reported by Cao et al (23), we constructed a SVD associated ceRNA network by integrating mRNA and lncRNA profiling, which contained 236,952 miRNA-mRNA interactions and 359 lncRNA-miRNA intersections. We used a hypergeometric test to evaluate enriched significance of miRNAs related with lncRNAs and mRNAs. A significant co-expressed lncRNA-mediated ceRNA network (LMCN) was constructed which contained 21 lncRNAs, 129 mRNA and 141 ceRNA interaction pairs. In the prediction of lncRNA functions, 4 significant gene oncology terms contained 7 genes were identified. In adition, we identified 29 mRNAs enriched pathways which were significantly related with lncRNAs. These results demonstrated that LMCNs were useful for identifying potential biomarkers associated with SVD.

Materials and methods

Data processing

miRNA-mRNA interactions and lncRNA-miRNA intersections which associated with SVD were downloaded from starBase 2.0 database (http://starbase.sysu.edu.cn/), whilie mRNA and lncRNA expression profiling were downloaded from ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). The mRNA and lncRNA expression profiling comprised of 30 pairs of human SVD and normal control samples from EMBL-EBI dataset whose dataset number was E-MTAB-3408, which employed mapping between probes and genes to perform.

Genes of expression profiling were intersected with mRNAs of miRNA-mRNA intersections and lncRNAs of lncRNA-miRNA intersections respectively. A new expression profiling which contained 7,649 mRNAs and 27 lncRNAs was identified.

Identification of potential ceRNA interactions

Hypergeometric test was used to evaluate enriched significance of miRNAs which were interacted both with lncRNAs and mRNAs. The P-values were subjected to false discovery rate (FDR) method, which were calculated as follows:

A FDR <0.01 was considered as the threshold.

Identification of sub-LMCN

Co-expression analysis of lncRNA-mRNA interactions was carried out to identify lncRNA-mRNA pairs. In the control and the disease groups, Pearson correlation coefficients of interaction pairs were calculated. The Pearson correlation coefficient was calculated as follows:

cov (X,Y) represents the covariance of variables X and Y. σ and σ are the standard deviations of X and Y, respectively. lncRNA-mRNA interactions with Pearson correlation coefficients >0.8 were considered as significant co-expression ceRNA interactions.

Gene Ontology (GO) prediction of mRNAs

Functional enrichment analysis of mRNAs, which utilized Gene Ontology Consortium (http://www.geneontology.org), significantly correlated with lncRNAs, which predicted the function of lncRNAs using the function of mature mRNAs. lncRNA functions were predicted by mature mRNA functions. P-value =0.05 was used as the threshold to perform functional enrichment analysis for mRNAs of sub-LMCN.

Pathway prediction of lncRNAs

The pathway database selected for this analysis is the Reactome (https://reactome.org), which uses Fisher's test to identify lncRNA enriched pathways. Pathways with P-values <0.05 were considered as mRNA enriched pathways which maybe mediated by lncRNAs. These pathways were putative lncRNA-mediated pathways.

Results

Identification of miRNA-mRNA interactions and lncRNA-miRNA intersections

SVD associated mRNA and lncRNA expression profiling data which were downloaded from ArrayExpress database contained 30 pairs of disease and control samples and we obtained a profile data of 12,628 genes were identified via mapping between probes and genes.

A new gene expression profile of 7,676 genes including 7,649 mRNA and 27 lncRNAs was obtained by crossing the genes in the expression profile data with the mRNAs and lncRNAs in the above miRNA-mRNA interactions and lncRNA-miRNA interactions, respectively. Furthermore, the interactions containing new-expressed gene (mRNA or lncRNA) were screen out from miRNA-mRNA interactions and lncRNA-miRNA interactions, and we obtained 236,952 pairs of miRNA-mRNA interactions and 359 pairs separately.

We used hypergeometric test to evaluate enrichment significance of miRNAs which were interacted with lncRNAs and mRNAs. The P-values were corrected by false discovery rate method and the LMCN was filtrated from the lncRNA-mRNA interactions with P-values <0.01. As a result, the LMCN contained 27 lncRNAs, 7,229 mRNA, and 28,871 lncRNAs-mRNA interrelationship pairs.

Identification of a highly competitive sub-LMCN

The partial sub-networks revealed a more detailed picture of how the lncRNAs synergized with competing mRNAs, while the LMCN could provide a global view of all possible competing ceRNA interactions. Therefore, we derived a high-competing lncRNA-related sub-network (sub-LMCN) from the LMCN by co-expressed analysis the lncRNA-mRNA interactions screened above. Pearson correlation coefficient threshold of lncRNA-mRNA interactions in the control and disease groups were calculated respectively. lncRNA-mRNAs interactions with Pearson correlation coefficient threshold >0.8 were considered as significant co-expressed ceRNA interactions. These ceRNAs constructed a network which contained 21 lncRNA, 129 mRNA and 141 interaction pairs (Fig. 1).

Figure 1.

Network of lncRNA-mRNA co-expression analysis. Orange nodes in the figure represent mRNA while the blue node indicates lncRNA and the larger the degree is, the larger the node size is. The larger the degree is, refers to the more lncRNA-regulated mRNA is, indicating that the lncRNA plays a greater role in regulating the network. lncRNA, long non-coding RNA.

What's more, a cooperative with competition network was identified from the sub-LMCN. SVD associated lncRNAs such as MIR155 host gene (MIR155HG) and metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) competed with other mRNAs in the sub-LMCN (Fig. 2).

Figure 2.

Collaborative competition module of sub-LMCN. The blue node represents lncRNA and the orange node represents mRNA. And the edge value refers to the Pearson correlation coefficient for each edge. LMCN, lncRNA-mediated ceRNA network; lncRNA, long non-coding RNA.

Prediction of lncRNA functions

The functions of mature mRNAs were used to predict lncRNA functions. mRNAs were enriched via GO analysis and we identified 4 GO terms, including snRNA metabolic process (P=0.0037), snRNA transcription (P=0.0058), snRNA transcription from RNA polymerase II promoter (P=0.0058) and ncRNA transcription (P=0.0280) (Fig. 3). The snRNA metabolic process including 7 genes (RPRD2, RPRD1A, NOP10, POLR2A, SNAPC3, POLR2H, TAF13) and each of the remaining three terms contain 6 genes (RPRD2, RPRD1A, POLR2A, SNAPC3, POLR2H, TAF13).

Figure 3.

GO predictions of mRNAs. The ordinate indicates the terms enriched by GO, the abscissa represents the number of genes enriched by GO, and the color of each column represents P-value. GO, Gene Ontology.

Fisher test was used to identify lncRNA enrichment pathways, which were selected from Reactome database. As results, 29 enriched pathways of mRNA enrichemt significantly correlated with lncRNA were identified, which maybe mediated by lncRNAs (Table I).

Table I.

Twenty-nine lncRNA-mediated pathways were identified by Fisher test.

Pathway names	Corrected P-values
Processing of capped intron-containing pre-mRNA	0.024321
Elongation arrest and recovery	0.033775
Gene expression	0.033775
Nuclear import of Rev protein	0.033775
Pausing and recovery of HIV elongation	0.033775
Pausing and recovery of Tat-mediated HIV elongation	0.033775
RNA polymerase II pre-transcription events	0.033775
Transcriptional regulation by small RNAs	0.033775
Transport of ribonucleoproteins into the host nucleus	0.033775
Transport of the SLBP dependant mature mRNA	0.033775
Transport of the SLBP independent mature mRNA	0.033775
Viral messenger RNA synthesis	0.033775
Formation of RNA Pol II elongation complex	0.034444
Influenza life cycle	0.034444
Interactions of Rev with host cellular proteins	0.034444
mRNA splicing - minor pathway	0.034444
RNA polymerase II promoter escape	0.034444
RNA polymerase II transcription elongation	0.034444
RNA polymerase II transcription initiation	0.034444
RNA polymerase II transcription initiation and promoter clearance	0.034444
RNA polymerase II transcription pre-initiation and promoter opening	0.034444
RNA polymerase III abortive and retractive initiation	0.034444
RNA polymerase III transcription	0.034444
Transcription	0.034444
Transport of mature mRNA derived from an intronless transcript	0.034444
Transport of mature mRNAs derived from intronless transcripts	0.034444
Antigen processing: ubiquitination and proteasome degradation	0.034552
Influenza infection	0.036016
Regulatory RNA pathways	0.049813

lncRNA, long non-coding RNA.

Discussion

ceRNAs mutually regulate miRNAs expression through competing mechanisms, which are important for physiological and pathological processes of a variety of cancers. It has developed several valuable ceRNA resources to promote the lncRNAs functional analysis such as DIANA-LncBase (24), starBase (25) and LncACTdb (26). Here, we have used a more comprehensive approach to construct a functional LMCN across 30 pairs of samples. In our study, we constructed a ceRNAs network by integrating lncRNA and mRNA co-expression profiling, which contained 27 lncRNAs, 7,649 mRNA and an expression data of 7,676 genes. And then, we obtained 236,952 pairs of miRNA-mRNA interactions and 359 pairs of lncRNA-miRNAs intersections through extracting from miRNA-mRNA interactions and lncRNA-miRNA interactions separately. Hypergeometric test was used to evaluate enrichment significance of miRNAs, which has interaction with both lncRNA and mRNA. The P-values were corrected by FDR method, and lncRNA-mRNA intersections (LMCN) with corrected P-values <0.01 were screened, which was containing 27 lncRNAs, 7,229 mRNAs and 28,871 pairs of LMCNs. In addition, co-expression analysis of lncRNA-mRNA interactions screened above, which was to calculate the Pearson correlation coefficients of these interactions in the control and disease groups. The interactions between Pearson correlation coefficients larger than 0.8 were screened for ceRNA interactions that were significantly co-expressed and network maps were constructed, which contained 21 lncRNAs, 129 mRNAs and 141 pairs of interactions. However, the lncRNA and miRNA with most nodes in the interaction network and the correlation with patents prognosis should be further analyzed, we will add the analysis in our next study.

By exploring the highly competitive sub-network of the LMCN, we identified a co-competing module from sub-LMCN that MALAT1 and MIR155HG comprised a ceRNA regulating module that competed with RPRD1A, RPRD2, and AGGF1. lncRNA MALAT1, a well-known lncRNA, is associated with pathogenesis and progress of variety of tumors (27). MIR155HG, the primary microRNA of miR-155, which plays an essential role in hematopoiesis, inflammation, and tumorigenesis (28). Blood-based biomarkers are critical for prediction of SVD patient survival. We observed MALAT1 expression in white blood cells in a tissue-specific RNA analysis (29) and in several blood-associated RNASeq profiles. These data indicated that MALAT1 can be secreted into the circulation and is therefore a potential blood-based biomarker for SVD.

Functional enrichment analysis of mRNAs significantly correlated with lncRNAs, which predicted the function of lncRNAs using the function of mature mRNAs. Seven genes, including RPRD2, RPRD1A, NOP10, POLR2A, SNAPC3, POLR2H, TAF13, were obtained through GO enrichment. In addition, 29 mRNA enriched pathways significantly associated with lncRNAs was obtained via Fisher test. The results demonstrated that the cancer-associated LMCN can be used to accelerate biomarker discovery and SVD therapeutic development. However, lack of KEGG analysis is the limitation of our study, we will offer KEGG analysis results in our future manuscript. What's more, we will evaluate our predictions using clinical examples in future study.

In conclusion, we identified 7 potential lncRNAs and 29 possible lncRNA-mediated pathways associated with SVD through enrichment. These findings may contribute to deepen understanding of SVD mechanisms.