Integrative Analyses of Whole-Transcriptome Sequencing Reveals CeRNA Regulatory Network of mRNAs, lncRNAs, miRNAs and circRNAs in Pulmonary Hypertension Treated with FGF21

Noncoding RNAs have been shown to play important roles in hypoxic pulmonary hypertension (HPH). Our preliminary data showed that HPH is attenuated by broblast growth factor 21 (FGF21) administration. Therefore, we further investigated the whole transcriptome RNA expression patterns and interactions in a mice HPH model treated with FGF21. By whole-transcriptome sequencing, differentially expressed mRNA, miRNA, lncRNA, and circRNA were successfully identied in normoxia (Nx) vs. hypoxia (Hx) and Hx vs. hypoxia + FGF21 (Hx + F21). Through intersection and predictive analysis, differentially co-expressed mRNA, miRNA, lncRNA, and circRNA were selected, followed by functional enrichment analysis. MAPK signaling pathway and epigenetic modication were enriched and may play fundamental roles in the therapeutic effects of FGF21. A ceRNA regulatory network was constructed with miR-7a-5p, miR-449c-5p, miR-676-3p and miR-674-3p as the core. Then the quantitative real time-PCR validation results were consistent with the results of whole-transcriptome sequencing. This study may provide potential biomarkers, pathway and ceRNA regulatory network in HPH treated with FGF21.


Introduction
As a subtype of pulmonary hypertension (PH) (group 3), hypoxic pulmonary hypertension (HPH) is an insidious and fatal disease with a poor long-term prognosis [1 2]. The reason for cure di culty in HPH patients is the lack of accurate and noninvasive biomarkers for early detection, and no effective medication with limited side effects [3]. Therefore, it is highly urgent to explore potential biomarkers and novel therapeutic targets for HPH.
Recently, high-throughput next-generation sequencing (NGS) technologies have enabled our research to be no longer limited to messenger RNA (mRNA) expression pro les, but from mRNA to non-coding RNA, such as microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA) [4 5]. With the deepening of research, non-coding RNAs have been found to play important roles in epigenetic regulation, transcriptional regulation, and posttranscriptional regulation, and participate in a variety of biological processes and diseases [6 7]. Competitive endogenous RNA (ceRNA) hypothesis mentioned by Salmena et al. states that mRNA, lncRNA, circRNA could act as miRNA sponge to modulate the stability or translational activity of target genes [8]. At present, the regulatory function of ceRNA has been widely accepted, and ceRNA pattern in HPH development has also been widely reported. As an example, Jiang Y et al. reported that circRNA Calm4 regulated hypoxia-induced pulmonary arterial smooth muscle cells (PASMCs) pyroptosis by ceRNA mechanism of sponging miR-124 [9]. However, few studies have delineated the transcriptomic landscape of HPH utilizing whole-transcriptome sequencing strategies.
Fibroblast growth factor 21 (FGF21) has been found to exert remarkable cardiovascular bene ts in recent years. Cong Liu et al. found that pharmacological treatment with FGF21 could reduce atherosclerosis by the improvement of plasma cholesterol metabolism [10]. Xuebo Pan et al. uncovered that replenishment of FGF21 alleviated angiotensin II-Induced hypertension and vascular dysfunction by angiotensinconverting enzyme 2 (ACE2) in mice [11]. Our preliminary data indicated that exogenous administration of FGF21 relieved HPH and right ventricular hypertrophy [12]. However, the underlying mechanism in PH treated with FGF21 is thought to be complex, and systematic elucidation of the mechanisms and e cacy of FGF21 in HPH is therefore necessary for wider acceptance and clinical application.
In this study, we focused on the perturbation of hypoxia intervention and FGF21 administration on the non-coding RNA panorama in HPH. We rst detected the expression of FGF21 in human and experimental HPH, and selected the appropriate dose of FGF21 treatment group mice (Hx + F21), together with the normoxia (Nx) and hypoxia mice (Hx), performed whole-transcriptome sequencing. Subsequently, differentially expressed mRNA, miRNA, lncRNA, and circRNA was screened between Nx and Hx groups, as well as Hx and Hx + F21 groups, followed by functional enrichment analysis. Through intersection and predictive analysis, we nally anchored the differentially co-expressed mRNA, miRNA, lncRNA, and circRNA regulated by hypoxia and FGF21, constructed a ceRNA regulatory network, and performed qRT-PCR veri cation. This research may provide potential biomarkers for the diagnosis and treatment of HPH, and provide a fundamental for clinical application of FGF21 in HPH.

Material And Methods
Human clinical blood samples 9 patients with high-altitude PH (HAPH) were enrolled between April 2019 and June 2019 at the First A liated Hospital of Wenzhou Medical University, and Qinghai Golmujianqiao Hospital. And 21 agematched healthy controls (HCs) were recruited. Blood samples were obtained from HAPH, or HCs in a quiet state. The blood was drawn into blood collection tubes, gently mixed and allowed to stand for 1 h at room temperature (RT). Then the tubes were centrifuged at 4°C within 2 h to separate the serum. The serum was immediately stored at -80°C until use to avoid degradation. All procedures performed in studies involving human participants were approved by the Ethics Committee of the First A liated Hospital of Wenzhou Medical University.

HPH models and drug applications
Male C57BL/6 mice were purchased from the Animal Center of the Chinese Academy of Sciences (Shanghai, China). All animals started modeling after a week of rest. The construction of HPH mouse model was described previously [13]. Brie y, mice were randomly divided into ve groups, each with 6 mice: normoxia (Nx, saline-treated), hypoxia (Hx, saline-treated), hypoxia plus FGF21 (Hx + F21, F21, 50, 100, 200 µg/kg/day). Mice in the hypoxic group were kept in a normobaric hypoxic animal chamber (9-11% oxygen concentration) for 3 weeks. The normoxia group was exposed to indoor air. All groups were given adequate water and food. All animal work was approved by the Animal Ethics Committee of Wenzhou Medical University. FGF21 was obtained from Prospec (CYT-281, Prospec, Rehovot, Israel). Brie y, FGF21 (1 mg) was dissolved in 10 ml 0·1% bovine serum albumin (BSA) solution. Then the medicine was injected intraperitoneally into mice at the above prescription dose per day before entering the hypoxia chamber.
After 3 weeks, all mice were anesthetized and subjected to invasive hemodynamic measurements; then, the mouse blood, and lungs were collected for further study.

Biochemical parameters
Mouse serum FGF21, and human serum FGF21 levels were detected with an FGF21 Quantikine Elisa kit (R&D Systems, Minneapolis, MN, USA), and a Human FGF21 immunoassay kit (ImmunoDiagnostics, lnc.), respectively. Blood glucose was measured using ACCU-CHEK active blood glucose meter.

Measurements of invasive hemodynamic and RV hypertrophy
Invasive hemodynamics measurements were measured as described [13]. We measured the mean right ventricular pressure (mRVP), and the mean carotid arterial pressure (mCAP) to evaluate the effectiveness of the model using PowerLab 8/35 multichannel biological signal recording system (AD Instruments, AUS). And the right ventricular wall, left ventricular wall and ventricular septum were separated to assess the hypertrophy of the right ventricle, and respond with two indicators: the ratio of right ventricle weight to body weight (RV/BW), and the ratio of right ventricle weight to left ventricle plus interventricular septum weight (RV/LV+S).
Tissue staining 4-µm-thick para n-embedded lung tissue blocks were sectioned for staining experiment. Hematoxylin-Eosin (HE) staining and Masson staining was both performed according to routine protocols. In bright eld, we use a microscope (Nikon, Tokyo, Japan) to select typical pulmonary arterioles (outer diameter 25-100 mm). The ratios of the wall thickness to the total thickness (WT/TT) and the ratios of the pulmonary artery wall area to the total area (WA/TA) were calculated to re ect pulmonary arterial remodeling. The ratio of the collagen area to the total area of the recording area was calculated to re ect the degree of collagen deposition. Immunohistochemical detection were used to determine FGF21 expression. Lung tissue sections were treated with citrate antigen repair buffer (pH 6·0) for antigen repair, then incubated with 0·3% hydrogen peroxide, and 5% BSA for blocking. Then, lung tissues were incubated with anti-FGF21 antibody (diluted 1:200, ab171941, Abcam) overnight at 4°C and incubated for 1 h with horseradish peroxidase-conjugated goat anti-rabbit IgG (Zhong Shan Jin Qiao) next day. Subsequently, the lung sections were observed under Nikon microscope and brown-yellow particles were regarded as positive signals of the immunohistochemical staining.

Western blot detection
A Pierce BCA protein assay kit (Thermo Fisher Scienti c Inc., Rockford, IL, USA) was used to quantify the total protein of mouse lung tissue homogenate or the total protein of rat pulmonary artery smooth muscle cells. After heat denaturation, it was prepared into 60-100 ug/ul protein standard and stored at 4°C. The protein was fully separated by 10% SDS-polyacrylamide gel electrophoresis (SDS-PAGE), and transferred to a 0.22 um pore PVDF membrane, and then incubated in 5% skim milk at 37°C for 90 min. After washing and cutting, the PVDG membrane was incubated in the primary antibody overnight at 4°C. the secondary antibody of rabbit or mouse with goat antiserum labeled with horseradish peroxidase (1: 10000, biosharp) was used to incubate at 37°C for 60 minutes, and then using ECL to exposure by (biorad exposure machine). The primary antibodies used were as follows: goat anti-rabbit lgG H&L (1:10000, ab205718, Abcam), rabbit anti-FGF21 (1:1000, ab171941, Abcam), rabbit anti-PCNA (1:1000, #13110, CST), rabbit anti-β-ACTIN (1:1000, #4970, CST). Experiments were repeated at least 3 times and β-ACTIN was used as an internal control.
For mRNA, lncRNA and circRNAs, total RNA was extracted from the lung tissue or rPASMCs using Trizol (Invitrogen), and chloroform. For miRNAs, total RNA was extracted from the lung tissue using SanPrep Column microRNA Extraction Kit. For mRNA, lncRNA and circRNAs, relative cDNA was synthesized from 1 to 3 mg of RNA using the iScript cDNA Synthesis Kit (Bio-Rad, USA). For miRNAs, relative cDNA was synthesized using miRNA First Strand cDNA Synthesis (Stem-loop Method) (Sangon Biotech, Shanghai, China). The mRNAs, lncRNAs and circRNAs expression were measured using qRT-PCR (CFX96 Real-Time System, Bio-Rad, USA). And the PCR for miRNAs was using microRNAs qPCR kit (Sangon Biotech, Shanghai, China). β-actin or U6 served as endogenous controls. The relative fold change was calculated using log2 FC. The primer sequences were shown in Table 1.

RNA sequencing
The whole transcriptome sequencing service of lung tissue of HPH animal models (5 mice in each group) were provided by RiboBio Co., Ltd (Guangzhou, China). Total RNA was extracted and prepared from lung tissue with Trizol reagent (Thermo Fisher Scienti c, USA), and its purity and concentration were determined by ND-1000 Nanodrop. All samples meeting the A260/A280>1.8 and A260/A230>2.0 standards were used for subsequent analysis. RNA integrity was evaluated on the Agilent 2200 TapeStation, and those with an average RNA integrity number (RIN) greater than 7.0 were used for sequencing. The NEBNext® Ultra™ RNA library was prepared using Illumina kits. RNA sequencing was performed by Illumina HiSeq 3000, using 2 × 150bp pair-end sequencing. The quality control of the original sequencing data was checked and processed by FASTQC software and Trimmomatic Tools. After comparing the sequencing data with the mouse reference genome mm10 using HISAT2 software, the reads mapped to each gene were calculated by HTSEQ. The raw data of RNA-seq was deposited in the NCBI Sequence Read Archive (SRA) under accession numbers, SRR SRR17180890-SRR17180934.

Differential expression analysis
For the whole transcriptome RNA sequencing results of three different treatments of mice, the "edgeR" R package [14] (3.34.1) was used for pairwise differential expression analysis (Nx vs. Hx; Hx vs. Hx + F21). All transcripts that meet the p.value <0.05, |log2 FC| >1 screening threshold are considered signi cant. Through the "ggplot" package, the volcano map visually displays dif-mRNAs, dif-miRNAs, dif-lncRNAs and dif-circRNAs, respectively. The "heatmap" R package was used to draw cluster heatmaps of the top 100 RNAs in |log2 FC| among the four transcripts in three different treatment groups.

Enrichment analysis
In this article, all gene enrichment analysis was performed by the "clusterpro ler" R package (v3.16.1).
Gene set enrichment analysis (GSEA) is a powerful method to enrich gene expression pro les based on the ranking of differential expression results to gain in-depth understanding of biological mechanisms. We download "c2.cp.kegg.v7.4.symbols.gmt" and "h.all.v7.4.symbols.gmt" from MSigDB as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and hallmark gene reference set. The log2 FC obtained by the differential expression analysis of Nx vs Hx and Hx vs Hx + F21 was used as the ranking list for GSEA analysis. All results of GSEA were visualized by "enrichplot" R package.
dif-genes that were up-regulated in Hx group compared with Nx group, and down-regulated in Hx + F21 group compared with Hx group were recognized as co-expressed dif-genes. Conversely, dif-genes with another co-regulatory relationship were also selected as co-expressing dif-genes. The co-expressed difgenes obtained above were used for KEGG and Gene Ontology (GO) enrichment analysis. Before analysis, org.Mm.eg.db(http://www.bioconductor.org/packages/release/data/annotation/html/org.Mm.eg.db.html, 3.14.0) was used to annotate and map dif-genes. The bubble charts drawn by the "ggplot2" R package were used to display the enrichment results of KEGG and GO, among which GO result displayed the top 10 enriched items in biological process (BP), cellular component (CC), and molecular function (MF). And p.value <0.05 was considered signi cant.
After the target genes predicted by the three databases were combined, they were further integrated with co-expressed dif-genes to obtain the regulatory relationship of dif-miRNA-dif-genes.
Finally, we focused on screening co-expressed dif-miRNAs which can simultaneously regulate dif-lncRNAs, dif-circRNAs, and dif-genes for subsequent construction of ceRNA networks. It is worth noting that |log2 FC|>1 was too harsh for the RNA-Seq results of circRNA, so we relaxed the standard for circRNA when screening dif-circRNA-dif-miRNA regulatory relationships. We only cared about p.value<0.05, and the threshold of |log2 FC| was not required [19]. The regulatory relationship will be experimentally veri ed in the follow-up qRT-PCR.

Enrichment analysis of target genes regulated by coexpression dif-miRNAs
In order to further explore the biological role of these co-expressed dif-miRNAs in the process of FGF21 treatment of HPH, we performed KEGG and GO enrichment analysis on the target genes they regulate. As before, the "clusterpro ler" R package was responsible for enrichment analysis, and the "ggplot2" R package was responsible for visualization, p.value <0.05 was considered as signi cant enrichment.
Construction of the dif-circRNA/dif-lncRNA-dif-miRNA-difgenes network To investigate whole transcriptome ceRNA regulation in HPH, the dif-circRNA/dif-lncRNA-dif-miRNA-difgenes ceRNA network was built based on the miRNA sponge hypothesis [20], which proposed that lncRNAs and circRNAs could act as miRNA sponges to interact with miRNAs and regulate downstream mRNAs expression. Thus, on the basis of previous selection, we further screen dif-lncRNAs/dif-circRNAs and dif-genes that have the same variation trend, their corresponding miRNAs which have the opposite trend were chosen simultaneously. RNAs that ful lled the above condition were combined to construct a dif-circRNA/dif-lncRNA-dif-miRNA-dif-genes ceRNA network in Nx vs. Hx and Hx vs. Hx + F21, respectively. In order to further re ect the therapeutic effect of the ceRNA network regulated by FGF21 on PH, we screened the ceRNA network relationships that exist simultaneously in Nx vs. Hx and Hx vs. Hx + F21. And these networks were integrated and constructed to obtain the nal co-expressed ceRNA network of the whole transcriptome. Cytoscape v3.7.1 [21] was applied to display the network results.

Statistical analysis
All Statistical analyses were performed with GraphPad Prism 9·0. All data were expressed as mean ± standard error mean (SEM). Comparisons between two groups were analyzed by unpaired two-tailed Student's t-test. Multiple comparisons were analyzed by ANOVA followed by Dunnett post-hoc test.
P.value of < 0.05 were considered statistically signi cant.

Result
Expression of FGF21 is decreased in human and experimental PH.
In this study, hypoxia-induced rPASMCs, and mouse PH models were established to determine the pathophysiological role of FGF21 in PH. Signi cantly lower FGF21 mRNA, and protein levels in rPASMC suspensions, and lung homogenates were detected after exposure to hypoxia ( Figure. 1A-D). FGF21 was diffusely distributed, mainly in the vascular cell periphery (Figure. 1E). Circulating FGF21 levels in HPH mice, and HAPH patients were further examined. Consistent with the decline of serum FGF21 levels in HPH mice, circulating FGF21 levels were notably decreased in subjects with HAPH ( Figure. 1F-G). These data support ndings from our previous research [13 22] and demonstrate that FGF21 is downregulated in HAPH patients and plays a role in PH disease.
PH amelioration in FGF21 treated mice.
To further investigate the role of FGF21 in HPH, the optimal dose of human recombinant FGF21 was determined by evaluating its effect in HPH mice administered different concentrations (50, 100, and 200 µg/kg, intraperitoneally). In hypoxia-exposed mice subjected to daily administration of FGF21 at a dose of 200 µg/kg, mRVP decreased signi cantly whereas the mCAP was unaffected, compared with that in hypoxia-exposed vehicle-treated mice ( Figure. Figure S3, DATA S5, S6, S7, S8). We then performed the hierarchical clustering heatmaps on the top 100 dif-mRNAs, dif-miRNAs, dif-lncRNAs, and dif-circRNAs from all three groups as shown in Figure 3. The results showed that Hx samples can be signi cantly separated from Nx and Hx + F21 samples, indicating the analysis results were reliable.
GSEA of differential expressed genes in Nx vs. Hx and Hx vs. Hx + F21.
GSEA was performed using the Hallmark and KEGG gene sets to explore the biological functions of dysregulated mRNAs. Comparing Hx with Nx, GSEA showed that a total of ve KEGG gene sets ( Figure   4A) and eight hallmark gene sets ( Figure 4C) were signi cantly enriched. Comparing Hx + F21 with Hx, GSEA showed that a total of three KEGG gene sets ( Figure 4B) and four hallmark gene sets ( Figure 4D) were signi cantly enriched. Overall, KEGG gene sets cell cycle and pyrimidine metabolism, Hallmark gene sets E2F targets, G2M checkpoint and MYC targets were downregulated in Hx, compared with Nx while these gene sets were all upregulated in Hx + F21, compared with Hx ( Figure 4A-B). Moreover, KEGG gene sets endocytosis, and Hallmark gene sets estrogen response early were upregulated in Hx, compared with Nx while these gene sets were all downregulated in Hx + F21, compared with Hx ( Figure 4C-D).
Screening of differentially co-expressed molecules and functional enrichment analysis.
To investigate the common and unique differential RNAs among all samples, we constructed venn diagrams. The results showed that 47 dif-genes, 5 dif-miRNAs, 113 dif-lncRNAs, and 11 dif-circRNAs were downregulated in Hx, compared with Nx; these co-expressed RNAs were upregulated in Hx + F21, compared with Hx ( Figure 5A-D, Data S9, S10, S11, S12). Moreover, 321 dif-genes, 6 dif-miRNAs, 520 dif-lncRNAs, and 7 dif-circRNAs were upregulated in Hx, compared with Nx; these co-expressed RNAs were all downregulated in Hx + F21, compared with Hx ( Figure 5A-D, Data S9, S10, S11, S12). Then we perform KEGG and GO analysis on the co-expressed dif-genes. KEGG pathway analysis showed that KEGG pathway such as cardiac muscle contraction, adrenergic signaling in cardiomyocytes and lysine degradation may be in volved in the effect of FGF21 in alleviating PH ( Figure 5E). Moreover, GO terms such as covalent chromatin modi cation, histone modi cation and peptidyl-lysine modi cation of BP; nuclear chromatin, transcription regulator complex and transmembrane transporter complex of CC; transcription coregulator activity, RNA polymerase II-speci c and DNA-binding transcription activator activity of MF were potentially associated with therapeutic effect of FGF21 ( Figure 5F), which indicating that FGF21 may have a potential role in modulating transcriptional activity and the epigenetic landscape.
MAPK signaling pathway may play an important role in HPH treated with FGF21.
The target genes for 11 co-expression dif-miRNAs (5 down and 6 up, in Hx, compared with Nx or Hx + F21) came from the intersection of three database prediction results and co-expressed dif-genes among three groups as shown in Figure 6A-B and Data S13. Comparing Hx with Nx, there were 243 target genes.
KEGG pathway analysis on these genes showed that MAPK signaling pathway was the most signi cantly enriched pathway ( Figure 6C). GO analysis revealed that regulation of ion transmembrane transport, positive regulation of cell projection organization and lung epithelial cell differentiation were enriched ( Figure 6D). Comparing Hx + F21 with Hx, there were 215 target genes. MAPK signaling pathway was also the most signi cant enriched pathway from KEGG analysis on these target genes ( Figure 6E). As for GO BP terms, histone modi cation, covalent chromatin modi cation and regulation of transporter activity were signi cantly enriched ( Figure 6F). These observations support the enrichment analysis results of the co-expressed dif-genes as shown in Figure 5E-F.
Construction of the ceRNA network in Nx vs. Hx and Hx vs. Hx + F21.
Construction of the co-expressed ceRNA network.
To explore the co-expressed ceRNA network regulated by FGF21, we next took the intersection of the upregulated-dif-miRNA-associated ceRNA network in Hx vs. Nx, and the downregulated-dif-miRNAassociated ceRNA network in Hx + F21 vs. Hx. There were 2 dif-miRNAs, 16 dif-genes, 16 dif-lncRNAs, and 2 dif-circRNAs to be found. Also, the downregulated-dif-miRNA-associated ceRNA network in Hx vs. Nx, and the upregulated-dif-miRNA-associated ceRNA network in Hx + F21 vs. Hx were taken to intersect. There were 2 dif-miRNAs, 79 dif-genes, 113 dif-lncRNAs, and 2 dif-circRNAs to be found. Cytoscape was applied to build the co-expressed ceRNA network (Figure 9, Data S18). These integrated core ceRNA network revealed potential co-expression relationships of genes, lncRNAs, and circRNAs regulated by these miRNAs in PH treated with FGF21.
QRT-PCR validation of differentially expressed RNAs in the core ceRNA network.

Discussion
FGF21 is a well-established vascular protective factor [10 23 24], and changes in circulating FGF21 levels are closely related to many diseases, such as subclinical atherosclerosis [25], alcoholic cardiomyopathy [26], and remyelination in the central nervous system [27]. Our research showed that the reduction of FGF21 level in the serum of HAPH patients, as well as in the serum and lung tissue of mice HPH model, which suggesting a strong correlation between the levels of circulating and lung FGF21 and HPH.
Administration of FGF21 alleviated the progression of HPH in our research, which is consistent with the vascular bene t of FGF21 reported by most researchers, and also broadens the application of FGF21 in the pulmonary circulation.
However, the mechanism underlying the therapeutic effect of FGF21 is complicated. We assumed that epigenetics may play an important role in its anti-vascular remodeling properties. By using wholetranscriptome sequencing, we identi ed a total of 1396 dif-mRNAs, 28 dif-miRNAs, 1263 dif-lncRNAs, and 48 dif-circRNAs in Hx compared with Nx, and a total of 1348 dif-mRNAs, 20 dif-miRNAs, 1027 dif-lncRNAs, and 13 dif-circRNAs in Hx + F21 compared with Hx. Enrichment analysis on the dif-mRNAs showed that cell cycle and G2M checkpoint might participate in the regulation of hypoxia, and the protective effect of FGF21 under hypoxic stimulation might be related to its mediation of cell endocytosis.
Shuang-Lan Xu et al. [28] and Jie Liu et al. [29] identi ed the differentially expressed circRNA/miRNA/mRNA and lncRNA/miRNA/mRNA in rat HPH model compared with controls, respectively. Their target genes were enriched in immune and in ammation-related pathways. Quite different from the rat model, in this study, we focus on mice models. Compared with immune and in ammatory effects, cell cycle may be more related to the pathology of vascular remodeling of HPH [30 31]. In addition, we used whole-transcriptome sequencing to construct a dif-circRNA/dif-lncRNA-dif-miRNA-dif-genes interaction ceRNA network landscape, which more comprehensively demonstrated the relationship between these transcripts in the ceRNA network of HPH mice, also provided potential biomarkers for predicting new therapeutic targets for HPH.
In addition to the disturbance of hypoxia, we paid more attention to the in uence of FGF21. Through the intersection of the three groups, we focused on identifying the differentially expressed RNA molecules and enriched signaling pathways of their target genes during the course of disease development and FGF21 treatment. The result of miRNA target gene function enrichment showed that the chromatin modi cation, histone modi cation and peptidyl-lysine modi cation was affected after FGF21 intervention, which veri ed our hypothesis that epigenetics might be involved in FGF21 regulation. In this regard, Sangwon Byun et al. [32] reported that FGF21 promoted Jumonji-D3 (JMJD3/KDM6B) histone demethylase, thereby activated hepatic autophagy and lipid degradation. Teresa Płatek et al. [33] showed that DNA methylation changes and miRNA pattern play an important role in FGF21 resistance in obesity.
Moreover, our data showed that MAPK signaling was involved in the process of HPH treated by FGF21, and the relationship between FGF21 and MAPK signaling has been reported in many articles [34][35][36].
In the ceRNA network constructed in our study, miRNAs are the bridges as well as the core factor. For the two RNAs upregulated by hypoxia, miRNA-674 could be induced by HIF-1α [37] and was upregulated in chronic obstructive pulmonary disease rats [38], which were closely related to HPH; miR-676 was suggested as a potential therapeutic target in heart failure [39], suggesting that miR-676 might affect ventricular remodeling. These all support our observations. For the other two RNAs downregulated by hypoxia, miR-449 suppressed cell proliferation and promoted apoptosis [40], suggesting that FGF21induced upregulation of miR-449 expression might promote the homeostasis of pulmonary vascular structure through similar regulation; the overexpression of miR-7-5p suppressed the growth and cancer metabolism of lung cancer [41], and HPH is characterized by enhanced proliferation of vascular smooth muscle cell that shares similar pathophysiological processes with cancer cell. Due to space limitations, we cannot discuss the RNA molecular in core ceRNA one by one. Subsequent qRT-PCR veri cation showed that our analysis was reliable. Future studies are warranted to ascertain further this association.
In conclusion, this study explored the potential molecular mechanism of HPH and the FGF21-treated HPH using whole-transcriptome sequencing. Epigenetic modi cation and MAPK signaling may play fundamental roles in the pathogenesis of HPH, and may account for the therapeutic effect of FGF21.
FGF21 may reduce hypoxia-induced pathological deterioration through a ceRNA network with miR-7a, miR-449, miR-676 and miR-674 as the core. These ndings suggest potential new targets in HPH characterization and FGF21 therapy. However, this study is limited by the small sample size, and further in vitro and in vivo veri cation is required to delineate the exact mechanism details.

Declarations Data availability statement
The raw data presented in this study will be made available by the authors, without undue reservation.
Sequencing data used in this manuscript can be found in the NCBI database SRA (PRJNA785350) and the SRR accession numbers are SRR17180890 to SRR17180934.

Declaration of interests
The authors have declared that no competing interest exists. Figure 1 Expression of FGF21 is decreased in human and experimental PH. Veri cation of the mRNA (A) and protein (B) levels of FGF21 in rPASMCs under normoxia and hypoxia by qRT-PCR and western blotting (WB), respectively. The mRNA (C) and protein (D) levels of FGF21 in lung tissue homogenates of adult C57BL/6 mice under normoxia and hypoxia were veri ed by qRT-PCR and WB. The localization and quanti cation of FGF21 in lung para n sections of adult C57BL/6 mice under normoxia and hypoxia (n = 6, respectively) are shown by immunohistochemical images (E). ELISA was used to detect the serum FGF21 expression level of adult C57BL/6 mice (F) under normoxia and hypoxia or the serum FGF21 expression level of healthy individuals (n = 21) and HAPH patients (n = 9) (G). *p < 0.05, **p < 0.01, ***p < 0.001. Scale bars, 20 μm. show the degree of collagen deposition around the pulmonary arteries. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Scale bars, 20 μm.      The competing endogenous RNA network in Nx vs. Hx. Triangles represent miRNAs, red indicates upregulation, green indicates down-regulation; squares represent genes, pink indicates up-regulation, and light green indicates down-regulation; diamonds represent lncRNAs, cyan indicates up-regulation, and lavender indicates down-regulation; rounds indicate circRNAs, purple indicates up-regulation and dark blue prompts a downward adjustment.  The competing endogenous RNA network in Nx vs. Hx vs.Hx + F21. The red triangles represent miRNAs, the pink squares represent genes, the cyan diamonds represent lncRNAs, and the purple circles represent circRNAs.