RNA interactions in right ventricular dysfunction induced type II cardiorenal syndrome

PAC induced RV dysfunction

Compared with the sham groups four weeks after surgery, the PAC group showed an enlargement of the right ventricle (Figure 1A, 1B). Masson staining revealed an increase in cardiac fibrosis in both the perivascular and intermuscular areas (Figure 1C), and expression of the fibrotic marker collagen I increased significantly (Figure 1D). Similarly, heart weight (Figure 1E) and right heart weight (Figure 1F) normalized to body weight increased significantly after PAC surgery.

Echocardiography revealed higher PV peak velocity (Figure 2A, 2D), higher PV max pressure (Figure 2E), thicker RV free wall (Figure 2B, 2F), and larger RV internal diameter (Figure 2C, 2G) in the PAC group than in the corresponding sham group. The significant differences in the tricuspid valve E/A ratio (TV E/A) and Tei index indicated RV diastolic dysfunction in the PAC group (Figure 2H, 2I). Moreover, a significantly lower tricuspid annular plane systolic excursion (TAPSE) in the PAC group indicated the existence of systolic dysfunction (Figure 2J). Taken together, PAC successfully led to RV dysfunction.

RV dysfunction induced renal dysfunction

In addition to RV dysfunction, we found that PAC mice showed signs of renal dysfunction. PAC mice displayed obvious renal atrophy with statistical significance (Figure 3A, 3B). Renal fibrosis in the PAC group was detected by Masson staining with an increase of approximately 15% of the fibrotic area compared to the sham group (Figure 3C). The fibrotic markers collagen I and III increased simultaneously (Figure 3D, 3E). We also examined renal function, noting a significant increase in plasma creatinine in the PAC group compared with the sham group (Figure 3F). Moreover, the level of neutrophil gelatinase-associated lipocalin (NGAL), an indicator of kidney injury [23], was significantly higher in the PAC group (Figure 3G). Thus, we generated an RV dysfunction-induced type II CRS mouse model.

Changes in coding and noncoding RNA expression patterns in response to PAC

Research to date has rarely addressed coding and noncoding RNA expression profiles in RV dysfunction-induced type II CRS models. To obtain a comprehensive understanding of transcriptome alterations in the heart and kidney, we performed whole transcriptome RNA sequencing on these tissues in sham and PAC mice.

A total of 741 DEmRNAs (556 upregulated and 185 downregulated) were detected in RV tissue and 86 DEmRNAs (49 upregulated and 37 downregulated) in kidney tissue according to the threshold of FDR < 0.05 and | log2(Fold Change) | ≥ 1. The top ten up- and downregulated DEmRNAs in the RV and kidney are shown in Tables 1, 2, respectively. Hierarchically clustered heat maps and volcano plots illustrated DEmRNAs in different tissues between the PAC and sham groups (Supplementary Figure 1A, 1B).

In addition, 159 DEmiRNAs were found in the RV, of which 123 were upregulated and 36 downregulated. In the kidney, there were 29 DEmiRNAs including 25 that were upregulated and four that were downregulated. The top ten up- and downregulated DEmiRNAs in the RV are shown in Table 3, and the top ten up- and six downregulated DEmiRNAs in the kidney are shown in Table 4. Visualization DEmiRNAs in different tissues between PAC and sham groups using hierarchical clustering heat maps and volcano plots is shown in Supplementary Figure 1C, 1D.

Based on mapping reads from filtered ribosomes to the reference genome, 13,648 circRNAs were identified, with 467 existing and 13,182 novel ones. In the RV, there were 130 upregulated circRNAs and 102 downregulated circRNAs, of which 59.27% were exonic, 16.33% were antisense, 7.35% were intronic and 1.6% were intergenic (Supplementary Figure 2A, 2B). A total of 41 and 62 DEcircRNAs in the kidney were detected, with exonic, antisense, intronic and intergenic circRNAs accounting for 66.44%, 10.18%, 6.28% and 1.65%, respectively (Supplementary Figure 2C, 2D). The top ten up- and downregulated DEcircRNAs in RV and kidney are provided in Tables 5, 6, respectively. These circRNAs are widely distributed among all chromosomes, including the mitochondrial chromosome and the sex chromosomes. Among all circRNAs identified, the top three chromosomes of circRNA distribution were found to be chromosomes 2 (chr 2), chr 1 and chr 11 in the sham and PAC groups for both the RV and kidney (Supplementary Figure 2E). The identified circRNAs are mainly between 201 to 700 bases in length (Supplementary Figure 2F).

Validation of RNA sequencing data by RT-PCR

To validate the sequencing results, we randomly selected five mRNAs, three miRNAs and two circRNAs in the two tissues for qRT-PCR. Aqp1 was downregulated in both tissues. By contrast, Acta1 was significantly upregulated in the RV but remained the same in the kidney. Although expression of Scp2 in the RV was not significantly different between the sham and PAC groups, it was significantly increased in the kidney. Thoc2 was only increased in the RV between the PAC and sham groups, but there was no significant change in the kidney. Supporting the existence of fibrosis in the RV and kidney, Tgfβ3 expression in both tissues was significantly higher in PAC mice than in sham mice (Figure 4A).

Regarding miRNA expression levels, miR-205-5p and miR-208b-3p were significantly upregulated in the RV tissue in the PAC group compared with the sham group; miR-21a-3p displayed an increasing trend with no significant difference between the two groups. MiR-202-5p and miR-122-3p were significantly upregulated in the kidneys of PAC mice compared with sham mice, whereas miR-215-5p was significantly downregulated (Figure 4B).

The qRT-PCR results showed mm9_circ_004072 and mm9_circ_007155 to be significantly upregulated in the RV of PAC mice. Mm9_circ_010239 showed a downregulation of approximately 3 times in the kidney, though expression of mm9_circ_013197 remained the same (Figure 4C).

All of the findings were consistent with the sequencing results, potently validating the reliability of our RNA sequencing.

Functional prediction of DEmiRNAs in RV tissue

To understand the potential molecular mechanisms of RV dysfunction, we first constructed pairs of miRNAs and their target genes with a negative correlation using a threshold of P < 0.05 and | log2(Fold Change) | ≥ 1. Functional analysis of the negatively regulated pairs was conducted according to the regulation direction of the target mRNAs (upregulated mRNAs with downregulated miRNAs group; downregulated mRNAs with upregulated miRNAs group).

A total of 2227 pairs containing 43 miRNAs and 210 target mRNAs were found in RV tissue. GO analysis of these upregulated mRNAs in the RV of PAC mice revealed significant enrichment in the biological process of positive regulation of the cell cycle and Ras protein signal transduction (Supplementary Figure 3A); and the downregulated mRNAs in RV showed enrichment mainly in metabolic processes, including low-density lipoprotein receptor particle metabolic process, carboxylic acid biosynthetic process and organic acid biosynthetic process (Supplementary Figure 3B). Based on KEGG analysis, the upregulated mRNAs in PAC mice were mostly enriched in the PI3K-Akt signaling pathway and focal adhesion pathway (Supplementary Figure 3C), whereas the downregulated mRNAs were mostly enriched in the PPAR signaling pathway (Supplementary Figure 3D).

Functional prediction in renal tissue

Renal dysfunction was observed in the PAC group, in order to uncover genetic alterations related to kidney dysfunction, GO and KEGG analyses were performed. As depicted in Supplementary Figure 1D, we noticed an abnormal opposite expression pattern of the first sample in the PAC group. To eliminate bias, the data of this sample were not included in the analysis.

A total of 140 miRNA-mRNA pairs consisting of 46 DEmRNAs and 31 DEmiRNAs were involved in the analysis. GO analysis of the upregulated mRNAs in renal tissue of the PAC group revealed significant biological process enrichment in positive regulation of cell migration, small GTPase-mediated signal transduction and positive regulation of neuron differentiation (Supplementary Figure 4A). and the downregulated mRNAs in the kidneys were also enriched mainly in metabolic processes including cofactor, coenzyme and lipid metabolic processes (Supplementary Figure 4B). KEGG analysis of upregulated mRNAs indicated enrichment in the PI3K-Akt signaling pathway, which was interestingly the same as for RV tissue (Supplementary Figure 4C), and the downregulated mRNAs were enriched in the phospholipase D signaling pathway and apelin signaling pathway (Supplementary Figure 4D).

Construction of regulation networks of miRNA-mRNA

To demonstrate the dysregulated relationship of miRNA-mRNA pairs, we first built miRNA-mRNA regulation network based on the expression levels of DEmRNAs and DEmiRNAs. In the RV tissue sample, to highlight the contribution of hub miRNA-mRNA pairs to the fibrotic and metabolic regulation, we selected DE pairs with the threshold of the abstract value of log₂(Fold Change) ≥ 10, with 34 DEmRNAs and 132 DEmiRNAs consisting of 333 pairs involved in the network (Supplementary Figure 5A). In the renal tissue sample, 46 DEmRNAs and 31 DEmiRNAs were involved in the network. The top three upregulated DEmiRNAs were miR-33-y, miR-7667-3p and miR-7668-5p. The top three downregulated DEmiRNAs were miR-6986-5p, miR-541-3p and miR-671-x (Supplementary Figure 5B).

Construction of ceRNA networks of circRNA-miRNA-mRNA

Then, based on the circRNA-miRNA and miRNA-mRNA targeting relationship, we predicted ceRNA pairs of circRNA-miRNA-mRNA. Interestingly, functional analysis of the ceRNA network showed similar to those of miRNA-mRNA: the most significantly changed biological processes were fibrosis, cell proliferation, and metabolism-related processes.

A total of 1104 upregulated ceRNA pairs were found in the RV tissue, including 32 circRNAs, 263 miRNAs, and 66 mRNAs. While 1,236 ceRNA pairs were downregulated, including 32 circRNAs, 295 miRNAs, and 58 target mRNAs. The network of the paired downregulated ceRNA is shown in Figures 5A, 6A. In right heart tissue, cardiac hypertrophy and metabolic processes were both dysregulated. The significantly upregulated ceRNA pairs were mainly enriched in the actin cytoskeleton regulatory pathway and the cGMP-PKG signaling pathway (Figure 5B). Next, according to the number of edges connected with each node, we selected hub RNAs located in key positions in the network (Figure 5C). In the upregulated ceRNA network for the RV, the RNAs with the highest number of edges were Rasal2, Ubn2 and novel_circ_002631. As illustrated in Figure 6B, the significantly downregulated ceRNA pairs were mainly enriched in metabolic pathways, including the thyroid hormone signaling pathway, adrenergic signaling pathway, valine, leucine and isoleucine degradation. Hub RNAs were Asxl2, Tmem245 and Strbp (Figure 6C).

In kidney tissue, 608 upregulated ceRNA pairs were selected, including 32 circRNAs, 175 miRNAs, and 20 mRNAs (Figure 7A). GO analysis suggested metabolism to be the most significantly dysregulated biological process, whereas KEGG analysis showed the significantly upregulated pairs to be enriched in pathways such as vitamin metabolism, sphingolipid metabolism and lysine degradation (Figure 7B). Hub RNAs were Akap7, Ubn2 and novel_circ_002631 (Figure 7C). A total of 834 ceRNA pairs were downregulated, as shown in Figure 8A, consisting of 32 circRNAs, 230 miRNAs and 41 mRNAs. Downregulated pairs were most enriched in the thyroid hormone signaling pathway and the FoxO signaling pathway (Figure 8B), with Elf1, Ercc6 and N4bp2l2 as hub RNAs (Figure 8C).

Study limitations

Few studies have been carried out to elucidate the molecular mechanisms of RV dysfunction-induced type II CRS. To achieve a comprehensive understanding of potential mechanisms and therapeutic targets in type II CRS, we focused on bioinformatics analysis in this study. Detailed research on several selected targets will be carried out in our future studies. However, there are several limitations of this study that should be acknowledged. First, the addition of a treatment group to improve type II CRS would provide further clues to determine what signaling pathways play more important roles in the pathogenesis of this disease model. To evaluate expression profiling changes due to intervention, it might be helpful to establish a targeted genetic approach for future studies. Second, we only performed a single time-point analysis, but whether dysfunction of the heart or kidney develops with time and the molecular signatures may vary at different time points. Further study at different time points after PAC should be conducted to better clarify the dynamic changes in circRNA-miRNA-mRNA interactions. Third, our results were obtained in a mouse model, and further exploration should be performed to determine whether they can be translated to humans. Accumulating evidence suggests that certain noncoding RNA expression patterns are conserved across species [46], and further human-specific studies are necessary to better understand the molecular mechanisms.

PAC model

Male C57BL/6 mice aged 7-8 weeks and weighing 22-25 grams were anesthetized by intraperitoneal injection of a mixture of xylazine (5 mg.kg⁻¹) and ketamine (100 mg.kg⁻¹), intubated with PE-50 tubing, and ventilated with room air using a mouse mini-ventilator. The respiration rate was set between 100 and 110 times per minute. PAC was performed as described previously [47]. After anesthesia, an incision was made in the skin parallel to the second rib, approximately 10 mm in length, with ophthalmic scissors, ensuring that the incision starts from the sternal angle and ends on the left anterior axillary line. Then, the second intercostal space was identified by counting the ribs from the sternal angle. The pectoralis major and pectoralis minor muscles were separated and cut above the second intercostal space to expose this space. The second intercostal space was bluntly penetrated and opened. The parenchyma and thymus were carefully separated until the pulmonary trunk was visible, and the PA and the ascending aorta were bluntly separated. A 6-0 braided silk suture was placed through the connective tissue between the PA and the ascending aorta, and after confirming that there was no bleeding in the PA and aorta, the PA together with a padding needle was ligated. The padding needle was removed immediately after the filling of the pulmonary conus was observed. After closing the chest and the skin, the skin was disinfected with 75% alcohol.

Echocardiography

Four weeks after surgery, echocardiography measurements (VEVO 2100, Visual Sonic, Toronto, ON, Canada) of the right heart were performed under anesthetic (1.5–2% isofluorane, 2 L/min oxygen flow rate) with mice intubated as reported elsewhere [48, 49]. Briefly, measurements of the RV outflow velocity were obtained from the parasternal short-axis view at the level of the pulmonary valve during end-diastole. After aligning the RV free wall perpendicularly to the transducer, the RV ejection time (ET), RV free wall thickness and RV internal diastolic diameter (RVIDd) were measured. The tricuspid closure-to-opening time (TCOT), tricuspid annular plane systolic excursion (TAPSE) and tricuspid valve E/A ratio (TV E/A) were measured from the apical four-chamber view. The RV myocardial performance index (or Tei index) was calculated as (TCOT-ET)/ET to evaluate RV function [50]. After data collection, the mice were sacrificed by an overdose of pentobarbital (150 mg/kg, i.p) and their hearts and kidneys were extracted for further analysis.

Histological examination and plasma creatinine measurement

Hearts and kidneys were fixed in 10% formalin and embedded in paraffin by standard protocols. Staining with hematoxylin and eosin (H&E) was performed to evaluate myocardial hypertrophy and Azan-Masson staining was used to assess the fibrosis area [48, 51]. Plasma creatinine (Cr) levels were assessed using a Mouse Cr ELISA Kit (Huabo Deyi Biotech., Beijing, China). ImageJ software was used for quantification.

qRT-PCR

Total RNA was extracted from mouse cardiac and kidney tissues with a total RNA isolation system (Omega, Norcross, GA, USA). One microliter of RNA was used to measure the quantity and quality with a NanoDrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA). Total RNA was converted to cDNA by using oligo (dT) primers with a PrimeScriptTM RT Master Mix (Takara Bio Inc., Shiga, Japan) or Mir-XTM miRNA First-Strand Synthesis Kit (Takara Bio USA Inc., CA, USA). Expression of genes was determined by Quantitect SYBR green real-time PCR (Takara Bio Inc., Shiga, Japan) using the primers listed in Table 7.

RNA purification, library construction and sequencing

After total RNA was extracted from mouse cardiac and renal tissues, RNA molecules in a size range of 18–30 nt were enriched by polyacrylamide gel electrophoresis (PAGE). The 3’ adapters were added, and 36-44 nt RNAs were enriched; 5’ adapters were then ligated to the RNAs. The ligation products were reverse transcribed by PCR amplification, and 140-160 bp PCR products were enriched to generate a cDNA library that was sequenced using Illumina HiSeqTM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China). Clean reads obtained were searched against the miRBase database (http://www.mirbase.org/. Release 21) to identify known miRNAs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses

MiRNAs and mRNAs with a fold change of ≥ 2 or ≤ -2 and false discovery rate (FDR)<0.05 were considered significantly differentially expressed, and chosen for further analysis. GO and KEGG pathway analyses were conducted to investigate the functions of all DEmRNAs and DEmiRNAs identified in this study. DERNAs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/) and significantly enriched GO terms were used to examine DERNAs from a biological process perspective. KEGG pathway analysis (http://www.genome.jp/kegg/) was performed to reveal pathway clusters related to the DERNAs. FDR was calculated to correct P values.

MiRNA target prediction and sponge analysis

RNAhybrid (v2.1.2) + svm_light (v6.01), Miranda (v3.3a) and TargetScan (Version:7.0) were used to predict miRNA targets. miRNA sequences and family information were obtained from TargetScan website (http://www.targetscan.org/). For circRNAs that have been annotated in circBase, the target relationship with miRNAs can be predicted by StarBase (v2.0). For novel circRNAs, three software programs Mireap, Miranda (v3.3a) and TargetScan (Version:7.0) were used to predict targets in animal samples.

Construction of the circRNA-miRNA-mRNA network

To identify interactions among DEcircRNAs, DEmiRNAs and DEmRNAs, we constructed a coexpression network based on a correlation analysis of the DERNAs. Correlation of expression between miRNA targets was evaluated using the Pearson correlation coefficient (PCC). Pairs with PCC<-0.7 and P<0.05 were selected as coexpressed negative miRNA-mRNA target pairs, and all RNAs were differentially expressed. For prediction of mRNAs interacting with circRNAs and miRNAs, miRTarBase (v6.1) was used to predict mRNAs targeted by miRNA sponges, and the circRNA-miRNA-mRNA coexpression network was visualized using Cytoscape software (v3.6.0) (http://www.cytoscape.org/).

Statistical analysis

Quantitative data are reported as the mean ± standard error. Statistical analysis for comparison between two groups was performed using two-tailed unpaired t-test. All analyses were performed with GraphPad Prism 7.0 software (GraphPad Software, Inc., CA, USA). A P value < 0.05 was deemed statistically significant.