Abstract

The incidence of hypoxic pulmonary hypertension (HPH) is increasing. Accumulating evidence suggests that long noncoding RNAs (lncRNAs) play an important role in HPH, but the functions and mechanism have yet to be fully elucidated. In the present study, we established a HPH rat model with 8 h of hypoxia exposure (10% O2) per day for 21 days. High-throughput sequencing identified 60 differentially expressed (DE) lncRNAs, 20 DE miRNAs and 695 DE mRNAs in rat lung tissue. qRT-PCR verified the accuracy of the results. The DE mRNAs were significantly enriched in immune response, inflammatory response, leukocyte migration, cell cycle, cellular response to interleukin-1, IL-17 signalling pathway, cytokine–cytokine receptor interaction and Toll-like receptor signalling pathway. According to the theory of competing endogenous RNA (ceRNA) networks, lncRNA–miRNA–mRNA network was constructed by Cytoscape software, 16 miRNAs and 144 mRNAs. The results suggested that seven DE lncRNAs (Ly6l, AABR07038849.2, AABR07069008.2, AABR07064873.1, AABR07001382.1, AABR07068161.1 and AABR07060341.2) may serve as molecular sponges of the corresponding miRNAs and play a major role in HPH.

Introduction

Hypoxic pulmonary hypertension (HPH) is generally characterized as atypical contraction and remodelling of the small pulmonary vessels caused by abnormal activation of the local pulmonary vasoconstriction system and proliferation of pulmonary smooth muscle cells [1]. HPH is common in almost all advanced stages of chronic lung disease, and the prevalence of HPH has increased in both adults and children in recent years [2]. Although many studies have investigated HPH, its pathogenesis has not been fully elucidated and must be explored from a new perspective.

Long noncoding RNAs (lncRNAs), which are more than 200 nucleotides in length, were once considered unrelated transcriptional waste with no biological function [3]; however, the competing endogenous RNA (ceRNA) hypothesis indicates that mRNA, transcriptional pseudogenes and lncRNAs interact with each other by competing for microRNA (miRNA) binding [4]. LncRNAs, including ceRNAs, have been found to be widely involved in the biological processes of many diseases by regulating the amount and function of miRNAs, thus affecting downstream signalling pathways [5,6]. In addition, many lncRNAs have been shown to be differentially expressed (DE) in HPH rats compared with healthy rats [7]. Subsequently, lncRNA-MEG3 and lncRNA-Tug1 have been shown to contribute to HPH as ceRNAs [8,9]. However, no studies have comprehensively analysed the correlation among lncRNAs, miRNAs and mRNAs in HPH. Thus, construction of lncRNA–miRNA–mRNA ceRNA networks is necessary. We screened for DE RNAs through high-throughput sequencing and constructed lncRNA–miRNA–mRNA ceRNA networks. Our study provides a novel theoretical basis for the pathogenesis of HPH from the perspective of ceRNA.

Materials and methods

Animal model

Twelve male specific pathogen-free Sprague-Dawley rats aged 6-8 weeks (150–200 g) were purchased from the Animal Centre of Kunming Medical University. The animal experiments were carried out in the Experimental Animal Center of Kunming Medical University and approved by the Ethics Committee of Kunming Medical University (KMMU2020207).

After a 3-day adaptation period, rats were randomly divided into the HPH group and the healthy group. According to previously reported methods [10], rats in the HPH group were exposed to normal pressure hypoxia (9.5–10.5% O2) 8 h a day for 21 days to construct the HPH rat model. Rats in the healthy group were housed under normoxic conditions.

RVSP, RVHI and WA (%) measurements

After 3 weeks in hypoxic or normoxic environments, rats were injected with sodium pentobarbital. When the anaesthetized rats had a decreased respiratory rate, complete muscle relaxation, disappearance of the eyelid reflex and a weak corneal reflex, the experiment was started. A polyethylene catheter was inserted from the right jugular vein to the right ventricle, then the right ventricle systolic pressure (RVSP) was recorded using the BL-420F biometric function test system. After the rats were killed by exsanguination under surgical levels of anaesthesia, lung and heart tissues were separated and collected. The right lung tissue was divided into two parts as follows: one part was used for sequencing, and one part was stored at −80°C. The left lung was saturated with 4% formaldehyde. Atrial tissue was removed from the collected heart, and the right ventricle, left ventricle and ventricular septum were separated. The index of right ventricular hypertrophy (RVHI) was measured using the following equation: RVHI = right ventricle / (left ventricle + ventricular septum). After fixation in 4% formaldehyde for 2 weeks, lung tissues were dehydrated with different concentrations xylene and then embedded with paraffin. Lung sections were stained with hematoxylin and eosin solutions, and then photographed using a light microscope. The ratio of wall area (WA%) of the pulmonary arterioles (approximate diameter of 100 µm and concomitant with terminal bronchiole) was calculated using the following formula to reflect the remodelling of the pulmonary microvasculature: WA% = (1 – lumen area / vessel total area).

Sequencing and identification of DE lncRNAs and DE mRNAs

Three right lung tissues from each of the two groups were randomly selected for sequencing. Total RNA was extracted from the lung tissues using TRIzol Reagent (Ambion Co., Carlsbad, U.S.A.), and the RNA concentration and purity were determined by oligosaccharide electrophoresis and an Agilent 2100 Bioanalyzer. The rRNA in the total RNA was removed using a ribosomal RNA removal kit, and 200–300 bp RNA fragments were generated by ion interruption. The first strand of cDNA was synthesized using 6-base random primers and reverse transcriptase, and the second strand of cDNA was synthesized using the first strand of cDNA as the template. After library construction, PCR amplification was used for library fragment enrichment, and the library selection was then performed according to the fragment size with a library size ranging from 300 to 400 bp. An Agilent 2100 Bioanalyser was used to conduct quality inspection of the library, and paired-end sequencing was then performed based on an Illumina sequencing platform. The Ensembl database (http://www.ensembl.org/) was selected as the reference genome.

The raw data were filtered, and the resulting clean data were aligned to the rat reference genome. We also used CPC, CNCI and Pfamscan software to determine whether these new transcripts that are formed by resplicing of exons based on sequencing principles and are different from known RNAs that have the protein-coding ability to obtain highly reliable lncRNAs. We considered that the new transcripts with no coding potential determined by all three software programs were highly reliable lncRNAs for subsequent analysis. HTSeq statistics were used to obtain the read count value for each gene, and fragments per kilobase of transcript per million mapped reads (FPKM) was used to standardize lncRNA and mRNA expression. EdgeR [11], a software package of R for the comprehensive analysis of count data, was used to identify differentially expressed (DE) lncRNAs and DE mRNAs by a threshold P-value <0.05 and log2 fold change (FC) >1 or <1.

Sequencing and identification of DE miRNA

Total RNA was extracted from the randomly selected lung tissues using TRIzol Reagent (Ambion Co., Carlsbad, U.S.A.). The TruSeq Small RNA Sample Prep Kit (Illumina) was used to construct the small RNA library after detection of total RNA concentration and purity. After PCR amplification and addition of the sequencing joint and index, gel electrophoresis was used to purify the library. An Agilent 2100 Bioanalyzer and Quant-iT PicoGreen dsDNA Assay Kit were used for quality detection and quantitative analysis of the library, respectively. Finally, the optimal sample size was selected and sequenced using an Illumina sequencing platform.

The raw data were processed to remove joints and quality filtering, and duplicated sequences were removed to obtain clean data. The clean data were mapped in the miRbase database to identify known miRNAs. According to the sequence number of miRNA maturation, the miRNA read count value was statistically measured. Counts per million (CPM) values were used to investigate the gene expression patterns of the samples as a whole. DE miRNAs were identified by using edgeR with a threshold P-value <0.05 and an absolute value of log2 (FC) <1.

Functional enrichment analysis

Metascape (http://metascape.org/) [12] was used to perform Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis for DE mRNAs. Biological processes (BP), cellular components (CC) and molecular functions (MF) in GO analysis were performed. The GO and KEGG analysis results were considered statistically significant with a threshold P-value <0.05. The results were saved as a text file, and ggplot2 packages in R were used for visual analysis.

Protein–protein interaction network

Protein–protein interaction (PPI) networks of DE mRNAs were constructed by STRING (https://string-db.org/). The minimum required interaction score was 0.40, and the disconnected nodes in the network were removed.

Construction of the ceRNA regulatory network

According to the regulation theory of competing endogenous RNA (ceRNA) networks, lncRNAs, including miRNA sponges, compete with mRNA through miRNA response elements (MREs) [4]. Therefore, in a ceRNA regulatory network, lncRNAs and mRNAs should show the same change trends, while miRNAs should show an opposing change. The MiRanda database was used to obtain the relationship between DE lncRNAs and DE miRNAs as well as between DE miRNAs and DE mRNAs. To analyse the interaction between lncRNAs and miRNAs, DE lncRNAs were annotated through the Ensemble database. Only known DE lncRNAs were included in the ceRNA network. The DE lncRNA–DE miRNA interactions and DE miRNA–DE mRNA interactions were then combined into the DE lncRNA–DE miRNA–DE mRNA network. Cytoscape software (version 3.7.0) was used to build this ceRNA network.

Quantitative polymerase chain reaction (qPCR) validation

The ceRNAs with large differential expression in the lncRNA-miRNA-mRNA network were selected for validation by qPCR. Primers were designed using the NCBI Primer-BLAST website (https://www.ncbi.nlm.nih.gov/tools/primer-blast/), and the sequences are shown in Table 1. U6 and GAPDH were used as internal references for noncoding RNA and mRNA, respectively. Total RNA from the lung tissue was extracted using TRIzol (Ambion Co., Carlsbad, U.S.A.). Complementary deoxyribose nucleic acid was synthetized using a reverse transcription kit (TIANGEN Biotech Co., Ltd., Beijing, China). Quantitative PCR was performed using an ABI Prism 7300HT Real Time PCR amplification apparatus. The results were analysed with the 2−ΔΔCt method.

Table 1
Primers for quantitative real-time PCR (qRT-PCR)
Gene namesCategoryPrimer sequence (5′-3′)
Ly61 lncRNA Forward-GGCCTCTTCTGCATCCTGTT 
  Reverse-CCTGGGTCCTAAGCAAAGGG 
AABR07068161.1 lncRNA Forward-ATGAGCAGATTACGGGCAGG 
  Reverse-CTGGATGAGATGGAGACCGC 
AABR07001382.1 lncRNA Forward-GAGCCATCCATGAACTCGGT 
  Reverse-TCATTGCCTTTCCTATGCTGC 
miR-21-5p miRNA Forward-TAGCTTATCAGACTGATGTTGA 
  Universal downstream primer of kit 
miR-485-5p miRNA Forward-GAGGCTGGCCGTGATGAATTC 
  Universal downstream primer of kit 
Cxcl11 mRNA Forward-CGCGGAGTGAAAGTGGTCAA 
  Reverse-AGCCTTCAGGGTAACAATCACT 
Hoxa1 mRNA Forward-CAGGAAGCAGACCCACCAAG 
  Reverse-GACCCACGTAGCCATACTCTC 
Trim29 mRNA Forward-CTGAAAGGCTATCCCTCCCTC 
  Reverse-TGGCCGGTAGTGAGACAGTA 
Myadml2 mRNA Forward-TGTTGCCGCGGAAGATAGAG 
  Reverse-CAGCGCCTAAGTGCAGGTAT 
Gene namesCategoryPrimer sequence (5′-3′)
Ly61 lncRNA Forward-GGCCTCTTCTGCATCCTGTT 
  Reverse-CCTGGGTCCTAAGCAAAGGG 
AABR07068161.1 lncRNA Forward-ATGAGCAGATTACGGGCAGG 
  Reverse-CTGGATGAGATGGAGACCGC 
AABR07001382.1 lncRNA Forward-GAGCCATCCATGAACTCGGT 
  Reverse-TCATTGCCTTTCCTATGCTGC 
miR-21-5p miRNA Forward-TAGCTTATCAGACTGATGTTGA 
  Universal downstream primer of kit 
miR-485-5p miRNA Forward-GAGGCTGGCCGTGATGAATTC 
  Universal downstream primer of kit 
Cxcl11 mRNA Forward-CGCGGAGTGAAAGTGGTCAA 
  Reverse-AGCCTTCAGGGTAACAATCACT 
Hoxa1 mRNA Forward-CAGGAAGCAGACCCACCAAG 
  Reverse-GACCCACGTAGCCATACTCTC 
Trim29 mRNA Forward-CTGAAAGGCTATCCCTCCCTC 
  Reverse-TGGCCGGTAGTGAGACAGTA 
Myadml2 mRNA Forward-TGTTGCCGCGGAAGATAGAG 
  Reverse-CAGCGCCTAAGTGCAGGTAT 

Statistical analysis

All data are presented as the mean ± standard deviation. An independent sample t-test was used to compare data between different groups with SPSS17.0 software. P<0.05 was considered statistically significant.

Results

Successful construction of animal models

In the six rats randomly selected for RNA sequencing, the rats in the HPH group showed significantly increased RVSP and RVHI compared with the healthy group (15.00 ± 2.65 mmHg versus 24.00 ± 2.65 mmHg and 0.12 ± 0.01 versus 0.18 ± 0.02, respectively; Figure 1A,B). At the same time, the wall of small pulmonary vessels in the HPH group was significantly thickened compared to that in the healthy group (Figure 1C). We next measured WA to quantitatively describe the changes in pulmonary small vessel remodelling. The results suggested that WA values were higher in the hypoxic group (Figure 1D). Together, these results demonstrated that the HPH rat model was successfully constructed.

Successful construction of the HPH rat model

Figure 1
Successful construction of the HPH rat model

(A) Right ventricular pressure (RVSP) in HPH and healthy rats. (B) Right ventricular hypertrophy (RVHI) was measured using the following formula: RVHI = right ventricle / (left ventricle + ventricular septum). (C) Haematoxylin and eosin (HE) staining of pulmonary arterial sections. (D) Wall area (WA) was calculated using the following equation: WA = 1 – lumen area / vessel total area; *P<0.05 and **P<0.01 versus healthy group.

Figure 1
Successful construction of the HPH rat model

(A) Right ventricular pressure (RVSP) in HPH and healthy rats. (B) Right ventricular hypertrophy (RVHI) was measured using the following formula: RVHI = right ventricle / (left ventricle + ventricular septum). (C) Haematoxylin and eosin (HE) staining of pulmonary arterial sections. (D) Wall area (WA) was calculated using the following equation: WA = 1 – lumen area / vessel total area; *P<0.05 and **P<0.01 versus healthy group.

Identification of DE mRNAs, DE miRNAs and DE lncRNAs

We obtained expression profile data of mRNAs, miRNAs and lncRNAs through high-throughput sequencing, and all data were uploaded to GSE159668. Through comparison to the healthy group, 60 DE lncRNAs (31 up-regulated and 29 down-regulated), 20 DE miRNAs (15 up-regulated and 5 down-regulated) and 695 DE mRNAs (534 up-regulated and 161 down-regulated) were screened with the threshold P-value <0.05 and the absolute value of log2 (FC) <1 (Supplementary Table S1). The heatmap and volcano map of noncoding RNAs (ncRNAs) and mRNAs are shown in Figure 2A,B.

Identification of the differentially expressed (DE) genes in the HPH rat model

Figure 2
Identification of the differentially expressed (DE) genes in the HPH rat model

(A) Heatmaps for the DE lncRNAs, DE miRNAs and DE mRNAs. Up-regulated genes are shown in red, and down-regulated genes are shown in green. (B) Volcano plots of the DE lncRNAs, DE miRNAs and DE mRNAs. The red, blue and black dots represent genes that are up-regulated, down-regulated and not significantly differentially expressed, respectively.

Figure 2
Identification of the differentially expressed (DE) genes in the HPH rat model

(A) Heatmaps for the DE lncRNAs, DE miRNAs and DE mRNAs. Up-regulated genes are shown in red, and down-regulated genes are shown in green. (B) Volcano plots of the DE lncRNAs, DE miRNAs and DE mRNAs. The red, blue and black dots represent genes that are up-regulated, down-regulated and not significantly differentially expressed, respectively.

Functional analysis of DE lncRNAs, DE miRNAs and DE mRNAs

Of the 695 DE mRNAs, many genes have been previously shown to play a role in different types of pulmonary hypertension. For example, MMP7 may participate in vascular remodelling in pulmonary hypertension [13]. Although there are weak associations between serum CXCL13 and disease severity, CXCL13 is overexpressed in both idiopathic pulmonary artery hypertension and chronic thromboembolic pulmonary hypertension [14]. Among the 20 DE miRNAs, miR-210-3p [15] and miR-21-5p [16] were associated with pulmonary hypertension. Of the 60 DE lncRNAs, only 16 DE lncRNAs were annotated, but their functions remain unclear. Metascape was used to perform functional enrichment of the DE mRNAs (Supplementary Table S2 and Figure 3A). GO analysis showed enrichment in chromosome segregation, humoural immune response, inflammatory response, leukocyte migration, cell cycle and cellular response to interleukin-1 (Figure 3B). In addition, KEGG pathway analysis revealed enrichment in the cell cycle, IL-17 signalling pathway, cytokine–cytokine receptor interaction and Toll-like receptor signalling pathway (Figure 3C).

Functional enrichment of DE mRNAs

Figure 3
Functional enrichment of DE mRNAs

(A) Top 20 clusters and their representative enriched terms (coloured by cluster ID). (B) Top 20 enriched terms on GO analysis. (C) Top 20 enriched terms on KEGG pathway analysis; DE, differentially expressed; GO, Gene Ontology; KEGG, Kyoto Encyclopaedia of Genes and Genomes.

Figure 3
Functional enrichment of DE mRNAs

(A) Top 20 clusters and their representative enriched terms (coloured by cluster ID). (B) Top 20 enriched terms on GO analysis. (C) Top 20 enriched terms on KEGG pathway analysis; DE, differentially expressed; GO, Gene Ontology; KEGG, Kyoto Encyclopaedia of Genes and Genomes.

Construction of the lncRNA–miRNA–mRNA ceRNA network

The MiRanda dataset (http://www.microrna.org/microrna/) was used to predict the targets of DE ncRNAs among the DE RNAs. The results indicated that two up-regulated lncRNAs were predicted to have three down-regulated DE miRNAs, and the three down-regulated miRNAs were predicted to target 97 up-regulated DE mRNAs. Combining the lncRNA–miRNA pairs and miRNA–mRNA pairs, the ceRNAs of up-regulated lncRNAs–down-regulated miRNAs–up-regulated mRNAs were constructed in Cytoscape software (version 3.7.0) (Figure 4A). Using the same method, another ceRNA of down-regulated lncRNAs–up-regulated miRNAs–down-regulated mRNAs (including 5 lncRNAs, 13 miRNAs and 47 mRNAs, Figure 4B) was constructed.

Construction of the lncRNA–miRNA–mRNA ceRNA regulatory network

Figure 4
Construction of the lncRNA–miRNA–mRNA ceRNA regulatory network

(A) The up-regulated lncRNA-down-regulated miRNA-up-regulated mRNA ceRNA regulatory network included 2 lncRNAs, 3 miRNAs and 97 mRNAs. (B) The down-regulated lncRNA-up-regulated miRNA-down-regulated mRNA ceRNA regulatory network included 5 lncRNAs, 13 miRNAs and 47 mRNAs. Red nodes indicate up-regulated RNAs, and green nodes indicate down-regulated RNAs. Triangles, ovals and rectangles indicate DE lncRNAs, DE miRNAs and DE mRNAs, respectively; ceRNA, competing endogenous RNA; DE, differentially expressed; lncRNAs, long noncoding RNAs; miRNAs, microRNAs.

Figure 4
Construction of the lncRNA–miRNA–mRNA ceRNA regulatory network

(A) The up-regulated lncRNA-down-regulated miRNA-up-regulated mRNA ceRNA regulatory network included 2 lncRNAs, 3 miRNAs and 97 mRNAs. (B) The down-regulated lncRNA-up-regulated miRNA-down-regulated mRNA ceRNA regulatory network included 5 lncRNAs, 13 miRNAs and 47 mRNAs. Red nodes indicate up-regulated RNAs, and green nodes indicate down-regulated RNAs. Triangles, ovals and rectangles indicate DE lncRNAs, DE miRNAs and DE mRNAs, respectively; ceRNA, competing endogenous RNA; DE, differentially expressed; lncRNAs, long noncoding RNAs; miRNAs, microRNAs.

Functional analysis and PPI network of DE mRNAs in the lncRNA–miRNA–mRNA ceRNA network

The 144 DE mRNAs in the ceRNA network were associated with the cell cycle, inorganic molecular entity transmembrane transporter activity, small molecule binding, cytokine–cytokine receptor interaction and p53 signalling pathway (Figure 5A). In addition, 19 hub genes were identified in the PPI network with MCODE score >15 (Figure 5B) and included Nusap1, Cenpf, Pbk, Ndc80, Kif2c, Ccnb1, Rrm2, Birc5, Cdca8, Ccna2, Cdc6, Racgap1, Kif18b, Exo1, Ska1, Kifc1, Hjurp, Incenp and Asf1b.

Functional enrichment and PPI networks of DE mRNAs in the ceRNA networks

Figure 5
Functional enrichment and PPI networks of DE mRNAs in the ceRNA networks

(A) Top 5 enriched terms in BP, CC, MF and KEGG analyses. (B) The PPI network of DE mRNAs in the ceRNA networks; BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 5
Functional enrichment and PPI networks of DE mRNAs in the ceRNA networks

(A) Top 5 enriched terms in BP, CC, MF and KEGG analyses. (B) The PPI network of DE mRNAs in the ceRNA networks; BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes.

qPCR validation of the ceRNA network genes

Three lncRNAs (Ly61, AABR07001382.1 and AABR07068161.1), two miRNAs (miR-21-5p and miR-485-5p) and four mRNAs (Cxcl11, Hoxa1, Trim29 and Myadml2) were selected for validation by qPCR to confirm the sequencing results. The expression of the nine selected RNAs was consistent with the sequencing results (Figure 6). Ly61, AABR07001382.1, miR-485-5p and Myadml2 showed down-regulated expression in the sequencing results, which was confirmed by qPCR. The expression of AABR07068161.1, miR-21-5p, Cxcl11, Hoxa1 and Trim29 in HPH rats was higher in both the high-throughput sequencing results and qPCR validation results.

RNA expression in the HPH group compared with the healthy group

Figure 6
RNA expression in the HPH group compared with the healthy group

*P<0.05, **P<0.01 and ***P<0.001 versus healthy group.

Figure 6
RNA expression in the HPH group compared with the healthy group

*P<0.05, **P<0.01 and ***P<0.001 versus healthy group.

Discussion

HPH occurs frequently in patients with chronic obstructive pulmonary disease or idiopathic pulmonary fibrosis and in people living at high altitudes [17], and its incidence is increasing. Some early changes in HPH are difficult to detect with traditional research tools and expertise. Much effort has been afforded to finding new therapeutic targets and biomarkers for HPH [18–21]. Most studies, however, have focused only on mRNAs [18,20,21]. In recent years, lncRNAs have attracted increasing attention due to their important roles in many diseases, but the functions and mechanisms of many lncRNAs in HPH are not yet clear. According to the interactions of lncRNAs, miRNAs and mRNAs through miRNA recognition elements [4], a lncRNA–miRNA–mRNA ceRNA network was therefore constructed in the present study to explore the mechanism and reveal novel biomarkers of HPH.

In the present study, we established a HPH rat model via 8 h of hypoxia exposure (10% O2) per day for 21 days. Successful construction of a rat HPH model was confirmed by detection of RVSP, RVHI and WA. High-throughput sequencing was performed on the lung tissue to investigate mRNA, lncRNA and miRNA expression level changes in HPH rats versus normal rats. The results revealed 31 up-regulated lncRNAs, 29 down-regulated lncRNAs, 15 up-regulated miRNAs, 5 down-regulated miRNAs, 534 up-regulated mRNAs and 161 down-regulated mRNAs. The expression levels of these RNAs were significantly different between the HPH group and the healthy group, indicating the potential use of these ncRNAs in distinguishing HPH from healthy controls.

To further explore the molecular functions and related pathways of the DE RNAs, GO and KEGG pathways were analysed. The immune response, inflammatory response, cellular response to interleukins, cell cycle, IL-17 signalling pathway, Toll-like receptor signalling pathway and leukocyte migration were highly enriched. These results suggested that many immune and inflammatory response pathways participate in HPH [22,23].

Based on the ceRNA regulatory network theory, we predicted the targets of the DE ncRNAs in the miRbase database to identify possible interactions among the screened DE RNAs. A lncRNA–miRNA–mRNA network was then constructed, and it included 7 lncRNAs, 16 miRNAs and 144 mRNAs. Almost none of the DE lncRNAs we screened have been studied, but many other lncRNAs have been confirmed to be important regulators of PH [24]. LncRNAs TUG1 and NONRATT015587.2, as ceRNAs, have been reported to have an impact on vascular remodelling in PH [25,26]. LncRNA CASC2 and Hoxaas3 regulate PASMC proliferation and migration in HPH by influencing the expression of miR-222 and homeobox a3, respectively [27,28]. LncRNA PAXIP1-AS1 plays an important role in PH by affecting PH-specific transcriptional programmes [29]. These previous results suggest that multiple lncRNAs are involved in HPH through different pathways. Thus, we propose that the seven DE lncRNAs (Ly6l, AABR07038849.2, AABR07069008.2, AABR07064873.1, AABR07001382.1, AABR07068161.1 and AABR07060341.2) in the ceRNA network may serve as molecular sponges of the corresponding miRNAs to regulate mRNA expression and participate in the occurrence and development of HPH by regulating the cell cycle, inorganic molecular entity transmembrane transporter activity, small molecule binding, cytokine–cytokine receptor interaction, p53 signalling pathway and other processes. At present, the role of the cell cycle in the pathogenesis of HPH has been confirmed by many studies [8,30], while the role of other biological processes in HPH has not been studied. Our future research will focus on exploring the roles of DE RNAs in the ceRNA network and the enriched biological processes in HPH pathogenesis.

The functions of some DE miRNAs have been explored. MiR-136-5p has been confirmed to regulate the inflammatory response and is up-regulated in many inflammation-related diseases [31,32]. Joshi et al. found that miR-21-5p is up-regulated in a pulmonary arterial hypertension rat model using an unbiased quantitative miRNA microarray analysis [16], which was consistent with our results. MiR-210-5p promotes epithelial–mesenchymal transformation, an important contributor to the development of HPH, by inhibiting the expression of PIK3R5 [33]. Yan et al. reported that miR-23a is highly expressed in hypoxia-treated pulmonary artery smooth muscle cells (PASMCs) and plays an important role in hypoxia-induced phenotypic transformation of PASMCs [34]. In the past 5 years, miR-23a has been found to be a biomarker for idiopathic PH [35] and a possible therapeutic target for pulmonary arterial hypertension because it promotes cell proliferation and migration by targeting the BMPR2/Smad1 signal in hypoxia-induced PASMCs [36].

In addition to these DE miRNAs, many others have been shown to be involved in the formation of HPH, especially hub genes. Nusap1 regulates the Wnt/β-catenin and TGF-β signalling pathways [37,38], which have been repeatedly reported to be related to HPH pathogenesis [39,40]. Ccnb1 has been confirmed to affect the cell cycle [41], which was indicated in our functional enrichment analysis. Cxcl11 is associated with various inflammatory responses [42,43] and is a critical contributor to HPH.

In general, our study showed differences in the expression levels of lncRNAs, miRNAs and mRNAs between HPH and healthy rats. By constructing a lncRNA–miRNA–mRNA ceRNA network, we identified a regulatory relationship between ncRNAs and mRNAs in addition to their roles in HPH development. Our study contributes to a better understanding of the pathogenesis of HPH and the role of ncRNAs in disease progression.

Although the role of lncRNAs and miRNAs in HPH was successfully analysed and discussed, the present study also had several limitations. First, the sample size was too small (3 HPH rats and 3 healthy rats) to draw conclusions. Second, we only performed qPCR verification, and the results of functional analysis were not experimentally verified. Further research is needed to determine the molecular mechanisms by which ncRNAs play a role in HPH generation and development.

Data Availability

All data has been uploaded to GSE159668. The data sets used and/or analyzed in this study are available from the corresponding author on reasonable request.

Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

Funding

This research was supported by a grant from the National Natural Science Foundation of China [grant numbers 81760015 and 82160016]; the “Special and Joint Program” of Yunnan Provincial Science and Technology Department and Kunming Medical University [grant number 2018FE001 (-206)]; the Yunnan Health Training Project of High-level Talents [grant number H-2018095]; the Young Academic and Technical Leaders of Yunnan Province [grant number 2017HB053]; and Famous Doctors of High-level Talent Training Support Program of Yunnan Province.

CRediT Author Contribution

Jie Liu: Data curation, Software, Formal analysis, Supervision, Writing—original draft, Writing—review & editing. Yishu Deng: Writing—original draft. Zeqin Fan: Software, Formal analysis. Shuanglan Xu: Data curation, Software, Formal analysis. Li Wei: Data curation, Formal analysis, Methodology. Xiaoxian Huang: Software, Formal analysis. Xiqian Xing: Software, Formal analysis, Funding acquisition, Writing—original draft, Project administration, Writing—review & editing. Jiao Yang: Funding acquisition, Project administration, Writing—review & editing.

Abbreviations

     
  • DE

    differentially expressed

  •  
  • GO

    Gene Ontology

  •  
  • HPH

    hypoxic pulmonary hypertension

  •  
  • KEGG

    Kyoto Encyclopaedia of Genes and Genomes

  •  
  • lncRNA

    long noncoding RNA

  •  
  • MRE

    miRNA response element

  •  
  • PPI

    protein–protein interaction

References

1.
Liu
M.-L.
,
Xing
S.-J.
,
Liang
X.-Q.
,
Luo
Y.
,
Zhang
B.
,
Li
Z.-C.
et al.
(
2020
)
Reversal of hypoxic pulmonary hypertension by hypoxia-inducible overexpression of angiotensin- (1-7) in pulmonary endothelial cells
.
Mol. Ther. Methods Clin. Dev.
17
,
975
985
[PubMed]
2.
Wijeratne
T.D.
,
Lajkosz
K.
,
Brogly
S.B.
,
Lougheed
M.D.
,
Jiang
L.
,
Housin
A.
et al.
(
2018
)
Increasing incidence and prevalence of WHO groups 1–4 pulmonary hypertension: a Population-Based Cohort Study in Ontario, Canada
.
Circulation Cardiovasc. Quality Outcomes
11
,
e003973
3.
Wang
S.
,
Zhang
C.
and
Zhang
X.
(
2020
)
Downregulation of long non-coding RNA ANRIL promotes proliferation and migration in hypoxic human pulmonary artery smooth muscle cells
.
Mol. Med. Rep.
21
,
589
596
[PubMed]
4.
Salmena
L.
,
Poliseno
L.
,
Tay
Y.
,
Kats
L.
and
Pandolfi
P.P.
(
2011
)
A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?
Cell
146
,
353
358
[PubMed]
5.
Wang
C.
,
Zi
H.
,
Wang
Y.
,
Li
B.
,
Ge
Z.
and
Ren
X.
(
2020
)
LncRNA CASC15 promotes tumour progression through SOX4/Wnt/β-catenin signalling pathway in hepatocellular carcinoma
.
Artific. Cell. Nanomed. Biotechnol.
48
,
763
769
6.
Yuan
Z.
and
Wang
W.
(
2020
)
LncRNA SNHG4 regulates miR-10a/PTEN to inhibit the proliferation of acute myeloid leukemia cells
.
Hematology
25
,
160
164
[PubMed]
7.
Wang
X.
,
Yan
C.
,
Xu
X.
,
Dong
L.
,
Su
H.
,
Hu
Y.
et al.
(
2016
)
Long noncoding RNA expression profiles of hypoxic pulmonary hypertension rat model
.
Gene
579
,
23
28
[PubMed]
8.
Xing
Y.
,
Zheng
X.
,
Fu
Y.
,
Qi
J.
,
Li
M.
,
Ma
M.
et al.
(
2019
)
Long noncoding RNA-maternally expressed gene 3 contributes to hypoxic pulmonary hypertension
.
Mol. Therapy: J. Am. Soc. Gene Ther.
27
,
2166
2181
9.
Yang
L.
,
Liang
H.
,
Shen
L.
,
Guan
Z.
and
Meng
X.
(
2019
)
LncRNA Tug1 involves in the pulmonary vascular remodeling in mice with hypoxic pulmonary hypertension via the microRNA-374c-mediated Foxc1
.
Life Sci.
237
,
116769
[PubMed]
10.
Zhang
R.
,
Shi
L.
,
Zhou
L.
,
Zhang
G.
,
Wu
X.
,
Shao
F.
et al.
(
2014
)
Transgelin as a therapeutic target to prevent hypoxic pulmonary hypertension
.
Am. J. Physiol. Lung Cell. Mol. Physiol.
306
,
L574
L583
[PubMed]
11.
Robinson
M.D.
,
McCarthy
D.J.
and
Smyth
G.K.
(
2010
)
edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
26
,
139
140
[PubMed]
12.
Zhou
Y.
,
Zhou
B.
,
Pache
L.
,
Chang
M.
,
Khodabakhshi
A.H.
,
Tanaseichuk
O.
et al.
(
2019
)
Metascape provides a biologist-oriented resource for the analysis of systems-level datasets
.
Nat. Commun.
10
,
1523
[PubMed]
13.
Arvidsson
M.
,
Ahmed
A.
,
Bouzina
H.
and
Rådegran
G.
(
2019
)
Matrix metalloproteinase 7 in diagnosis and differentiation of pulmonary arterial hypertension
.
Pulm Circ.
9
,
2045894019895414
[PubMed]
14.
Olsson
K.M.
,
Olle
S.
,
Fuge
J.
,
Welte
T.
,
Hoeper
M.M.
,
Lerch
C.
et al.
(
2016
)
CXCL13 in idiopathic pulmonary arterial hypertension and chronic thromboembolic pulmonary hypertension
.
Respir. Res.
17
,
21
[PubMed]
15.
Herrera-Rivero
M.
,
Zhang
R.
,
Heilmann-Heimbach
S.
,
Mueller
A.
,
Bagci
S.
,
Dresbach
T.
et al.
(
2018
)
Circulating microRNAs are associated with Pulmonary Hypertension and Development of Chronic Lung Disease in Congenital Diaphragmatic Hernia
.
Sci. Rep.
8
,
10735
[PubMed]
16.
Joshi
S.R.
,
Dhagia
V.
,
Gairhe
S.
,
Edwards
J.G.
,
McMurtry
I.F.
and
Gupte
S.A.
(
2016
)
MicroRNA-140 is elevated and mitofusin-1 is downregulated in the right ventricle of the Sugen5416/hypoxia/normoxia model of pulmonary arterial hypertension
.
Am. J. Physiol. Heart Circ. Physiol.
311
,
H689
H698
[PubMed]
17.
Rowan
S.C.
,
Keane
M.P.
,
Gaine
S.
and
McLoughlin
P.
(
2016
)
Hypoxic pulmonary hypertension in chronic lung diseases: novel vasoconstrictor pathways
.
Lancet Respir. Med.
4
,
225
236
[PubMed]
18.
Bouvard
C.
,
Genet
N.
,
Phan
C.
,
Rode
B.
,
Thuillet
R.
,
Tu
L.
et al.
(
2020
)
Connexin-43 is a promising target for pulmonary hypertension due to hypoxaemic lung disease
.
Eur. Respir. J.
55
,
[PubMed]
19.
Fan
J.
,
Fan
X.
,
Guang
H.
,
Shan
X.
,
Tian
Q.
,
Zhang
F.
et al.
(
2020
)
Upregulation of miR-335-3p by NF-κB transcriptional regulation contributes to the induction of pulmonary arterial hypertension via APJ during hypoxia
.
Int. J. Biol. Sci.
16
,
515
528
[PubMed]
20.
Liu
M.
,
Galli
G.
,
Wang
Y.
,
Fan
Q.
,
Wang
Z.
,
Wang
X.
et al.
(
2020
)
Novel therapeutic targets for hypoxia-related cardiovascular diseases: the role of HIF-1
.
Front. Physiol.
11
,
774
[PubMed]
21.
Wang
J.
,
Wu
J.
,
Zhu
X.
,
Chen
J.
,
Zhao
J.
,
Xu
Y.
et al.
(
2020
)
Absence of the MFG-E8 gene prevents hypoxia-induced pulmonary hypertension in mice
.
J. Cell. Physiol.
22.
Hudalla
H.
,
Michael
Z.
,
Christodoulou
N.
,
Willis
G.R.
,
Fernandez-Gonzalez
A.
,
Filatava
E.J.
et al.
(
2019
)
Carbonic anhydrase inhibition ameliorates inflammation and experimental pulmonary hypertension
.
Am. J. Respir. Cell Mol. Biol.
61
,
512
524
[PubMed]
23.
Ogoshi
T.
,
Tsutsui
M.
,
Kido
T.
,
Sakanashi
M.
,
Naito
K.
,
Oda
K.
et al.
(
2018
)
Protective role of myelocytic nitric oxide synthases against hypoxic pulmonary hypertension in mice
.
Am. J. Respir. Crit. Care Med.
198
,
232
244
[PubMed]
24.
Jin
Q.
,
Zhao
Z.
,
Zhao
Q.
,
Yu
X.
,
Yan
L.
,
Zhang
Y.
et al.
(
2020
)
Long noncoding RNAs: emerging roles in pulmonary hypertension
.
Heart Fail. Rev.
25
,
795
815
[PubMed]
25.
Bonnet
S.
,
Boucherat
O.
,
Provencher
S.
and
Paulin
R.
(
2019
)
Early evidence for the role of lncRNA TUG1 in vascular remodelling in pulmonary hypertension
.
Can. J. Cardiol.
35
,
1433
1434
[PubMed]
26.
Sun
Z.
,
Liu
Y.
,
Yu
F.
,
Xu
Y.
,
Yanli
L.
and
Liu
N.
(
2019
)
Long non-coding RNA and mRNA profile analysis of metformin to reverse the pulmonary hypertension vascular remodeling induced by monocrotaline
.
Biomed. Pharmacother.
115
,
108933
27.
Han
Y.
,
Liu
Y.
,
Yang
C.
,
Gao
C.
,
Guo
X.
and
Cheng
J.
(
2020
)
LncRNA CASC2 inhibits hypoxia-induced pulmonary artery smooth muscle cell proliferation and migration by regulating the miR-222/ING5 axis
.
Cell. Mol. Biol. Lett.
25
,
21
[PubMed]
28.
Zhang
H.
,
Liu
Y.
,
Yan
L.
,
Wang
S.
,
Zhang
M.
,
Ma
C.
et al.
(
2019
)
Long noncoding RNA Hoxaas3 contributes to hypoxia-induced pulmonary artery smooth muscle cell proliferation
.
Cardiovasc. Res.
115
,
647
657
[PubMed]
29.
Jandl
K.
,
Thekkekara Puthenparampil
H.
,
Marsh
L.M.
,
Hoffmann
J.
,
Wilhelm
J.
,
Veith
C.
et al.
(
2019
)
Long non-coding RNAs influence the transcriptome in pulmonary arterial hypertension: the role of PAXIP1-AS1
.
J. Pathol.
247
,
357
370
[PubMed]
30.
Zhang
J.
,
Li
Y.
,
Qi
J.
,
Yu
X.
,
Ren
H.
,
Zhao
X.
et al.
(
2020
)
Circ-calm4 serves as an miR-337-3p sponge to regulate Myo10 (Myosin 10) and promote pulmonary artery smooth muscle proliferation
.
Hypertension (Dallas, Tex: 1979)
75
,
668
679
[PubMed]
31.
Deng
G.
,
Gao
Y.
,
Cen
Z.
,
He
J.
,
Cao
B.
,
Zeng
G.
et al.
(
2018
)
miR-136-5p regulates the inflammatory response by targeting the IKKβ/NF-κB/A20 pathway after spinal cord injury
.
Cell. Physiol. Biochem.: Int. J. Exp. Cell. Physiol. Biochem. Pharmacol.
50
,
512
524
32.
Zhong
Y.
,
Yu
C.
and
Qin
W.
(
2019
)
LncRNA SNHG14 promotes inflammatory response induced by cerebral ischemia/reperfusion injury through regulating miR-136-5p /ROCK1
.
Cancer Gene Ther.
26
,
234
247
[PubMed]
33.
Liu
W.
,
Jiang
D.
,
Gong
F.
,
Huang
Y.
,
Luo
Y.
,
Rong
Y.
et al.
(
2020
)
miR-210-5p promotes epithelial-mesenchymal transition by inhibiting PIK3R5 thereby activating oncogenic autophagy in osteosarcoma cells
.
Cell Death Dis.
11
,
93
[PubMed]
34.
Yan
L.
,
Gao
H.
,
Li
C.
,
Han
X.
and
Qi
X.
(
2017
)
Effect of miR-23a on anoxia-induced phenotypic transformation of smooth muscle cells of rat pulmonary arteries and regulatory mechanism
.
Oncol. Lett.
13
,
89
98
[PubMed]
35.
Sarrion
I.
,
Milian
L.
,
Juan
G.
,
Ramon
M.
,
Furest
I.
,
Carda
C.
et al.
(
2015
)
Role of circulating miRNAs as biomarkers in idiopathic pulmonary arterial hypertension: possible relevance of miR-23a
.
Oxid. Med. Cell. Longev.
2015
,
792846
36.
Zhang
Y.
,
Peng
B.
and
Han
Y.
(
2018
)
MiR-23a regulates the proliferation and migration of human pulmonary artery smooth muscle cells (HPASMCs) through targeting BMPR2/Smad1 signaling
.
Biomed. Pharmacother.
103
,
1279
1286
37.
Gao
S.
,
Yin
H.
,
Tong
H.
,
Zhan
K.
,
Yang
G.
,
Hossain
M.A.
et al.
(
2020
)
Nucleolar and spindle associated protein 1 (NUSAP1) promotes bladder cancer progression through the TGF-β signaling pathway
.
OncoTargets Ther.
13
,
813
825
[PubMed]
38.
Sun
L.
,
Shi
C.
,
Liu
S.
,
Zhang
E.
,
Yan
L.
,
Ji
C.
et al.
(
2020
)
Overexpression of NuSAP1 is predictive of an unfavourable prognosis and promotes proliferation and invasion of triple-negative breast cancer cells via the Wnt/β-catenin/EMT signalling axis
.
Gene
747
,
144657
[PubMed]
39.
Alapati
D.
,
Rong
M.
,
Chen
S.
,
Lin
C.
,
Li
Y.
and
Wu
S.
(
2013
)
Inhibition of LRP5/6-mediated Wnt/β-catenin signaling by Mesd attenuates hyperoxia-induced pulmonary hypertension in neonatal rats
.
Pediatr. Res.
73
,
719
725
[PubMed]
40.
Zhang
N.
,
Dong
M.
,
Luo
Y.
,
Zhao
F.
and
Li
Y.
(
2018
)
Danshensu prevents hypoxic pulmonary hypertension in rats by inhibiting the proliferation of pulmonary artery smooth muscle cells via TGF-β-smad3-associated pathway
.
Eur. J. Pharmacol.
820
,
1
7
[PubMed]
41.
Wang
F.
,
Chen
X.
,
Yu
X.
and
Lin
Q.
(
2019
)
Degradation of CCNB1 mediated by APC11 through UBA52 ubiquitination promotes cell cycle progression and proliferation of non-small cell lung cancer cells
.
Am. J. Transl. Res.
11
,
7166
7185
[PubMed]
42.
Kochumon
S.
,
Madhoun
A.A.
,
Al-Rashed
F.
,
Azim
R.
,
Al-Ozairi
E.
,
Al-Mulla
F.
et al.
(
2020
)
Adipose tissue gene expression of CXCL10 and CXCL11 modulates inflammatory markers in obesity: implications for metabolic inflammation and insulin resistance
.
Therap. Adv. Endocrinol. Metab.
11
,
2042018820930902
43.
Kumaravel
S.
,
Singh
S.
,
Roy
S.
,
Venkatasamy
L.
,
White
T.K.
,
Sinha
S.
et al.
(
2020
)
CXCL11-CXCR3 axis mediates tumor lymphatic cross talk and inflammation-induced tumor, promoting pathways in head and neck cancers
.
Am. J. Pathol.
190
,
900
915
[PubMed]
This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY).

Supplementary data