1.
Introduction
Cardiac fibrosis is regarded as the characteristic pathological end-stage of cardiovascular disease and remains the most common chronic disease worldwide [1]. Cardiac fibrosis leads to high morbidity and mortality, as well as a substantially increased global medical burden [2]. Defined as the accumulation of extracellular matrix (ECM) proteins, cardiac fibrosis results from an imbalance in its production and degradation, thus contributing to cardiac dysfunction [3,4]. However, the specific mechanism of cardiac fibrosis is not clear, and its biomedical effects have never been satisfactorily determined [5,6]. Therefore, efforts should be made to promote the translation of mechanistic knowledge on cardiac fibrosis into potential clinical therapeutic applications.
As indispensable non-coding RNAs (ncRNAs), both long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) play important roles in the regulation of messenger RNA (mRNA) expression levels [7,8,9]. Regulation of the lncRNA–miRNA–mRNA–competitive endogenous RNA (ceRNA) network is actively involved in the onset and progression of multiple diseases, including cancers [8,10,11]. An increasing number of studies have demonstrated that massive mRNAs and ncRNAs (lncRNAs and miRNAs) could serve as crucial biomarkers for cardiac fibrosis [12,13,14]. Various ceRNA networks have been reported to mediate the pathological process of cardiac fibrosis [14,15]. Gao et al. explored a cardiac fibroblast-related ceRNA network, which provided potential target genes in this field [12]. A pro-fibrotic lncRNA (PFL, also known as NONMMUT022555) has been reported to act as a ceRNA of the cardioprotective miRNA, let-7d [16]. The ceRNA network of NORAD/miR-125a-3p/Fyn has also been shown to improve cardiac fibrosis and reduce inflammatory responses [17]. LncRNA H19 has been reported to act as a ceRNA to mediate the connective tissue growth factor (CTGF) expression by sponging miR-455 in cardiac fibrosis [18], whereas another study found that H19 alleviates cardiac fibrosis by targeting miR-22-3p/lysine (K)-specific demethylase 3A (KDM3A) in myocardial infarction [19]. Dai et al. revealed that lncRNA nuclear-enriched abundant transcript 1 (NEAT 1) regulates atrial fibrosis via the miR-320/neuronal per Arnt-Sim domain protein 2 (NPAS2) axis in atrial fibrillation [20]. In addition, lncRNA taurine upregulation gene 1 (TUG1) mediates CTGF expression by sponging miR-133b in myocardial fibrosis after myocardial infarction [21]. However, there is limited knowledge from the perspective of RNA sequencing regarding the roles of the ceRNA regulatory network in the pathogenesis of cardiac fibrosis, which requires further clarification.
In this study, we constructed a rat model of cardiac fibrosis induced using the chronic intermittent hypoxia (CIH) method. Whole transcriptome sequencing was performed to identify differentially expressed RNAs (DERs) and the related functional pathways. Protein–protein interaction (PPI) and ceRNA networks were constructed by integrating the whole transcriptome analysis. The present study aimed to reveal the crucial regulators and functional pathways related to cardiac fibrosis in rats and to further elucidate the underlying mechanisms of cardiac fibrosis.
2.
Materials and methods
2.1. Experimental animals and model of cardiac fibrosis
All the animal experiments were conducted on Sprague–Dawley rats (adult, male) using an experimental model of cardiac fibrosis induced by the CIH method [22]. The experimental rats were obtained and maintained as previously described [23]. The rat model of CIH adopted in this study was described in our previous study, where the establishment of cardiac fibrosis was successfully verified [22,23]. All the animal experiments and study protocols were approved by the hospital's Medical Ethics Committee.
2.2. Whole transcriptomics analysis
After successful establishment of the cardiac fibrosis model, the right atrial tissue samples of the control and CIH groups (n = 3/group) were collected and subjected to whole transcriptome sequencing as previously described [23]. Preprocessing, quality control, and normalization of the raw data was also described in our previous study [23]. Subsequently, the whole transcriptome data of lncRNAs, miRNAs, and mRNA were used for bioinformatics analyses of differential expression, functional and pathway enrichments, and construction of PPI and ceRNA regulatory networks in R.
2.3. DERs selection
The expression levels of lncRNA, miRNA, and mRNA between the control and CIH groups were compared to screen the DERs, including differentially expressed mRNAs (DEmRNAs), differentially expressed lncRNAs (DElncRNAs), and differentially expressed miRNAs (DEmiRNAs), by using the "limma" package in R3.6.1 [24]. |log2(fold change)| > 1 and P < 0.05 were selected as the screening thresholds. All three types of DERs are shown in volcano plots. Moreover, the expression of DERs was presented using bidirectional hierarchical clustering heatmaps based on the Pearson correlation algorithm [25,26].
2.4. Functional and pathway enrichment analysis
The screened DEmRNAs were then subjected to functional and pathway enrichment analyses using the Database for Annotation, Visualization, and Integrated Discovery (DAVID), version 6.8 [27]. Biological processes of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted [28]. A false discovery rate (FDR) value < 0.05 was considered as the threshold for enrichment significance.
2.5. Construction of PPI network
A PPI network was built to obtain the interactions and relationships among proteins encoded by DEmRNAs using the STRING version 11.0 database [29]. Interacting pairs were retained when the interaction scores were higher than 0.4. The PPI network was visualized using Cytoscape version 3.6.1 [30]. Next, functional and pathway enrichment analyses were performed for the nodes in the PPI network. The set threshold of enrichment significance was P < 0.05.
2.6. Construction of ceRNA regulatory network
The interactions of these DERs were predicted to construct a ceRNA regulatory network. In brief, the DElncRNA–DEmiRNA interactions were predicted using miRanda software. The parameters were set as Gap Extend = 0, Score Threshold = 80, Energy Threshold = −20, and Matched Seq% Threshold = 80%. Only connections with opposite directions between the DElncRNAs and DEmiRNAs were retained.
The target genes of DEmiRNAs in the DElncRNA–DEmRNA network were predicted using the miRWalk 3.0 database [31]. The DEmRNAs among the target genes were filtered. Only connections with opposite directions between the DEmiRNAs and DEmRNAs were retained.
The regulatory relationships in lncRNA–miRNA and miRNA–mRNA interactions were integrated to construct the ceRNA regulatory network, which was visualized using Cytoscape software Version 3.6.1 [30]. In addition, functional and pathway enrichment analyses were conducted based on the DEmRNAs in the ceRNA regulatory network. The threshold for enrichment significance was set at P < 0.05.
2.7. Identification of the crucial regulators in the related ceRNA regulatory network
According to the DEmRNAs involved in the KEGG pathways of the ceRNA regulatory network and miRNA–mRNA interactions of the ceRNA regulatory network, a Sankey diagram was generated to illustrate the regulatory relationships and relevant pathways. Moreover, the Comparative Toxicogenomics Database (CTD) was used to search for disease-relevant KEGG pathways and genes [32], and "atrial fibrillation" was used as the keyword. By comparing the screened genes and the important KEGG pathways involved in them, which were obtained from the previously constructed ceRNA regulatory network, the overlapping crucial regulators and relevant KEGG pathways were identified.
2.8. Quantitative reverse transcription PCR (qRT-PCR)
Total RNA was extracted from atrial specimens of rats using TRIzol reagent (Promega, Beijing, China), which was then reverse-transcribed routinely to cDNA (Transgen Biotech, Beijing, China). Next, qRT-PCR analysis was conducted using real-time PCR software with the SYBR Green PCR kit (Transgen Biotech, Beijing, China). The quantification of mRNA was measured using the 2-ΔΔCT method, and the level of lncRNA/genes and miRNA was separately normalized to GAPDH and U6. The primer sequences used in this study are listed in Table 1.
2.9. Statistical analysis
R software (R version 3.6.1) and Statistical Product and Service Solutions (SPSS 21.0) were used to conduct statistical analysis. P < 0.05 or FDR value < 0.05 was considered as a significant difference.
3.
Results
3.1. Identification of DERs and functional enrichment analysis
A total of 25,332 lncRNAs, 1204 miRNAs, and 22,601 mRNAs were identified. After the routine analysis for comparison and screening, 268 lncRNAs, 20 miRNAs, and 436 mRNAs were identified as DERs. All DERs are displayed in volcano plots (Figure 1A–C) and clustering heatmaps (Figure 1D–F). As shown in Figure 1D–F, the expression patterns of DERs were distinctly separated based on the different sample groups, which indicated that the identified DERs were expressed characteristically and specifically in cardiac fibrosis.
Next, DEmRNAs were subjected to functional and pathway enrichment analyses. As shown in Figure 2, a total of 18 relevant biological processes such as "chromosome segregation, " "cell division, " and "microtubule-based movement" were significantly enriched, and 6 KEGG signaling pathways such as "cell cycle, " "alcoholism, " and "viral carcinogenesis" were significantly enriched.
3.2. Construction of PPI network
The interactions and relationships among proteins encoded by DEmRNAs were predicted using the STRING database to construct a PPI network. A total of 3936 interacting pairs were obtained when interaction scores were higher than 0.4. A PPI network with 306 nodes was established (Figure 3A). Accordingly, the topological structural properties of all the 306 nodes in the PPI network were analyzed. The top 20 nodes are listed in Table 2 in the order of their degrees (from high to low), such as Cdk1, Ccnb1, Cdca8, and Aurkb. Functional and pathway enrichment analyses were also performed. As presented in Figure 3B, C, the relevant biological processes, such as "chromosome segregation, " "cell division, " and "microtubule-based movement, " were significantly enriched; in addition, KEGG signaling pathways, such as "pathways in cancer, " "alcoholism, " and "oocyte meiosis" were significantly enriched.
3.3. Construction of a ceRNA regulatory network
The communication between DERs was predicted by constructing a ceRNA regulatory network. Consequently, 262 DElncRNA–DEmRNA pairs and 125 DEmiRNA–DEmRNA pairs were obtained. By integrating these pairs, two ceRNA regulatory networks were established (Figure 4A, B). lncRNA–miRNA–mRNA and ceRNA interactions, such as NONRATT012985.2-rno–miR-3577–Arnt2, NONRATT000861.2-rno–miR-3577–BIRC5, and NONRATT019720.2-rno–miR-3577–GDF6, were identified from these networks.
Functional and pathway enrichment analyses were also performed. The results indicated that the relevant biological processes, such as "chromosome segregation, " "cell division, " and "spermatogenesis" were significantly enriched; in addition, KEGG signaling pathways, such as "microRNAs in cancer, " "pathways in cancer, " and "metabolic pathways" were significantly enriched (Figure 4C, D).
3.4. Identification of the crucial regulators in the related ceRNA regulatory network
According to the DEmRNAs involved in the KEGG pathways and miRNA–mRNA interactions, a Sankey diagram was drawn (Figure 5). Regulatory relationships and relevant pathways, such as miR-3577–Arnt2-pathways in cancer and the miR-3577–BIRC5–Hippo signaling pathway, were identified.
By searching the CTD database for "atrial fibrillation, " 154 KEGG pathways and 168 relevant genes were identified. After comparing the genes and KEGG pathways in the ceRNA regulatory network, 8 overlapping disease pathways and 1 directly involved gene, Arnt2, were identified. Arnt2 is involved in the disease pathway rno05200: Pathways in cancer. We speculated that the genes involved in the eight disease pathways might be important regulators of cardiac fibrosis. The relevant regulatory lncRNAs and miRNAs, possibly upstream of these genes, might also be closely associated with cardiac fibrosis. In particular, other genes, such as WNT2B, GNG7, LOC100909750, Cyp1a1, E2F1, BIRC5, and LPAR4, which are also involved in rno0520: pathways in cancer along with Arnt2, might be closely related to the pathogenesis of cardiac fibrosis. Finally, three mRNAs, Arnt2, BIRC5, and GDF6, one miRNA, miR-3577, and two lncRNAs, NONRATT011877.2 and NONRATT019720.2, were identified as crucial regulators in the ceRNA regulatory network.
3.5. Verification of the crucial regulators in the related ceRNA regulatory network by qRT-PCR
The above crucial regulators in the ceRNA regulatory network were selected for verification based on qRT-PCR experiments. As shown in Figure 6, the relative levels of mRNAs, Arnt2, BIRC5, and GDF6, were all remarkably upregulated in the CIH group compared with the control group (P < 0.05), the relative level of miRNA, miR-3577, was notably downregulated in the CIH group (P < 0.05), and the relative levels of lncRNAs, NONRATT011877.2 and NONRATT019720.2, were both significantly upregulated in the CIH group (P < 0.05). These results validated the consistency of the qRT-PCR experiments and RNA sequencing analyses.
4.
Discussion
Recent studies have shown that mRNAs and ncRNAs, including lncRNAs and miRNAs, play important roles in cardiac fibrosis [13]. Nevertheless, integrated whole-transcriptome analysis of cardiac fibrosis is still indispensable for the lncRNA–miRNA–mRNA-mediated ceRNA regulatory network [33].
In this study, we identified DERs, including 268 lncRNAs, 20 miRNAs, and 436 mRNAs, as well as the related functional pathways involved in cardiac fibrosis. Subsequently, a PPI network containing 306 nodes was constructed and the top 5 nodes were Cdk1, Ccnb1, Cdca8, Aurkb, and Plk1, all of which were upregulated in the myocardial tissues of cardiac fibrosis. A study in 2014 already proved that targeting Cdk1 may inhibit fibrosis and subsequently confer protection against cardiac fibrosis [34]. Chen et al. found that Cdk1 could promote atrial fibrosis by phosphorylating paxillin at Ser244 CDK1 and plays a key role in fibroblast differentiation [35]. However, other proteins of the top5 nodes have not yet been studied in cardiac fibrosis.
The results of the ceRNA regulatory network demonstrated specific lncRNA–miRNA–mRNA interactions in this study. For example, BIRC5 could be regulated by rno–miR-3577 as well as its upstream lncRNAs, including NONRATT012985.2 and NONRATT000861.2; Cyp1a1 could be regulated by rno–miR-673-5p and its upstream lncRNAs, including NONRATT006306.2 and TCONS_00006085; and E2F1 could be regulated by rno–miR-184 and its upstream lncRNAs, including NONRATT008981.2 and NONRATT019720.2. For these crucial regulators in the present study, previous studies have revealed that BIRC5 is expressed in cardiac progenitor cells (CPCs), and survivin (encoded by BIRC5) can directly induce CPCs proliferation and enhance cardiomyocyte survival. Transgenic Cyp1a1 is closely associated with hypertension and cardiac fibrosis. Liao et al. found that E2F1 may be a potential therapeutic target for cardiac fibrosis [36,37,38,39]. However, the underlying mechanisms have not been elucidated. The ceRNA relationships identified in this study may provide hints for further mechanistic investigations.
Meanwhile, the enriched functional pathways including "cell division" and "cell cycle" have also been reported by international studies as potential pathological factors of cardiac fibrosis [40,41,42]. For instance, the overlapping disease pathways in this study, including "pathways in cancer" and "Hippo signaling pathway, " have also been previously reported to be related with cardiac regeneration and proliferation, arbitrating the cardiac fibroblast identity and activation [43,44,45,46,47]. Finally, the identified crucial regulators in the related ceRNA regulatory network, three mRNAs, Arnt2, BIRC5, and GDF6, one miRNA, miR-3577, and two lncRNAs—NONRATT011877.2 and NONRATT019720.2—were screened out, and their differential expression was successfully validated by the qRT-PCR experiment. Nevertheless, only BIRC5 plays a role in the biological regulation of cardiac fibrosis. Therefore, this evidence could be regarded as the experimental foundation for these potential regulators of the ceRNA regulatory network of cardiac fibrosis.
However, there are still some limitations in this study, such as the lack of in vitro data, depletion of specific genes, and the need of further experiments for mechanism exploration and pathway validation. More importantly, how the ceRNA regulatory network impacts the advancement of knowledge in clinical patients subjected to cardiac fibrosis, how to prove medically the effectiveness of these crucial regulators and the functional pathways, and how to practice the significance of clinical transformation would need more investigations and efforts in future medical practice.
Compared with the results of a previous study, this work provides a comprehensive understanding of the ceRNA network-based regulation of cardiac fibrosis. In this study, different crucial regulators and functional pathways were identified as closely related to cardiac fibrosis in a newly established ceRNA regulatory network based on an experimental rat model. In conclusion, the present study mainly indicates the integration of whole-transcriptome data to construct the ceRNA network and identify the crucial regulators, including Arnt2, BIRC5, and GDF6, to construct networks such as NONRATT012985.2-rno–miR-3577–Arnt2, NONRATT000861.2-rno–miR-3577–BIRC5, and NONRATT019720.2-rno–miR-3577–GDF6, and trace the related functional pathways in cardiac fibrosis in rats, thus suggesting the potential underlying mechanisms of cardiac fibrosis.
Acknowledgments
This work was supported by the Science & Technology Development Fund of Tianjin Education Commission for Higher Education [No. 2020KJ168].
Conflict of interest
The authors declare that they have no conflicts of interest.