Transcriptome-wide map of m6A circRNAs identified in a rat model of hypoxia mediated pulmonary hypertension

doi:10.21203/rs.2.15953/v4

PDF

Research article

Transcriptome-wide map of m6A circRNAs identified in a rat model of hypoxia mediated pulmonary hypertension

https://doi.org/10.21203/rs.2.15953/v4

This work is licensed under a CC BY 4.0 License

Journal Publication

published 13 Jan, 2020

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background: Hypoxia mediated pulmonary hypertension (HPH) is a lethal disease and lacks effective therapy. CircRNAs play significant roles in physiological process. Recently, circRNAs are found to be m 6 A-modified. The abundance of circRNAs was influenced by m 6 A. Furthermore, the significance of m 6 A circRNAs has not been elucidated in HPH yet. Here we aim to investigate the transcriptome-wide map of m 6 A circRNAs in HPH. Results: Differentially expressed m 6 A abundance was detected in lungs of HPH rats. M 6 A abundance in circRNAs was significantly reduced in hypoxia in vitro . M 6 A circRNAs were mainly from protein-coding genes spanned single exons in control and HPH groups. Moreover, m 6 A influenced the circRNA–miRNA–mRNA co-expression network in hypoxia. M 6 A circXpo6 and m 6 A circTmtc3 were firstly identified to be downregulated in HPH. Conclusion: Our study firstly identified the transcriptome-wide map of m 6 A circRNAs in HPH. M 6 A can influence circRNA–miRNA–mRNA network. Furthermore, we firstly identified two HPH-associated m 6 A circRNAs: circXpo6 and circTmtc3. However, the clinical significance of m 6 A circRNAs for HPH should be further validated. Key words: m 6 A circRNAs; hypoxia mediated pulmonary hypertension; m 6 A circXpo6; m 6 A circTmtc3

Epigenetics & Genomics

m6A circRNAs

hypoxia mediated pulmonary hypertension

m6A circXpo6

m6A circTmtc3

Pulmonary hypertension (PH) is a lethal disease and defined as an increase in the mean pulmonary arterial pressure ³ 25 mmHg at rest, as measured by right heart catheterization [1]. Hypoxia mediated pulmonary hypertension (HPH) belongs to group III PH according to the comprehensive clinical classification of PH, normally accompanied by severe chronic obstructive pulmonary disease (COPD) and interstitial lung diseases [2]. HPH is a progressive disease induced by chronic hypoxia [1]. Chronic hypoxia triggers over-proliferation of pulmonary artery endothelial cells (PAECs) and pulmonary artery smooth muscle cells (PASMCs), and activation of quiescent fibroblasts, the hallmark of HPH [1, 3]. The pathological characteristics of HPH are pulmonary vascular remolding, pulmonary hypertension, and right ventricular hypertrophy (RVH) [4]. So far there is no effective therapy for HPH [2]. More effective therapeutic targets are needed to be discovered.

Circular RNAs (circRNAs) were firstly found abundant in eukaryotes using RNA-seq approach [5-7]. Pre-mRNA is spliced with the 5’ and 3’ ends, forming a ‘head-to-tail’ splice junction, then circRNAs are occurred [5]. According to the genome origin, circRNAs may be classified into four different subtypes: exonic circRNA, intronic circRNA, exon–intron circRNA and tRNA introns circRNA [5]. CircRNAs are reported to play crucial roles in miRNA binding, protein binding, regulation of transcription, and post-transcription [5, 8]. Recent reports indicated that circRNAs can translate to proteins [8, 9]. Moreover, circRNAs are widely expressed in human umbilical venous endothelial cells when stimulated by hypoxia [10, 11]. Up to date, only a few reports mentioned PH-associated circRNAs. CircRNAs expression profile is demonstrated in HPH and chronic thromboembolic pulmonary hypertension [12]. However, the post-transcript modification of circRNAs in HPH is still unknown.

N⁶-methyladenosine (m⁶A) is regarded as one part of “epitranscriptomics’’ and identified as the most universal modification on mRNAs and noncoding RNAs (ncRNAs) in eukaryotes [13, 14]. DRm⁶ACH (D denotes A, U or G; R denotes A, G; H denotes A, C, or U) is a consensus motif occurred in m⁶A modified RNAs [15-17]. M⁶A modification is mainly enriched around the stop codons, at 3’untranslated regions and within internal long exons [17-19]. Several catalyzed molecules act as “writers’’, ‘‘readers’’, and ‘‘erasers’’ to regulate the m⁶A modification status [14]. The methyltransferase complex is known as writers, including methyltransferase-like-3, -14 and -16 (METTL3/METTL14/METTL16), Wilms tumor 1-associated protein (WTAP), RNA binding motif protein 15 (RBM15), vir like m⁶A methyltransferase associated (KIAA1429) and zinc finger CCCH-type containing 13 (ZC3H13), appending m⁶A on DRACH [17, 20, 21]. METTL3 is regarded as the core catalytically active subunit, while METTL14 and WTAP play a structural role in METTL3’s catalytic activity [18, 22]. The ‘‘erasers”, fat mass and obesity related protein (FTO) and alkylation repair homolog 5 (ALKBH5), catalyze the N-alkylated nucleic acid bases oxidatively demethylated [22]. The ‘‘readers’’, the YT521-B homology (YTH) domain-containing proteins family includes YTHDF (YTHDF1, YTHDF2, YTHDF3), YTHDC1, and YTHDC2, specifically recognizes m⁶A and regulates splicing, localization, degradation and translation of RNAs [14, 22, 23]. The YTHDF1 and YTHDF2 crystal structures forms an aromatic cage to recognize m⁶A sites in cytoplasm [24]. YTHDC1 is the nuclear reader and YTHDC2 binds m⁶A under specific circumstances or cell types [24]. Hypoxia may alter the balance of writers-erasers-readers and induce tumor growth, angiogenesis, and progression [25, 26].

Interestingly, circRNAs can be m⁶A-modified. M⁶A circRNAs displayed cell-type-specific methylation patterns in human embryonic stem cells and HeLa cells [14]. CircRNAs contained m⁶A modifications are likely to promote protein translation in a cap-independent pattern [9]. However, m⁶A circRNAs has not been elucidated in HPH yet. Here we are the first to identify the expression profiling of m⁶A circRNAs in HPH.

M⁶A level of circRNAs was reduced in HPH rats and most circRNAs contained one m⁶A 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 m⁶A circRNAs expression profiling in N and HPH. M⁶A abundance in 166 circRNAs was significantly upregulated. Meanwhile, m⁶A 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 m⁶A abundance in purified circRNAs. The m⁶A 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 m⁶A peak either in lungs of N or HPH rats (Figure 1I).

M⁶A 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, m⁶A-circRNAs, and non-m⁶A circRNAs in N and HPH, respectively. N and HPH groups showed a similar genomic distribution of m⁶A circRNAs and non-m⁶A circRNAs (Figure 2A and 2B). Moreover, about 80% of m⁶A circRNAs and non-m⁶A 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, m⁶A circRNAs and non-m⁶A circRNAs were mostly encoded by single exons. Therefore, it was indicated that m⁶A 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 m⁶A peaks

The length of differentially-expressed m⁶A circRNAs was mostly enriched in 1-10000 bps (Figure 3A). The host genes of upregulated m⁶A 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 m⁶A peaks. In the GO analysis (Figure 3C, left), the parent genes of circRNAs with upregulated m⁶A 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 m⁶A 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 m⁶A peaks were mainly involved in the tight junction and lysine degradation in the KEGG pathway analysis (Figure 3D, right).

Hypoxia can influence the m⁶A level of circRNAs and circRNAs abundance

360 m⁶A circRNAs were shared in N and HPH groups. 49% of m⁶A circRNAs detected in N group were not detected in HPH group, and 54% of m⁶A circRNAs detected in HPH group were not detected in N group (Figure 4A). To explore whether m⁶A methylation would influence circRNAs expression level, expression of the 360 common m⁶A circRNAs were identified. More circRNAs tended to decrease in HPH compared to N (Figure 4B). Moreover, expression of m⁶A circRNAs was significantly downregulated compared with non-m⁶A circRNAs in hypoxia, suggesting that m⁶A 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 m⁶A abundance, and 107 downregulated circRNAs with decreased m⁶A 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 m⁶A abundance influence the availability of miRNAs to target genes, we selected differentially-expressed circRNAs with increased or decreased m⁶A 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).

M⁶A circXpo6 and m⁶A circTmtc3 were downregulated in PASMCs and PAECs in hypoxia

M⁶A 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). M⁶A 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 m⁶A, we performed M⁶A 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). m⁶A circXpo6 and m⁶A 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)

In this study, we identified the transcriptome-wide map of m⁶A circRNAs in hypoxia mediated pulmonary hypertension. On the whole, we found that m⁶A level in circRNAs was reduced in lungs when exposed to hypoxia. M⁶A circRNAs were mainly derived from single exons of protein-coding genes in N and HPH. M⁶A abundance in circRNAs was downregulated in hypoxia in vitro. M⁶A influenced the circRNA–miRNA–mRNA co-expression network in hypoxia. Moreover, circXpo6 and circTmtc3 were the novel identified circRNAs modified by m⁶A in hypoxia mediated pulmonary hypertension.

M⁶A plays important roles in various biological processes. M⁶A is associated with cancer progression, promoting the proliferation of cancer cells and contributing to the cancer stem cell self-renewal [18, 21]. Lipid accumulation was reduced in hepatic cells when m⁶A abundance in peroxisome proliferator-activator (PPaR) was decreased [34]. Enhanced m⁶A level of mRNA contributed to compensated cardiac hypertrophy [35]. Also, m⁶A modification of large intergenic noncoding RNA 1281 was necessary for mouse embryonic stem cells differentiation [36].

Although it has been reported that m⁶A mRNAs were influenced by hypoxia, there is no report about m⁶A circRNAs in HPH yet. Up to now, no consistent conclusion was reached about the link between m⁶A and hypoxia. Previous reports found that the m⁶A abundance in mRNA was increased under hypoxia stress in HEK293T cells and cardiomyocytes [37, 38]. The increased m⁶A level stabilized the mRNAs of Glucose Transporter 1 (Glut1), Myc proto-oncogene bHLH transcription factor (Myc), Dual Specificity Protein Phosphatase 1 (Dusp1), Hairy and Enhancer of Split 1 (Hes1), and Jun Proto-Oncogene AP-1 Transcription Factor Subunit (Jun) without influencing their protein level [37]. In contrast, another reported that m⁶A level of total mRNA was decreased when human breast cancer cell lines were exposed to 1% O₂[26]. Hypoxia increased demethylation by stimulating hypoxia-inducible factor (HIF)-1α- and HIF-2α–dependent over-expression of ALKBH5 [26]. In addition, transcription factor EB activates the transcription of ALKBH5 and downregulates the stability of METTL3 mRNA in hypoxia/reoxygenation-induced autophagy in ischemic diseases [38]. Our study found that m⁶A abundance in total circRNAs was decreased in hypoxia exposure. Moreover, our study indicated that circXpo6 and circTmtc3 were the novel identified circRNAs modified by m⁶A in HPH. M⁶A abundance in circXpo6 and circTmtc3 was decreased in hypoxia. It is probably because of HIF-dependent and ALKBH5-mediated m⁶A demethylation [26].

Previous reports indicated that m⁶A methylation close to 3’UTR and stop codon of mRNA is inversely correlated with gene expression [14, 39]. Low m⁶A level is negatively associated with circRNAs expression, while high m⁶A level is not linked to circRNAs expression in human embryonic stem cells and HeLa cells [14]. Consistent with the previous reports [14, 39], our study found that m⁶A reduced the total circRNAs abundance in hypoxia. The association between m⁶A level and specific gene abundance is remained as an open question. Some previous reports indicated that m⁶A level was positively associated with long non-coding RNA (lncRNA) or mRNA expression [40, 41]. M⁶A was positively associated with RP11-138 J23.1 (RP11) expression when ALKBH5 was overexpressed in colorectal cancer [40]. mRNAs were downregulated after METTL14 deletion in b-cells [41]. On the contrary, another reports insisted that m⁶A level was negatively associated with RNA expression [42-44]. the mRNA lifetime of Family with Sequence Similarity 134, Member B (FAM134B) was prolonged when the m⁶A site was mutant [42]. The decreased m⁶A level resulted in the increased expression of N-methyl-D-aspartate receptor 1 (NMDAR1) in Parkinson's disease [43]. Forkhead Box protein M1 (FOXM1) abundance was increased when ALKBH5 was upregulated in glioblastoma [44]. Our study indicated that the expression of circXpo6 and circTmtc3 was decreased with the downregulated m⁶A level. The association between m⁶A level and circRNAs abundance was not determined yet. We suspected that m⁶A may influence the expression of circXpo6 and circTmtc3 through similar manners as before [40, 41]. But it needs further validation.

Competing endogenous RNA (CeRNA) mechanism was proposed that mRNAs, pseudogenes, lncRNAs and circRNAs interact with each other by competitive binding to miRNA response elements (MREs) [45, 46]. M⁶A acts as a post-transcript regulation of circRNAs and influences circRNAs expression, thus we suggested that m⁶A could also regulate the circRNA–miRNA–mRNA co-expression network. When the circRNAs were classified, we found that these downstream targets regulated by circRNA–miRNA of interest were mostly enriched in PH-associated Wnt and FoxO signaling pathways [30, 31]. The Wnt/β-catenin (bC) pathway and Wnt/ planar cell polarity (PCP) pathway are the two most critical Wnt signaling pathways in PH [30]. As known, the two important cells associated with HPH are PASMCs and PAECs [1, 3]. The growth of PASMCs was increased when Wnt/bC and Wnt/PCP pathways were activated by platelet derived growth factor beta polypeptide b (PDGF-BB) [30, 47]. In addition, the proliferation of PAECs was enhanced when Wnt/bC and Wnt/PCP pathways were activated by bone morphogenetic protein 2 (BMP2). Furthermore, the FoxO signaling pathway is associated with the apoptosis-resistant and hyper-proliferative phenotype of PASMCs [31]. Reactive oxygen species is increased by hypoxia and activates AMPK-dependent regulation of FoxO1 expression, resulting in increased expression of catalase in PASMCs [48]. Our study firstly uncovered that m⁶A influenced the stability of circRNAs, thus affecting the binding of circRNAs and miRNA, resulting in the activation of Wnt and FoxO signaling pathways.

In conclusion, our study firstly identified the transcriptome-wide map of m⁶A circRNAs in HPH. M⁶A level in circRNAs was decreased in lungs of HPH and in PASMCs and PAECs exposed to hypoxia. M⁶A level influenced circRNA–miRNA–mRNA co-expression network in HPH. Moreover, we firstly identified two downregulated m⁶A circRNAs in HPH: circXpo6 and circTmtc3. CircRNAs may be used as biomarkers because it is differentially enriched in specific cell types or tissues and not easily degraded [6]. Also, the aberrant m⁶A methylation may contribute to tumor formation and m⁶A RNAs may be a potential therapy target for tumor [17].

Limitations still exist in the study. First, we did not analyze the m⁶A level between circRNAs and the host genes. Second, the exact mechanism of hypoxia influences m⁶A was not demonstrated. Thirdly, the function of circXpo6 and circTmtc3 in HPH was not elaborated. Lastly, besides hypoxia mediated pulmonary hypertension, many other significant PH models should also be noted, such as monocrotaline mediated PH, monocrotaline + pneumonectomy mediated PH, and so on. It is insufficient that we explored the expression profiling of m⁶A circRNAs only in hypoxia mediated pulmonary hypertension. We plan to explore the expression profiling of m⁶A circRNAs in monocrotaline-induced PH and other PH models. Moreover, the clinical significance of m⁶A circRNAs for HPH should be further validated.

Hypoxia mediated PH rat model and measurement of RVSP and RVH

Sprague-Dawley rats (SPF, male, 180-200 g, 4 weeks) were obtained from the Animal Experimental Center of Zhejiang University, China. Rats were maintained in a normobaric normoxia (FiO₂ 21%, n = 6) or hypoxic chamber (FiO₂ 10%, n = 6) for 3 weeks [3, 49]. Rats were then anesthetized with pentobarbital sodium (130 mg/kg, ip) [50]. Then, rats were fixed in supine position on the board. All of the operations were performed after rats were anesthetized and became unconscious. RVSP was measured as below. Right ventricle catheterization was performed through the right jugular using a pressure-volume loop catheter (Millar) as the previous reports [49, 51]. After measurement of RVSP, rats were sacrificed and heart tissues were removed and segregated. The ratio of [RV/ (LV + S)] was used as an index of RVH. Lung were removed and immediately frozen at liquid nitrogen or fixed in 4% buffered paraformaldehyde solution. All experimental procedures were conducted in line with the principles approved by the Institutional Animal Care and Use Committee of Zhejiang University.

Histological analysis

Lung tissues were embedded in paraffin, sectioned at 4 μm and stained with hematoxylin and eosin (H&E) and a-smooth muscle actin (a-SMA, 1:100, ab124964, Abcam, USA). The ratio of pulmonary small artery wall thickness and muscularization were calculated [3].

Isolation and hypoxia-treatment of PASMCs and PAECs

PASMCs and PAECs were isolated using the methods according to previous reports [32, 50, 52]. PASMCs and PAECs were cultured in Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum (FBS) and 20% FBS for 48h, respectively [32, 53]. The cells were incubated in a 37°C, 21% O₂ or 1% O₂–5% CO₂ humidified incubator. PASMCs at 70–80% confluence in 4 to 7 passages were used in experiments. PAECs at 80–90% confluence in 4 to 5 passages were used in experiments [54].

RNA isolation and RNA-seq analysis of circRNAs

Total RNA (10 mg) was obtained using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) from lungs (1 g) of control and HPH rats. The extracted RNAs were purified with Rnase R (RNR07250, Epicentre) digestion to remove linear transcripts. Paired-end reads were harvested from Illumina Hiseq Sequence after quality filtering. The reads were aligned to the reference genome (UCSC RN5) with STAR software. CircRNAs were detected and annotated with CIRI software [55]. Raw junction reads were normalized to per million number of reads mapped to the genome with log2 scaled.

MeRIP and Library Preparation

Total RNA was extracted as the methods described above. Then, rRNA was depleted following DNase I treatment. RNase R treatment (5 units/mg) was performed in duplicate with 5 mg of rRNA-depleted RNA input. High-throughput m⁶A and circRNAs sequencing were performed by Cloudseq Biotech Inc. (Shanghai, China). Fragmented RNA was incubated with anti-m⁶A polyclonal antibody (Synaptic Systems, 202003) in IPP buffer for 2 hours at 4°C. The mixture was then incubated with protein A/G magnetic beads (88802, Thermo Fisher) at 4°C for an additional 2 hours. Then, bound RNA was eluted from the beads with N⁶-methyladenosine (PR3732, BERRY & ASSOCIATES) in IPP buffer and extracted with Trizol reagent (15596026, Thermo Fisher). NEBNext® Ultra™ RNA Library Prep Kit (E7530L, NEB) was used to construct RNA-seq library from immunoprecipitated RNA and input RNA. The m⁶A-IP and input samples were subjected to 150 bp paired-end sequencing on Illumina HiSeq sequencer. Methylated sites on circRNAs were identified by MetPeak software.

Construction of circRNA–miRNA–mRNA co-expression network

The circRNA–miRNA–mRNA co-expression network was based on the ceRNA theory that circRNA and mRNA shared the same MREs [45, 46]. Cytoscape was used to visualize the circRNA–miRNA–mRNA interactions based on the RNA-seq data. The circRNA-miRNA interaction and miRNA–mRNA interaction of interest were predicted by TargetScan and miRanda.

Measurement of Total m⁶A, MeRIP-RT-PCR and MeRIP-qRT-PCR

Total m⁶A content was measured in 200 ng aliquots of total RNA extracted from PASMCs and PAECs exposed to 21% O₂and 1% O₂for 48 h using an m⁶A RNA methylation quantification kit (P-9005, Epigentek). MeRIP (17-701, Millipore) was performed according to the manufacturer’s instruction. A 1.5 g aliquot of anti-m⁶A antibody (ABE572, Millipore) or anti-IgG (PP64B, Millipore) was conjugated to protein A/G magnetic beads overnight at 4^◦C. A 100 ng aliquot of total RNA was then incubated with the antibody in IP buffer supplemented with RNase inhibitor and protease inhibitor. The RNA complexes were isolated through phenol-chloroform extraction (P1025, Solarbio) and analyzed via RT-PCR or qRT-PCR assays. Primers sequences are listed in Table 1.

Data analysis

3’ adaptor-trimming and low quality reads were removed by cutadapt software (v1.9.3). Differentially methylated sites were identified by the R MeTDiff package. The read alignments on genome could be visualized using the tool IGV. Differentially expressed circRNAs were identified by Student’s t-test. GO and KEGG pathway enrichment analysis were performed for the corresponding parental mRNAs of the DE circRNAs. GO enrichment analysis was performed using the R topGO package. KEGG pathway enrichment analysis was performed according to a previous report [56]. GO analysis included BP analysis, CC analysis, and MF analysis. MicroRNAs sponged by the target genes were predicted by TargetScan and microRNA websites. P values are calculated by DAVID tool for GO and KEGG pathway analysis. The rest statistical analyses were performed with SPSS 19.0 (Chicago, IL, USA) and GraphPad Prism 5 software (La Jolla, CA). N refers to number of samples in figure legends. The statistical significance was determined by Student’s t-test (two-tailed) or two-sided Wilcoxon-Mann-Whiteney test. P < 0.05 was considered statistically significant. All experiments were independently repeated at least three times.

PH Pulmonary hypertension

HPH Hypoxia mediated pulmonary hypertension

COPD Chronic obstructive pulmonary disease

PAECs Pulmonary artery endothelial cells

PASMCs Pulmonary artery smooth muscle cells

RVH Right ventricular hypertrophy

circRNAs Circular RNAs

m⁶A N6-methyladenosine

ncRNAs Noncoding RNAs

METTL3, METTL14, METTL16 methyltransferase-like-3, -14 and -16

WTAP Wilms tumor 1-associated protein

RBM15 RNA binding motif protein 15

KIAA1429 Vir like m6A methyltransferase associated

ZC3H13 Zinc finger CCCH-type containing 13

FTO Fat mass and obesity related protein

ALKBH5 Alkylation repair homolog 5

YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2 YT521-B homology domain-containing proteins family

RVSP Right ventricular systolic pressure

GO Gene ontology analysis

KEGG Kyoto Encyclopedia of Genes and Genomes

BP Biological process analysis

CC Cellular component analysis

MF Molecular functions analysis

MeRIP M⁶A RNA Immunoprecipitation

PPaR Peroxisome proliferator-activator

Glut1 Glucose Transporter 1

Myc Myc proto-oncogene bHLH transcription factor

Dusp1 Dual Specificity Protein Phosphatase 1

Hes1 Hairy and Enhancer of Split 1

Jun Jun Proto-Oncogene AP-1 Transcription Factor Subunit

HIF-1α, HIF-2α Hypoxia-inducible factor-1α,-2α

LncRNA Long non-coding RNA

RP11 RP11-138 J23.1

FAM134B Family with Sequence Similarity 134, Member B

NMDAR1 N-methyl-D-aspartate receptor 1

FOXM1 Forkhead Box protein M1

CeRNA Competing endogenous RNA

MREs MiRNA response elements

PDGF-BB Platelet derived growth factor beta polypeptide b

BMP2 Bone morphogenetic protein 2

MeRIP-qRT-PCR MeRIP-quantitative RT-PCR

Ethics approval and consent to participate

The protocol of this research was submitted to and approved by the Institutional Animal Care and Use Committee of Zhejiang University, China (permit number: SYXK 2017-0006).

Consent to publish

Not Applicable.

Availability of data and materials

The datasets generated during the current study are available in the Zenodo repository, the persistent web links are https://doi.org/10.5281/zenodo.3596039, https://doi.org/10.5281/zenodo.3596178, https://doi.org/10.5281/zenodo.3596016, and https://doi.org/10.5281/zenodo.3596276.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was supported by the National Natural Science Foundation of China (81570043), the National Natural Science Foundation of China (81270107), the Zhejiang Provincial Natural Science Foundation of China (LY19H010003). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Authors’ Contributions

HS performed the experiments, analyzed the data, and made the figures; GWW, LFW and XQM analyzed the results; KJY and RFZ designed the study, analyzed the data, explained the findings and wrote the paper. All authors read and approved the final manuscript.

Acknowledgements

We thanked all subjects who participated in this study. We thank Cloud-Seq Biotech Ltd. Co. (Shanghai, China) for the MeRIP-Seq service and the subsequent bioinformatics analysis.

Pugliese SC, Poth JM, Fini MA, Olschewski A, El Kasmi KC,and Stenmark KR: The role of inflammation in hypoxic pulmonary hypertension: from cellular mechanisms to clinical phenotypes. Am J Physiol Lung Cell Mol Physiol 2015, 308:L229–L252.
Galiè N, Humbert M, Luc Vachiery J, Gibbs S, Lang I, Torbicki A, et al: 2015 ESC/ERS Guidelines for the diagnosis and treatment of pulmonary hypertension. The European respiratory journal 2015, 46:879–882.
Zhang RF, Shi LH, Zhou L, Zhang GS, Wu XH, Shao FC, et al: Transgelin as a therapeutic target to prevent hypoxic pulmonary hypertension. Am J Physiol Lung Cell Mol Physiol 2014, 306:L574–L583.
Ball MK, Waypa GB, Mungai PT, Nielsen JM, Czech L, Dudley VJ, et al: Regulation of hypoxia-induced pulmonary hypertension by vascular smooth muscle hypoxia-inducible factor-1alpha. American journal of respiratory and critical care medicine 2014, 189(3):314-324.
Meng X, Li X, Zhang P, Wang J, Zhou Y, Chen M: Circular RNA: an emerging key player in RNA world. Briefings in bioinformatics 2017, 18(4):547-557.
Salzman J: Circular RNA Expression: Its Potential Regulation and Function. Trends in genetics : TIG 2016, 32(5):309-316.
Chen LL, Yang L: Regulation of circRNA biogenesis. RNA biology 2015, 12(4):381-388.
Pamudurti NR, Bartok O, Jens M, Ashwal-Fluss R, Stottmeister C, Ruhe L, et al: Translation of CircRNAs. Molecular cell 2017, 66(1):1-13.
Yang Y, Fan X, Mao M, Song X, Wu P, Zhang Y, et al: Extensive translation of circular RNAs driven by N(6)-methyladenosine. Cell research 2017, 27(5):626-641.
Boeckel JN, Jae N, Heumuller AW, Chen W, Boon RA, Stellos K, et al: Identification and Characterization of Hypoxia-Regulated Endothelial Circular RNA. Circulation research 2015, 117(10):884-890.
Dang RY, Liu FL, Li Y: Circular RNA hsa_circ_0010729 regulates vascular endothelial cell proliferation and apoptosis by targeting the miR-186/HIF-1alpha axis. Biochemical and biophysical research communications 2017, 490(2):104-110.
Wang J, Zhu MC, Kalionis B, Wu JZ, Wang LL, Ge HY, et al: Characteristics of circular RNA expression in lung tissues from mice with hypoxiainduced pulmonary hypertension. International journal of molecular medicine 2018, 42(3):1353-1366.
Chen K, Lu Z, Wang X, Fu Y, Luo GZ, Liu N, et al: High-resolution N(6) -methyladenosine (m(6) A) map using photo-crosslinking-assisted m(6) A sequencing. Angewandte Chemie 2015, 54(5):1587-1590.
Zhou C, Molinie B, Daneshvar K, Pondick JV, Wang J, Van Wittenberghe N, et al: Genome-Wide Maps of m6A circRNAs Identify Widespread and Cell-Type-Specific Methylation Patterns that Are Distinct from mRNAs. Cell reports 2017, 20(9):2262-2276.
Zheng Y, Nie P, Peng D, He Z, Liu M, Xie Y, et al: m6AVar: a database of functional variants involved in m6A modification. Nucleic acids research 2018, 46(D1):D139-D145.
Xuan JJ, Sun WJ, Lin PH, Zhou KR, Liu S, Zheng LL, et al: RMBase v2.0: deciphering the map of RNA modifications from epitranscriptome sequencing data. Nucleic acids research 2018, 46(D1):D327-D334.
Wang S, Chai P, Jia R, Jia R: Novel insights on m(6)A RNA methylation in tumorigenesis: a double-edged sword. Molecular cancer 2018, 17(1):101.
Dai D, Wang H, Zhu L, Jin H, Wang X: N6-methyladenosine links RNA metabolism to cancer progression. Cell Death Dis 2018, 9(2):124.
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR: Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell 2012, 149(7):1635-1646.
Wang X, Huang J, Zou T, Yin P: Human m(6)A writers: Two subunits, 2 roles. RNA biology 2017, 14(3):300-304.
Deng X, Su R, Weng H, Huang H, Li Z, Chen J: RNA N(6)-methyladenosine modification in cancers: current status and perspectives. Cell research 2018, 28(5):507-517.
Wu B, Li L, Huang Y, Ma J, Min J: Readers, writers and erasers of N(6)-methylated adenosine modification. Current opinion in structural biology 2017, 47:67-76.
Liao S, Sun H, Xu C: YTH Domain: A Family of N(6)-methyladenosine (m(6)A) Readers. Genomics, proteomics & bioinformatics 2018, 16(2):99-107.
Patil DP, Pickering BF, Jaffrey SR: Reading m(6)A in the Transcriptome: m(6)A-Binding Proteins. Trends in cell biology 2018, 28(2):113-127.
Panneerdoss S, Eedunuri VK, Yadav P, Timilsina S, Rajamanickam S, Viswanadhapalli S, et al: Cross-talk among writers, readers, and erasers of m6A regulates cancer growth and progression. SCIENCE ADVANCES 2018, 4:eaar8263.
Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, et al: Hypoxia induces the breast cancer stem cell phenotype by HIF-dependent and ALKBH5-mediated m(6)A-demethylation of NANOG mRNA. Proc Natl Acad Sci 2016, 113(14):E2047-2056.
Rowan SC, Keane MP, Gaine S, McLoughlin P: Hypoxic pulmonary hypertension in chronic lung diseases: novel vasoconstrictor pathways. The Lancet Respiratory Medicine 2016, 4(3):225-236.
Chan JJ, Tay Y: Noncoding RNA:RNA Regulatory Networks in Cancer. International journal of molecular sciences 2018, 19(5):1310.
Baarsma HA, Konigshoff M: 'WNT-er is coming': WNT signalling in chronic lung diseases. Thorax 2017, 72(8):746-759.
de Jesus Perez V, Yuan K, Alastalo TP, Spiekerkoetter E, Rabinovitch M: Targeting the Wnt signaling pathways in pulmonary arterial hypertension. Drug discovery today 2014, 19(8):1270-1276.
Savai R, Al-Tamari HM, Sedding D, Kojonazarov B, Muecke C, Teske R, et al: Pro-proliferative and inflammatory signaling converge on FoxO1 transcription factor in pulmonary hypertension. Nature medicine 2014, 20(11):1289-1300.
Su H, Xu X, Yan C, Shi Y, Hu Y, Dong L, et al: LncRNA H19 promotes the proliferation of pulmonary artery smooth muscle cells through AT1R via sponging let-7b in monocrotaline-induced pulmonary arterial hypertension. Respir Res 2018, 19(1):254.
Caruso P, MacLean MR, Khanin R, McClure J, Soon E, Southgate M, et al: Dynamic changes in lung microRNA profiles during the development of pulmonary hypertension due to chronic hypoxia and monocrotaline. Arteriosclerosis, thrombosis, and vascular biology 2010, 30(4):716-723.
Zhong X, Yu J, Frazier K, Weng X, Li Y, Cham CM, et al: Circadian Clock Regulation of Hepatic Lipid Metabolism by Modulation of m(6)A mRNA Methylation. Cell reports 2018, 25(7):1816–1828.
Dorn LE, Lasman L, Chen J, Xu X, Hund TJ, Medvedovic M, et al: The N(6)-Methyladenosine mRNA Methylase METTL3 Controls Cardiac Homeostasis and Hypertrophy. Circulation 2019, 139(4):533-545.
Yang DD, Qiao J, Wang GY, Lan YY, Li GP, Guo XD, et al: N⁶-Methyladenosine modification of lincRNA 1281 is critically required for mESC differentiation potential, vol. 46; 2018.
Nate J. Fry BAL, Olga R. Ilkayeva, Christopher L. Holley, and Kyle D. Mansfield: N6-methyladenosine is required for the hypoxic stabilization of specific mRNAs. Rna 2017, 23(9):1444-1455.
Song H, Feng X, Zhang H, Luo Y, Huang J, Lin M, et al: METTL3 and ALKBH5 oppositely regulate m(6)A modification of TFEB mRNA, which dictates the fate of hypoxia/reoxygenation-treated cardiomyocytes. Autophagy 2019:1-19.
Luo GZ, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, et al: Unique features of the m6A methylome in Arabidopsis thaliana. Nature communications 2014, 5:5630.
Wu Y, Yang X, Chen Z, Tian L, Jiang G, Chen F, et al: m(6)A-induced lncRNA RP11 triggers the dissemination of colorectal cancer cells via upregulation of Zeb1. Molecular cancer 2019, 18(1):87.
Yao Y, Bi Z, Wu R, Zhao Y, Liu Y, Liu Q, et al: METTL3 inhibits BMSC adipogenic differentiation by targeting the JAK1/STAT5/C/EBPbeta pathway via an m(6)A-YTHDF2-dependent manner. FASEB journal : official publication of the Federation of American Societies for Experimental Biology 2019:fj201802644R.
Cai M, Liu Q, Jiang Q, Wu R, Wang X, Wang Y: Loss of m(6) A on FAM134B promotes adipogenesis in porcine adipocytes through m(6) A-YTHDF2-dependent way. IUBMB life 2019, 71(5):580-586.
Chen X, Yu C, Guo M, Zheng X, Ali S, Huang H, et al: Down-Regulation of m6A mRNA Methylation Is Involved in Dopaminergic Neuronal Death. ACS chemical neuroscience 2019.
Fitzsimmons CM, Batista PJ: It's complicated... m(6)A-dependent regulation of gene expression in cancer. Biochimica et biophysica acta Gene regulatory mechanisms 2019, 1862(3):382-393.
Chen L, Zhang S, Wu J, Cui J, Zhong L, Zeng L, Ge S: circRNA_100290 plays a role in oral cancer by functioning as a sponge of the miR-29 family. Oncogene 2017, 36(32):4551-4561.
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP: A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell 2011, 146(3):353-358.
J West JD, Austin ED, Gaskill C, Baskir R, Bilousova G, Jean JC, et al: Identification of a common Wnt-associated genetic signature across multiple cell types in pulmonary arterial hypertension. Am J Physiol Cell Physiol 2014, 307:C415–C430.
Awad H, Nolette N, Hinton M, Dakshinamurti S: AMPK and FoxO1 regulate catalase expression in hypoxic pulmonary arterial smooth muscle. Pediatric pulmonology 2014, 49(9):885-897.
Cowburn AS, Crosby A, Macias D, Branco C, Colaco RD, Southwood M, et al: HIF2alpha-arginase axis is essential for the development of pulmonary hypertension. Proceedings of the National Academy of Sciences of the United States of America 2016, 113(31):8801-8806.
Yun X, Jiang HY, Lai N, Wang J, and Shimoda LA: Aquaporin 1-mediated changes in pulmonary arterial smooth muscle cell migration and proliferation involve β-catenin. Am J Physiol Lung Cell Mol Physiol 2017, 313:L889–L898.
Yamazato Y, Ferreira AJ, Hong KH, Sriramula S, Francis J, Yamazato M, et al: Prevention of pulmonary hypertension by Angiotensin-converting enzyme 2 gene transfer. Hypertension 2009, 54(2):365-371.
Xiao R, Su Y, Feng T, Sun M, Liu B, Zhang J, et al: Monocrotaline Induces Endothelial Injury and Pulmonary Hypertension by Targeting the Extracellular Calcium-Sensing Receptor. Journal of the American Heart Association 2017, 6(4):e004865.
Zhang HY, Liu Y, Yan LX, Du W, Zhang XD, Zhang M, et al: Bone morphogenetic protein-7 inhibits endothelial-mesenchymal transition in pulmonary artery endothelial cell under hypoxia. Journal of Cellular Physiology 2018, 233(5):4077-4090.
Omura J, Satoh K, Kikuchi N, Satoh T, Kurosawa R, Nogi M, et al: Protective Roles of Endothelial AMP-Activated Protein Kinase Against Hypoxia-Induced Pulmonary Hypertension in Mice. Circulation research 2016, 119(2):197-209.
Thomas B. Hansen MTV, Christian K. Damgaard and Jørgen Kjems: Comparison of circular RNA prediction tools. Nucleic acids research 2015, 44(6):e58.
Tian L, Greenberg SA, Kong SW, Altschuler J, Kohane IS, Park PJ: Discovering statistically significant pathways in expression profiling studies. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(38):13544-13549.

Table 1. Primers for RT-PCR or qRT-PCR

Name	Sequence	Product size(bp)
linear Xpo6	Sense: 5’CTGTGTTTTGGGTCAGGAGC3’ Antisense: 5’ATCGAGTTCCTCTAGCCTGC3’	199
linear Tmtc3	Sense: 5’ACTCTGCTGTGATTGGACCA3’ Antisense:5’AGAAGAGGTTTGATGCGGGA3’	203
circXpo6	Sense: 5’TCTGGGAGACAAGGAAGCAG3’ Antisense: 5’CAGGATGGGGATGGGCTG3’	101
circTmtc3	Sense: 5’TACCCATGTTCAGCCAGGTT3’ Antisense:5’GAAGCCAAGCATTCACAGGA3’	100

Additional file 1:

Data S1. Differentially expressed m⁶A abundance in circRNAs.

Additional file 2:

Data S2. Differentially expressed m⁶A abundance linked with differentially expressed circRNAs abundance.

Additional file 3:

Differentially expressed m⁶A circRNAs in the lungs of HPH rat model.

Additional file 4:

Differentially expressed circRNAs in the lungs of HPH rat model.

PDF

Journal Publication

published 13 Jan, 2020

Read the published version in BMC Genomics →

Editor assigned by journal
03 Jan, 2020
Submission checks completed at journal
02 Jan, 2020
Editor invited by journal
02 Jan, 2020

You are reading this latest preprint version

Transcriptome-wide map of m6A circRNAs identified in a rat model of hypoxia mediated pulmonary hypertension

Status:

Journal Publication

Version 4

Abstract

Figures

Background

Results

Discussion

Conclusion

Methods

Abbreviations

Declarations

References

Table

Additional File Legends

Supplementary Files

Status:

Journal Publication

Version 4