Pulmonary Hypertension Remodels the Genomic Fabrics of Major Functional Pathways

Pulmonary hypertension (PH) is a serious disorder with high morbidity and mortality rate. We analyzed the right-ventricular systolic pressure (RVSP), right-ventricular hypertrophy (RVH), lung histology, and transcriptomes of six-week-old male rats with PH induced by (1) hypoxia (HO), (2) administration of monocrotaline (CM), or (3) administration of monocrotaline and exposure to hypoxia (HM). The results in PH rats were compared to those in control rats (CO). After four weeks exposure, increased RVSP and RVH, pulmonary arterial wall thickening, and alteration of the lung transcriptome were observed in all PH groups. The HM group exhibited the largest alterations, as well as neointimal lesions and obliteration of the lumen in small arteries. We found that PH increased the expression of caveolin1, matrix metallopeptidase 2, and numerous inflammatory and cell proliferation genes. The cell cycle, vascular smooth muscle contraction, and oxidative phosphorylation pathways, as well as their interplay, were largely perturbed. Our results also suggest that the upregulated Rhoa (Ras homolog family member A) mediates its action through expression coordination with several ATPases. The upregulation of antioxidant genes and the extensive mitochondrial damage observed, especially in the HM group, indicate metabolic shift toward aerobic glycolysis.


Introduction
Pulmonary hypertension (PH), characterized by increasing pulmonary artery pressure and vascular resistance [1], is a serious complication of a number of unrelated cardiopulmonary, inflammatory, and autoimmune diseases, as well as of drug toxicity. PH is characterized by endothelial dysfunction, enhanced vasoconstrictor reactivity, activation of proliferative and anti-apoptotic pathways, vascular remodeling, elevated pulmonary artery pressure, and right-ventricular hypertrophy (RVH), leading to right-ventricular failure and premature death.
Recently updated PH classification maintains the original five groups with some modification [2]. The term pulmonary arterial hypertension (PAH) is assigned to Group1. It includes idiopathic (IPAH) and heritable PAH (HPAH), PAH associated with congenital heart diseases, drug toxicity, connective tissue disorders, portal hypertension, human immunodeficiency virus (HIV) infection, and schistosomiasis. The average survival time in patients with PAH without treatment was reported in 1991 to be around 2.8 years [3]. Despite the improvement in quality of life, the current therapy fails to reverse or halt the progression of the vascular disease [4]. Progressive pulmonary vascular changes lead to neointimal lesions resulting in irreversibility of the PH [5]. The five-year survival rate remains As previously described [36], rats were anesthetized with an intraperitoneal injection of xylazine (6 mg/kg) and ketamine (60 mg/kg). Through an incision in the neck, the trachea was exposed and cannulated with PE 240 tubing, and the rat ventilated in room air (~70-80 breaths/min). The chest was opened, PE 50 tubing with a small needle was inserted into the right ventricle, and the pressure was recorded on a Grass polygraph (model 7E). At the end of the pressure measurements, the lungs were perfused with autoclaved normal saline to remove blood. The left lung and the heart were placed in 10% buffered formaldehyde. The right lung was quickly frozen in liquid nitrogen and stored at −80 • C for transcriptomic analysis at a later date.

Assessment of Right-Ventricular Hypertrophy (RVH)
One week later, the heart was removed from formaldehyde and the atria were trimmed. The free wall of the right ventricle (RV) was separated from the left ventricle (LV) and the septum (S) and weighed. The RVH was calculated for all groups as follows: Right-ventricular hypertrophy occurs when the average RV/(LV + S)(PH)) is significantly (p < 0.05) larger than RV/(LV + S)(CO) (i.e., RVH significantly larger than one).

Histology
Formalin-preserved lung tissue was processed for paraffin block. Then, 5-6-µm sections were cut and stained with hematoxylin and eosin (H&E) for evaluation of the pulmonary vasculature.

Microarray
Lungs were quickly removed, frozen in liquid nitrogen, and stored at −80 • C. Our optimized protocol [37,38] and the "multiple yellow" strategy [39] were used to profile each "condition" in four biological replicas. Briefly. total RNA was extracted with Qiagen RNeasy minikit, concentration was determined with a NanoDrop ND-2000 Spectrophotometer, and purity determined with an Agilent RNA 6000 Nano kit in an Agilent 2100 Bioanalyzer. Total RNA was then reverse-transcribed in the presence of Cy3/Cy5 dUTP and the incorporation of fluorescent tags was determined again with the nanodrop. Differently labeled samples of biological replicates were hybridized overnight with Agilent 60 mer 4 × 44k whole genome rat V2 arrays (#G2519F). The arrays were scanned with an Agilent G2539A dual-laser scanner, and primary analysis was performed with (Agilent) Feature Extraction 11.1 software. All spots affected by local corruption or with foreground fluorescence less than twice the background were disregarded, and data were normalized and filtered following our standard procedure [40].

Transcriptomic Analysis
A gene was considered as significantly regulated if the absolute fold-change exceeded the cut-off calculated for that gene [41]. The cut-off (CUT, Equation (1)) accounts for the observed expression variability of that gene in the compared conditions and the existence of several spots probing it redundantly in the microarray. The non-uniform expression variability among the genes results from both the biological variability within biological replicas and the inherent technical noise. This strategy aims to substantially reduce both the positive and the negative false hits when using an arbitrary cut-off (e.g., 1.5×) for the fold-change. CUT (Equation (1)) adds correction for the groups of microarray spots redundantly probing the same gene.
, where : (PH) = (HO), (CM), (HM) r i = 4R i − 1 = number of degrees of freedom, R i = number of spots probing the same gene i χ 2 = chi − square score for r i degrees of freedom and probability ε (= 0.05) in this report s (condition) ik = standard deviation of gene i probed by spot k in the specified condition µ (condition) ik = average expression of gene i probed by spot k in the specified condition (2) Profiling four biological replicates allows independent assessment of (i) average expression level, (ii) expression variability, and (iii) expression coordination of each gene with each individual gene. Combination of these measures was used to determine the remodeling of the topology and interplay of the functional genomic fabrics. The genomic fabric was defined by us as the most inter-coordinated and stably expressed gene network responsible for selected biological processes [42].

Pathway Analysis
The genes pertaining to important functional pathways were selected from the maps developed by Kyoto Encyclopedia of Genes and Genomes (KEGG) [43]. Particular attention was given to the genes involved in oxidative phosphorylation (map 00190, ATPases, cytochrome c oxidases, NADH dehydrogenases), chemokine signaling pathway (map 04062), vascular smooth muscle contraction (map 04270), mitochondrion ribosomal proteins, and cell cycle (map 04110).
The corrected weighted pathway regulation (WPR) [44] was used to quantify the transcriptomic effects of the three PH models on the selected functional pathways.
WPR is more informative than the traditional percentage of regulated genes. WPR weights the contribution of each gene proportional to its normal (here, in the CO group) average expression, absolute fold-change (until the cut-off accounting for the combined contributions of the technical noise and biological variability), and confidence of the expression regulation.
Like in previous experiments with kidneys [45] and hearts of mice subjected to chronic hypoxia [39,46,47], we computed the pair-wise Pearson correlation coefficient ρij between the expression levels of all pairs of expressed genes i and j in biological replicas of each condition. As speculated in previous papers [48,49], strong expression correlation of two genes may represent a kind of "transcriptomic stoichiometry" aiming to optimize the pathway involving the encoded proteins. Expressions of (p < 0.05) statistically significantly synergistically expressed genes (ρij > 0.950) oscillate in phase (both genes increase and decrease their expression simultaneously), while expressions of antagonistically expressed genes (ρij < −0.950) manifest opposite tendencies (when one increases, the other decreases). Genes with |ρij| < 0.025 are considered independently expressed. Different from the traditional cluster analysis that correlates the genes according to their similar behavior across various conditions or time points, we compute the correlation only between the expression levels in biological replicas of the same condition to determine "stoichiometric gene networks" [49].
The influence of a given gene i toward a particular pathway Γ was evaluated by its average coordination power (CP) against the expressions of the pathway Γ genes.
CP takes values from −100% (all Γ-genes perfectly negatively correlated with gene i) to 100% (all Γ-genes perfectly positively correlated with gene i). A neutral gene toward that pathway has the CP close to zero, indicating that either almost all correlations are null or the positives are balanced by the negatives.
Here, we introduce a new measure, termed gene pair prominence (GPP), to quantify the relevance of gene interplays in each condition. The two genes can be from the same or from two different functional pathways. GPP takes values from −100% for the most prominent, negatively correlated genes to 100% for the most prominent, positively correlated genes in that condition. In order to simplify the landscape, GPP was considered zero for gene pairs whose absolute Pearson product-moment correlation coefficient between their expression levels in the four biological replicas of that condition was not significant (|ρ ij | < 0.950).

Decreased Weight Gain
The four-week weight gain in the control group was (204 ± 23) g. All PH groups had significantly lower weight gains as presented in Figure 1a; the lowest was observed in MCT-treated rats exposed to hypoxia, i.e., group HM (WG = 33.6 ± 16.0) g.

Decreased Weight Gain
The four-week weight gain in the control gou

Right-Ventricle Hypertrophy
Figure 1b presents the significant increase in ratio between the right-ventricle weight and the sum of the left-ventricle and septum weights in the three PH models as normalized to the RV/(LV + S) in the control group. Note that right-ventricle hypertrophy was statistically significant in all three experimental models.

Increased Right-Ventricle Systolic Pressure
As depicted in Figure 1c, RVSP was significantly increased in all PH groups (HO, CM, HM) compared to the controls (CO). While the difference between HO and CM groups was not statistically significant (p = 0.364), the differences between the HM and HO groups (p = 0.000127) and between HM and CM groups (p = 0.00078) were highly significant.  Figure 1b presents the significant increase in ratio between the right-ventricle weight and the sum of the left-ventricle and septum weights in the three PH models as normalized to the RV/(LV + S) in the control group. Note that right-ventricle hypertrophy was statistically significant in all three experimental models.

Increased Right-Ventricle Systolic Pressure
As depicted in Figure 1c, RVSP was significantly increased in all PH groups (HO, CM, HM) compared to the controls (CO). While the difference between HO and CM groups was not statistically significant (p = 0.364), the differences between the HM and HO groups (p = 0.000127) and between HM and CM groups (p = 0.00078) were highly significant.

Histological Alterations
Figure 1d-i present H&E-stained pulmonary arteries (external diameter, 42-101 µm) from the controls (D) and the experimental groups: HO (e), CM (f), and HM (g, h, i). Pulmonary arteries from CM and HO show substantial medial wall thickening compared with the control. The arteries from the HM group show further medial wall thickening, and as seen in (h) and (i), exhibit neointima formation and occlusion of the lumen.

Overview of the Transcriptomic Alterations
Raw and processed gene expression results were deposited and are publicly available at [50]. Figure 2a,b indicate that all three experimental PH conditions regulated large numbers of genes and produced large weighted pathway regulation (WPR) scores. The overall transcriptomic alterations (% of ALL REG and WPR) were consistent with the alterations in weight gain, RVH, and RVSP data. Thus, MCT alone had a slightly larger effect than hypoxia alone; when MCT-treated rats were exposed to hypoxia, the alterations were substantially larger. Figure 2c presents the WPR scores of some interesting groups of genes involved in respiration. Figure 2d presents the significant fold-changes of genes whose regulation is often associated with PH. Remarkably, in all three rat models of PH, Cav1 was overexpressed, by 13.18× in HO, 39.14× in CM, and 79.06× in HM, with a substantial increase in the HM group. We found also that upregulation of Cav1 was accompanied by the upregulation of Nos3 (by 4.21× in HM) and upregulation of Mmp2 (by 9.14× in HO, 33.59× in CM, and 81.47× in HM group). The vast majority of the regulated genes by separate exposure to either hypoxia or MCT were altered in the same direction. However, as presented in Figure 2e, few of them were oppositely regulated. Interestingly, when hypoxia and MCT were combined (as in HM), the regulation mostly followed the effects of the MCT exposure alone but with a larger fold-change (except for Plcb4 whose regulation in the HM group was opposite to regulations in both HO and CM groups).

Regulation of the Immune-Inflammatory Response
We found that PH regulated the expression of numerous genes involved in the (KEGG-determined map 04062) chemokine signaling pathway (Table 1). Table S1 in the Supplementary presents other immune-inflammatory response genes (chemokines, cytokines, cytokine receptors, interferons, interleukins, and tumor necrosis factors) that were significantly regulated in at least one of the three rat PH models. Importantly, compared to the traditional standard of uniform 1.5× absolute fold-change, our procedure to attach fold-change cut-offs to every single gene in each comparison identified additional regulated genes and eliminated false hits.

Alteration of Cellular Respiration
Particular attention was given to the regulation of the major groups of genes involved in the active ionic transport across the plasma membrane (ATPases) and oxidative phosphorylation: cytochrome c oxidases and NADH dehydrogenases ( Table 2). As illustrated In Figure 2c, the weighted pathway regulation (WPR) score revealed that these groups of genes were altered much more than the entire transcriptome. Remarkably, Atp1b2, one of the major players responsible for establishing and maintaining the transmembrane electrochemical gradients of Na + and K + , was the most downregulated in all three PH models (−146× in HO, −257× in CM, and −111× in HM). Table 2. Significantly up-(positive fold-change in bold numbers) and downregulated (negative fold-change in bold numbers) ATPases, cytochrome c oxidases, and NADH dehydrogenases in the three PH rat models. Note the range of the CUTs, from 1.24 (for the regulation of Ndufv3 in HM vs. CO) to 3.68 (for the regulation of Cox6a1 in HO vs. CO). Our procedure to determine the cut-off for the absolute fold-change for every gene in each comparison instead of using a fixed cut-off (like 1.5×) eliminated the false regulations (expression ratio in Italics, e.g., Atp11b in HO/CO), whose absolute (although over 1.5×) fold-change was below the CUT (2.28) computed for that gene in the compared conditions.   Alteration of the genes encoding mitochondrial ribosomal proteins is also presented for comparison in the Supplementary Table S2.

Remodeling of the Relationship between the Genes Involved in the Vascular Smooth Muscle Contraction and Oxidative Phosphorylation
We studied how much the interaction of the vascular smooth muscle contraction genes (VSMC, rno04270) (presented in Figure 3) with genes involved in oxidative phosphorylation (OP), particularly with ATPases, is remodeled by pulmonary hypertension. For this, we computed the Pearson correlation coefficients for the 1248 pairs that can be formed with the 32 ATPases and 39 VSMC genes adequately quantified in each of the four experimental conditions. Interestingly, we found that each PH condition significantly reduced the net positive correlation (i.e., number of positively (p) − number of negatively (n) correlated gene-pairs): from 450 (= 472p − 22n) in CO to 236 (= 280p − 44n) in HO, 340 (= 455p − 115n) in CM and 208 (= 262p − 54n) in HM. Figure 5 presents the significantly synergistically, antagonistically, and independently expressed pairs of VSMC genes with ATPases in CO, and the pairs whose expression correlation was significantly altered in the three PH models. Interestingly, the negative correlation of Rhoa with 15 ATPases (Atp11b, Atp13a2, Atp1a1, Atp1a3, Atp2a3, Atp2b1, Atp2c1, Atp4a, Atp5e, Atp5g3, Atp6v0a2, Atp6v0b, Atp6v0c, Atp6v0e1, Atp6v1g2) and no positive correlation in control rats was cancelled in the PH models, with several correlations even reversed.
Genes 2020, 11, 126 17 of 27 Figure 5. Synergistically (red squares), antagonistically (blue squares), and independently (yellow squares) expressed pairs of VSMC genes and ATPases in CO (a) whose expression correlation was significantly altered by hypoxia (b), monocrotaline (c), or the combined action of hypoxia and monocrotaline (d). A blank square in (a) indicates that the correlation is not (p < 0.05) statistically significant. Figure 5. Synergistically (red squares), antagonistically (blue squares), and independently (yellow squares) expressed pairs of VSMC genes and ATPases in CO (a) whose expression correlation was significantly altered by hypoxia (b), monocrotaline (c), or the combined action of hypoxia and monocrotaline (d). A blank square in (a) indicates that the correlation is not (p < 0.05) statistically significant.
We then used the new measure gene pair prominence (GPP, Equation (4)) to rank the gene pairs of the interplay between the two groups of genes. The positive and negative GPPs were plotted separately (as CO+ and CO−, HO+ and HO−, CM+ and CM−, HM+ and HM−) to emphasize the type and importance of their expression correlation. The GPPs of all genes in each condition were expressed as the percentage of the maximal absolute GPP in that condition, regardless of whether it was positive or negative. Except for the CM condition, where the genes of the dominant pair (Cyp4a8-Atp5e) were negatively correlated, the positively correlated pairs had much larger GPPs than the negatively correlated ones, meaning that expression of most ATPases oscillates in phase with most VSMC genes.
The analysis pointed out three remarkable H + transporters (Atp5e, Atp6v0a2, Atp6v0e2) in strong partnership with three VSMC genes: Cyp4a8 (cytochrome P450 family 4, subfamily a, polypeptide 8), Gnas (GNAS complex locus), and Mapk3 (mitogen-activated protein kinase 3). Mapk3, involved in a wide range of pathways (listed in [51]), provides the most prominent positive pairs with Atp6v0a2 in normal atmospheric conditions and with Atp6v0e2 in hypoxia (obtained by reducing the atmospheric pressure in half).

Discussion
We examined the transcriptomic alterations in the right lungs of male rats with pulmonary hypertension induced by exposure to hypoxia (HO group), administration of monocrotaline (CM group), or exposure to hypoxia of monocrotaline-treated animals (HM group). There are numerous common features of the three experimental groups, making them valuable models for the human pulmonary hypertension, including lower weight gain but right-ventricle hypertrophy and increased right-ventricle systolic pressure. However, there are also significant differences as revealed by the present and previous studies.
At four weeks after hypobaric hypoxia, there is significant PH and RVH. The hypoxia-induced PH in rats is not associated with any disruption of EC, loss of endothelial Cav1 or Nos3, or enhanced expression of Cav1 in VSMC. Despite the presence of endothelial Cav1 protein, pro-proliferative and anti-apoptotic pathways are activated, indicating Cav1 dysfunction [30,52]. In addition, the lung sections from infants with respiratory distress syndrome show that the presence of increased pulmonary artery pressure unaccompanied by endothelial damage does not result in the loss of endothelial CAV1 or enhanced expression of CAV1 in VSMC [29].
In the CM model, the loss of endothelial Cav1 is accompanied by the loss of endothelial proteins such as Pecam1 and soluble guanylate cyclase, and the activation of pro-proliferative and anti-apoptotic pathways leading to PH. At two weeks post monocrotaline administration, progressive loss of endothelial Cav1 is associated with PH and RVH in rats [12,40]. As the disease progresses, a further loss of endothelial Cav1 is accompanied by a loss of von Willebrand factor (vWF), indicative of extensive endothelial damage. Importantly, some of the arteries with vWF loss start to exhibit enhanced expression of CAV1 in VSMC, accompanied by increased Mmp2 expression and activity [51]. At four weeks, the significant loss of endothelial Cav1 is accompanied by enhanced expression of Cav1 in VSMC. The total Cav1 expression in the lungs increases by 39.14× compared to controls.
Exposing monocrotaline-treated rats to hypoxia (HM group) accelerates the disease process, resulting in significantly higher RVSP and RVH, extensive EC damage, and endothelial Cav1 loss, accompanied by enhanced expression of Cav1 in VSMC and neointima lesions [36,53,54]. Furthermore, 61% of the arteries exhibited enhanced expression of Cav1 in VSMC; Cav1 expression in the lungs increased to 81% of the controls. However, the neointimal lesions in the HM group exhibited loss of Cav1 and normal Nos3 expression. In the absence of Cav1, Nos3 gets uncoupled and produces nitrosative and oxidative stress. Cav1 expression is significantly increased (79.06×). It is worthy of note here that the neointimal cells have low Cav1 and near normal Nos3 expression, which is a set-up for oxidative and nitrosative stress. These results strongly support a dual role of Cav1 in the pathogenesis and the progression of PH, similar to what was reported in cancer [55,56].
We found that Mmp2 was overexpressed by 9.14× in HO, by 39.59× in CM, and by 81.47× in HM. MMP2 is known to degrade extracellular matrix and facilitate cell proliferation and migration. Increased MMP2 expression activity was reported in human PAH [57,58] and in monocrotaline-induced PH in rats [53].
Many of the upregulated inflammatory genes from Table 1 were also reported by other authors as related to PH. For instance, Ccl5 (also known as RANTES; 4.56× in CM and 4.40× in HM), an important chemoattractant for monocytes and T cells, was shown to be increased in PAH [59]. Interestingly, deletion of Ccl5 was shown to attenuate the Sugen-hypoxia model of PH via Cav1-dependent amplification of Bmpr2 signaling [60]. As presented in Figure 2d, we found Bmpr2 as upregulated by 1.78× in HM, but it stayed practically unchanged in HO and CM.
Increased expression of CXCR4, a receptor for the chemokine stromal cell-derived factor 1 (SDF1), was reported in the lungs of patients with IPAH, HPAH, and PAH associated with congenital heart defect [61]. Increased Cxcr4 was also reported in Sugen + hypoxia and monocrotaline + hypoxia rodent models of PH [62]. Inhibition of Cxcr4 moderately attenuated pulmonary vascular remodeling in the Sugen + hypoxia model [63], while overexpression of Cxcr4 participates in the repair of tissue injury [64]. In the present study (Table S1), Cxcr4 was found to be upregulated in CM (4.56×) and HM (2.84×). The lower amplification in HM as compared to CM may be explained by the (even not statistically significant) negative effect of hypoxia (−1.21× in HO). Of note is the upregulation of the tumor necrosis factors and their receptors, confirming the upstream regulator role of the TNF [1].
Ciapin1 (upregulated by 3.16× in CM and by 5.16× in HM, Table S1) promotes anti-apoptosis and cell proliferation via the cyclins D1 (Ccnd1) and E1 (Ccne1), and cyclin-dependent kinases (Cdk) 2 and 4 [65]. Our study confirmed the upregulation for Ccnd1 (1.93× in HM) and Cdk4 (1.49× in HO, 4.26× in CM, and 5.31× in HM). Yet, we found Cdk2 as not regulated and Ccne1 as downregulated (−2.78× in HO). Interestingly, in human VSMC cell culture, CIAPIN1 siRNA was shown to inhibit cell proliferation and enhance apoptosis by increasing Bcl2 and Bax [66]. These results in cell culture suggest that upregulation of Ciapin1 should correlate with downregulation of Bcl2 and Bax. Indeed, we found Bcl2 as downregulated (−2.47× in HM), but Bax was found as upregulated (2.59× in CM and 2.69× in HM). However, the upregulation of Bax is not necessarily a contradiction because Ciapin1 knockdown by siRNA and its upregulation by MCT are triggered by different mechanisms and, hence, not symmetrical phenomena, involving the opposite regulation of the same set of genes. Moreover, as we proved for oligodendrocytes [67] and neurons [68], the heterogeneous cellular environment and external stimuli (missing in the VSMC monoculture but present in the lung) are potent transcriptome modulators.
Among others, we found that Nfkb1 (nuclear factor kappa B subunit 1), one of the controllers of the cytokine production, was upregulated by 3.50× in CM and by 3× in HM (Figure 1d). Icam1 (intercellular adhesion molecule 1) was significantly upregulated by 6.79× in HO, 23.65× in CM, and 43.97× in HM, while Vcam1 (vascular cell adhesion molecule 1) was upregulated by 2.12× in HO, by 6.46× in CM, and by 11.25× in HM. Lgals3 (galectin3), with a critical role in vascular inflammation and fibrosis and heart failure, was also found as upregulated by 10.89× in CM and 17.83× in HM. Importantly, increased expression of Lgals3 and Icam1 was shown to be present in patients with IPAH and connective tissue disease [69].
Anxa 1 (annexin A1, increased by 9.64× in HO, 56.48× in CM, and 106.00× in HM) plays an essential role in cell invasion and migration [70]. It is also an anti-inflammatory protein that controls pro-inflammatory mediator release. It promotes leukocyte detachment from EC and serves as a negative regulator of the transmigratory processes [71,72], and it provides protection from neointima formation in atherosclerosis [73]. Eng (endoglin, 3.03× in CM and 9.00× in HM), a transmembrane receptor for TGFβ signaling, plays a key role in the balance of Alk1 and Alk5 signaling that regulates EC proliferation. Increased expression of Alk1/Eng was reported in EC in IPAH. Furthermore, endoglin deficiency protects mice from hypoxic PH [74]. Nfat5 (nuclear factor activated T cells 5; 3.59× in HM), is implicated in the regulation of genes associated with migration and proliferation [75,76]. Our genomic analysis also revealed increased expression of Edn1 (endothelin 1, 1.88× in HO, 1.83× in CM, and 3.11× in HM), Myc (myelocytomatosis oncogene, 5.92× in CM and 6.42× in HM), Pdgfa (platelet-derived growth factor α, 2.77× in HO, 4.63× in CM, and 6.39× in HM).
The complementary effects of the administration of monocrotaline and exposure to hypoxia are perfectly illustrated by the regulation of the six quantified cyclin-dependent kinases, the key regulatory enzymes of the cell cycle. Thus, Cdkn2a was downregulated in HO and HM groups but not affected in CM, while Cdkn1a was downregulated and Cdkn1b, Cdkn1c, and Cdkn2b were upregulated in CM and HM groups but not in the HO group. Cyclin-dependent kinases are among the targeted genes for the treatment of pulmonary arterial hypertension [77].
As illustrated in panels (d), (e), and (f) in both Figure 3 (VSMC pathway) and Figure 4 (CC pathway), we found significant correlations between the gene expression regulations in the three models. However, the correlation between the regulations in HM and CM models was much stronger than between CM and HO or between HM and HO. We found 23 genes (Figure 2e) whose regulation in HO was opposite to regulation in both CM and HM, and only one gene (Plcb4) whose regulation in HM was opposite to regulation in both CM and HM. Since the HM rats were exposed to both hypoxia and MCT treatment, each upregulating Plcb4, the downregulation of this phospholipase C in HM is surprising and deserves further study.
Mitochondria are the major sources of reactive oxygen species (ROS) and highly susceptible to oxidative stress. Increased ROS production and mitochondrial dysfunction occur in a number of diseases including cardiovascular diseases. A significant increase (3.46×) in the expression of Parp1 (poly (ADP-ribose) polymerase-1) was found in the lungs of rats exposed to both MCT and hypoxia (HM group). Most of the activity of Parp1 is localized in the nucleus. Recent studies showed DNA damage and associated increased expression of Pparp1 to be important aspects of human PAH. Parp1 is implicated in DNA repair, allowing cell proliferation during stress. Parp1 plays an important role in cellular functions during health and disease [78]. Although, in humans, PARP1 was reported to upregulate HIF1A, IL6, and NFAT5 [79,80], we found only Nfat5 as upregulated (see above). Hif1α was downregulated by −1.72× in HO, by −2.04× in CM, and by −1.59× in HM, and Il6 was downregulated by −3.18× in HM.
Recent studies in PAH revealed a metabolic shift toward aerobic glycolysis ("Warburg effect") similar to what was observed in cancer, utilized by proliferating cells for survival. Cav1 was shown to stabilize mitochondria, prevent aerobic glycolysis [35], and negatively regulate NADPH-derived reactive oxygen species [82]. Furthermore, Cav1 modulates Nrf2 expression [83]. It is an important observation that neointimal cells have low Cav1 and normal eNOS expression, resulting in oxidative and nitrosative stress [30]. Compared with the HO group, in the CM and HM groups, there was significantly increased expression of several antioxidants, such as Nfe2l2, Prdx2 (peroxiredoxin 2; 2.66× in HO, 33.94× in CM, and 46.01× in HM), Gpx1, Gpx2 (4.92× in HO, 23.70× in CM, and 24.47× in HM), Txn1, Sod2, Hmox1, and Nqo1. These alterations provide an anti-oxidative and glycolytic milieu for the survival of the proliferating cells.
The analysis of the expression coordination between respiratory and VSMC genes ( Figure 5) revealed several interesting interactions. For instance, in CO, Rhoa has no positive partners but has 15 negative partners and one independent partner (i.e., 0/15/1) among the 32 ATPases, with its coordination power being CP = −74.47%. The situation is reversed in all the three PH models: nine positives, no negatives, no independents (9/0/0) with a CP of +61.07% in HO, 23/0/1, a CP of +66.37 in CM, 6/0/3, and a CP of +63.93 in HM. Importantly, 10 of the positively correlated ATPase partners of Rhoa in CM were negative partners in CO. This switch from negative to positive correlation of Rhoa with the ATPases justifies the reported "role of the Nox4-derived ROS-mediated RhoA/Rho kinase pathway in rat hypertension induced by chronic intermittent hypoxia" [84]. Confirming other reports of increased Rhoa in PH [83], we found Rhoa as upregulated by 7.00× in HO, 6.78× in CM, and 9.21× in HM. Our results complete the results of a previous report [85] by showing that Rhoa is upregulated and plays an important role regardless of what caused the PH. Moreover, we found also that Rhoa's action is mediated by the ATPases, a result that needs further investigation, especially in the light of the recent discovery of pathogenic mutations of the (not quantified in our study) ATP13A3 gene in IPAH [86].
The gene pair prominence (GPP) analysis demonstrated the plasticity of the genomic fabrics and their interplay. Figure 6 illustrates how the GPP landscape of the interaction between the vascular smooth muscle contraction (VSMC) genes and the ATPases is remodeled in the three PH conditions with respect to control. Interestingly, while, in control, a few VSMC genes interact with almost all ATPases, in all three PH models, almost all VSMC genes interact with a few ATPases, making the vascular smooth muscle contraction more vulnerable to certain point alteration within the respiratory chain. This new type of analysis can be used to identify the most important gene pairs in specific interactions.
The gene pair prominence (GPP) analysis demonstrated the plasticity of the genomic fabrics and their interplay. Figure 6 illustrates how the GPP landscape of the interaction between the vascular smooth muscle contraction (VSMC) genes and the ATPases is remodeled in the three PH conditions with respect to control. Interestingly, while, in control, a few VSMC genes interact with almost all ATPases, in all three PH models, almost all VSMC genes interact with a few ATPases, making the vascular smooth muscle contraction more vulnerable to certain point alteration within the respiratory chain. This new ty

Conclusions
All PH groups exhibited abnormal expression of genes involved in inflammation, cell proliferation, and vascular smooth muscle contraction. Importantly, monocrotaline-treated rats that were exposed to hypoxia, which shows similar changes to PAH, displayed evidence of mitochondrial damage and activation of antioxidants, which are likely to cause metabolic switch and provide a survival milieu for the proliferating cells. Our finding of the reversed correlation of Rhoa with ATPases from negative in control rats to positive in PH suggests a potential therapeutic avenue by targeting Rhoa and/or its partners, as suggested in 2008 [87].

Conflicts of Interest:
The authors declare no conflict of interest