M6A level of circRNAs was reduced in HPH rats and most circRNAs contained one m6A peak
3 weeks treatment by hypoxia resulted in right ventricular systolic pressure (RVSP) elevating to 42.23 ± 1.96 mmHg compared with 27.73 ± 1.71 mmHg in the control (P < 0.001, Figure 1A and 1B). The ratio of the right ventricle (RV), left ventricular plus ventricular septum (LV + S) [RV/ (LV + S)] was used as an index of RVH. RVH was indicated by the increase of RV/ (LV + S) compared with the control (0.25 ± 0.03 vs. 0.44 ± 0.04, P = 0.001, Figure 1C). The medial wall of the pulmonary small arteries was also significantly thickened (19.28 ± 2.19% vs. 39.26 ± 5.83%, P < 0.001, Figure 1D and 1E). Moreover, in the normoxia group, 53.82 ± 3.27% of the arterioles were non-muscularized (NM) vessels, and 25.13 ± 1.83% were fully muscularized (FM) vessels. In contrast, partially muscularized vessels (PM) and FM vessels showed a greater proportion (32.88 ± 3.15% and 41.41 ± 3.35%) in HPH rats, while NM vessels occupied a lower proportion (25.71 ± 2.55%) (Figure 1F). Figure 1G displayed the heatmap of m6A circRNAs expression profiling in N and HPH. M6A abundance in 166 circRNAs was significantly upregulated. Meanwhile, m6A abundance in 191 circRNAs was significantly downregulated (Additional file 1: Data S1, filtered by fold change ³ 4 and P £ 0.00001). Lungs of N and HPH rats were selected to measure m6A abundance in purified circRNAs. The m6A level in total circRNAs isolated from lungs of HPH rats was lower than that from controls (Figure 1H). Moreover, over 50% circRNAs contained only one m6A peak either in lungs of N or HPH rats (Figure 1I).
M6A circRNAs were mainly from protein-coding genes spanned single exons in N and HPH groups
We analyzed the distribution of the parent genes of total circRNAs, m6A-circRNAs, and non-m6A circRNAs in N and HPH, respectively. N and HPH groups showed a similar genomic distribution of m6A circRNAs and non-m6A circRNAs (Figure 2A and 2B). Moreover, about 80% of m6A circRNAs and non-m6A circRNAs were derived from protein-coding genes in both groups. A previous report indicated that most circRNAs originated from protein-coding genes spanned two or three exons [14]. While in our study, over 50% and 40% of total circRNAs from protein-coding genes spanned one exon in N and HPH groups, respectively (Figure 2C and 2D). Similarly, m6A circRNAs and non-m6A circRNAs were mostly encoded by single exons. Therefore, it was indicated that m6A methylation was abundant in circRNAs originated from single exons in N and HPH groups.
The distribution and functional analysis for host genes of circRNAs with differentially expressed m6A peaks
The length of differentially-expressed m6A circRNAs was mostly enriched in 1-10000 bps (Figure 3A). The host genes of upregulated m6A circRNAs were located in chromosome 1, 2 and 10, while the downregulated parts were mostly located in chromosome 1, 2 and 14 (Figure 3B).
Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to explore the host genes of circRNAs with differentially-expressed m6A peaks. In the GO analysis (Figure 3C, left), the parent genes of circRNAs with upregulated m6A peaks were enriched in the protein modification by small protein conjugation or removal and macromolecule modification process in the biological process (BP). Organelle and membrane-bounded organelle were also the two largest parts in the cellular component (CC) analysis. Binding and ion binding were the two main molecular functions (MF) analysis. The top 10 pathways from KEGG pathway analysis were selected in the bubble chart (Figure 3C, right). Among them, the oxytocin signaling pathway, protein processing in endoplasmic reticulum and cGMP-PKG signaling pathway were the top 3 pathways involved. In addition, vascular smooth muscle contraction pathway was the most associated pathway in PH progression [27].
In Figure 3D left, the parent genes of circRNAs with downregulated m6A peaks were mainly enriched in the cellular protein modification process and protein modification process in BP. Organelle and membrane-bounded organelle made up the largest proportion in the CC classification. The MF analysis was focused on receptor signaling protein activity and protein binding. The parent genes of circRNAs with decreased m6A peaks were mainly involved in the tight junction and lysine degradation in the KEGG pathway analysis (Figure 3D, right).
Hypoxia can influence the m6A level of circRNAs and circRNAs abundance
360 m6A circRNAs were shared in N and HPH groups. 49% of circRNAs modified by m6A were only detected in N group, and 54% of circRNAs modified by m6A were only detected in HPH group (Figure 4A). To explore whether m6A methylation would influence circRNAs expression level, expression of the 360 common m6A circRNAs were identified. More circRNAs tended to decrease in HPH compared to N (Figure 4B). Moreover, expression of m6A circRNAs was significantly downregulated compared with non-m6A circRNAs in hypoxia, suggesting that m6A may downregulate the expression of circRNAs in hypoxia (Figure 4C, P = 0.0465).
Construction of a circRNA–miRNA–mRNA co-expression network in HPH
We found 76 upregulated circRNAs with increased m6A abundance, and 107 downregulated circRNAs with decreased m6A abundance (Figure 5A, Additional file 2: Data S2). As known, circRNAs were mostly regarded as a sponge for miRNAs and regulated the expression of corresponding target genes of miRNAs [28]. To explore whether circRNAs with differentially-expressed m6A abundance influence the availability of miRNAs to target genes, we selected differentially-expressed circRNAs with increased or decreased m6A abundance. GO enrichment analysis and KEGG pathway analysis were also performed to analyze target mRNAs. Target mRNAs displayed similar GO enrichment in the two groups (Figure 5B and 5C). Two main functions were determined in BP analysis: positive regulation of biological process and localization. Intracellular and intracellular parts make up the largest proportion in CC part. Target mRNAs were mostly involved in protein binding and binding in MF part. In the KEGG pathway analysis, the top 10 most enriched pathways were selected (Figure 5D and 5E). Wnt and FoxO signaling pathways were reported to be involved in PH progression [29-31]. Then, we analyzed the target genes involved in these two pathways. SMAD4 was associated with PH and involved in Wnt signaling pathways. MAPK3, SMAD4, TGFBR1, and CDKN1B were involved in FoxO signaling pathways. To explore the influence of circRNA-miRNA regulation on PH-associated genes expression, we constructed a circRNA-miRNA-mRNA network, integrating matched expression profiles of circRNAs, miRNAs and mRNAs (Figure 5F and 5G). MicroRNAs sponged by the target genes of interest were analyzed. MiR-125a-3p, miR-23a-5p, miR-98-5p, let-7b-5p, let-7a-5p, let-7g-5p, and miR-205 were analyzed because they were reported to be associated with PH [32, 33]. We filtered the key mRNAs and miRNAs, and founded that the two circRNAs were the most enriched, which were originated from chr1:204520403-204533534- (Xpo6) and chr7:40223440-40237400- (Tmtc3).
M6A circXpo6 and m6A circTmtc3 were downregulated in PASMCs and PAECs in hypoxia
M6A abundance was significantly reduced in PASMCs and PAECs when exposed to hypoxia (0.107% ± 0.007 vs. 0.054% ± 0.118, P = 0.023 in PASMCs; 0.114% ± 0.011 vs. 0.059% ± 0.008, P = 0.031 in PAECs, Figure 6A). M6A abundance in circRNAs was lower than it in mRNAs (0.1–0.4%) [17, 18]. Next, we confirmed the back-splicing of circXpo6 and circTmtc3 by CIRI software. The sequence of linear Xpo6 and Tmtc3 mRNA was analyzed. Then we identified that circXpo6 was spliced form exon 7, 8, and 9 of Xpo6. CircTmtc3 was spliced form exon 8, 9, 10, and 11 (Figure 6B). Using cDNA and genomic DNA (gDNA) from PASMCs and PAECs as templates, circXpo6 and circTmtc3 were only amplified by divergent primers in cDNA, while no product was detected in gDNA (Figure 6C). To identify whether circXpo6 and circTmtc3 were modified by m6A, we performed M6A RNA Immunoprecipitation (MeRIP)-RT-PCR and MeRIP-quantitative RT-PCR (MeRIP-qRT-PCR) to detect the expression of circXpo6 and circTmtc3 (Figure 6D and 6E). m6A circXpo6 and m6A circTmtc3 were significantly decreased in PASMCs and PAECs when exposed to hypoxia (P = 0.002, and P = 0.015 in PASMCs and P = 0.02, and P = 0.047 in PAECs)