Two common databases, ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), were used to search datasets of IPAH patients. We used the following search terms: pulmonary hypertension or pulmonary arterial hypertension or PAH or PH or idiopathic pulmonary arterial hypertension or IPAH. Two researchers independently screened the two databases and selected relevant datasets based on the same following criteria:

Our workflow for bioinformatics analysis of publicly available datasets is illustrated in Figure 1.

The gene expression profiling was annotated through the corresponding annotation package and the R software (version 3.6.3). The expression data were further normalized by using the “limma” package (Supplementary Figure 1). Then, we performed RRA analysis to rank the upregulated and downregulated expression genes through the R package of “Robust Rank Aggregation” and listed the genes according to their fold changes (11). In the final list, a Bonferroni correction was performed to cut down the false-positive results, and p-values that represented the possibility of ranking high were calculated for each gene.

We performed GO and KEGG pathway enrichment analyses by using the cluster profile function of R software to match the biological themes of the gene clusters, with a cutoff p-value of 0.05. The STRING database (https://string-db.org/cgi/input.pl) was used to establish a PPI network. In the network, the significant Molecular Complex Detection (MCODE) and hub genes were present by using Cytoscape software (version 3.8.0).

To perform weighted gene coexpression network analysis (WGCNA), genes with both p < 0.05 and logarithmic fold changes (logFCs) > 0.25 were selected from RRA results. To improve the number of samples and avoid generating less reliable results, we integrated and normalized the four datasets by batch normalization using the “sva” and “limma” package in the R computing environment. We transformed the adjacency matrix into a topological overlap matrix (TOM) and split the genes into different modules based on the TOM-based dissimilarity measure. Key modules were identified by setting the soft-thresholding power of 6 (scale-free R² = 0.9), cut height of 0.25, and minimal module size of 20. We further set gene significance (GS) > 0.3 and module membership (MM) > 0.8 to define the hub genes in WGCNA.

We downloaded the microRNA expression dataset GSE67597 and identified pulmonary differentially expressed miRNAs (DE-miRNAs) in IPAH patients by using the GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) with the thresholds of both p < 0.05 and |log2FC| > 2. The miRNet database, a useful tool for miRNAs analysis, was used to predict the target genes of DE-miRNAs (12). To visualize the relationship between DE-miRNAs and predicted target mRNAs, we established a miRNA–mRNA network by using Cytoscape (version 3.8.0).

Twelve male C57BL/6 mice (8–10 weeks old) were randomly assigned to the following two groups (n = 6 per group): CH+SU group: mice were exposed to chronic hypoxia (10% O₂) for 28 days, weekly receiving a subcutaneous injection of 20 mg/kg of SU5416 (an angiogenesis inhibitor), which was dissolved in a carboxymethylcellulose solution; the normoxic control group (Nor): mice were cultured under normoxic condition and weekly received a subcutaneous injection with an equivalent volume of the dissolvent solution according to their weights, as shown in Supplementary Figure 2. Two groups were provided with food and water ad libitum. At the end of the experiment, we used isoflurane to anesthetize the mice and measured their right ventricular systolic pressure (RVSP) by using a Millar pressure transducer catheter through right-sided heart catheterization, as we previously described (13). To evaluate the RV remodeling, we weighed the wall of the RV and the left ventricle plus septum (LV+S) to calculate the ratio of (RV/LV+S).

The whole heart and left lung were embedded in paraffin, sliced into 5-μm sections, and then stained with hematoxylin and eosin (H&E) and Masson; and immunofluorescence followed examination with a light microscope (Nikon).

Total RNA from lungs was isolated using TRIzol, and the extracted RNA was reverse transcribed into complementary DNA (cDNA) using the PrimeScript RT reagent kit (Takara Bio). RT-PCR was carried out using SYBR Green Premix Ex Taq (Takara Bio) to detect the expressions of the genes, which were selected from our analysis. The primer sequences were provided in supplements.

We used R software and GraphPad Prism 7 to perform statistical analyses, presented data as mean ± standard error, and calculated statistical significance by the Student t-test. p < 0.05 were considered statistically significant.

After filtering ArrayExpress and GEO based on the eligibility criteria, five microarray datasets of PH were finally selected. The information of these datasets was listed, including study country, types of RNA source, GEO accession ID, detection platforms, experiment type, and sample information. The number of IPAH patients in each study ranged from 6 to 32, and the controls ranged from 8 to 25. Finally, a total of 64 IPAH patients and 58 controls were involved in our study (Table 1).

A total of 183 upregulated and 77 downregulated differentially expressed genes (DEGs) were successfully identified by integrating four genome-wide gene expression datasets through the RRA method (Supplementary Table 1). Among the upregulated genes, PKP2 was ranked as the first one (p = 4.71E−07, adjusted p = 1.36E−02), followed by ALAS2 (p = 1.11E−06, adjusted p = 3.19E−02). Meanwhile, RNASE2 (p = 8.94E−10, adjusted p = 2.57E−05) and PROK2 (p = 2.20E−08, adjusted p = 6.33E−04) were ranked as the first and second downregulated genes in RRA analysis. The top 30 upregulated and downregulated genes in IPAH were illustrated by a heatmap (Figure 2). Among these genes, the roles of some ones had been well-explored in PAH, such as interleukin 13 receptor alpha 2 (IL13RA2), angiotensin I converting enzyme 2 (ACE2) (14), and vascular cell adhesion molecule 1 (VCAM1) (15). Notably, some genes were novel, and their functions in PAH had not been researched in published literature, such as hemoglobin alpha 2 (HBA2), frizzled-related protein 2 (SFRP2), and ribonucleases 2 (RNASE2).

To research the potential biological pathways involved in IPAH, we performed GO and KEGG to analyze the DEGs. Several enrichment biological processes in GO terms were identified, such as leukocyte migration, cell chemotaxis, and myeloid leukocyte chemotaxis (Figure 3A). In terms of molecular function, receptor-ligand activity was listed as the most significantly GO terms (Figure 3B). Moreover, some cellular component terms were enriched, such as the external side of plasma membrane and collagen-containing extracellular matrix (ECM) (Figure 3C). According to KEGG analysis, we found that the DEGs were mostly associated with cytokine–cytokine receptor interaction, chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, fluid shear stress and atherosclerosis, and IL-17 signaling pathway (Figure 3D).

Next, PPI network of these DEGs was constructed by using the STRING database. The visualization was carried out using Cytoscape software (Figure 4A). MCODE plugin was used to find the top hub genes; the top 3 closely connected modules were identified. The genes in MCODE 1 include CCL21, SAA1, CCL4L1, ADORA3, NMUR1, PPBP, CXCR2, CXCL10, CXCR1, CXCL9, CCR1, AGTRL, PF4, CX3CL1, C5, and CX3CR1, which are associated with Class A/1 (Rhodopsin-like receptors), peptide ligand-binding receptors, and G alpha (i) signaling events; the genes in MCODE 2 are related with positive regulation of leukocyte activation, allograft rejection, and cellular response to interferon-gamma; and the genes in MCODE 3 are associated with oxygen transport and gas transport (Figures 4B–D, Supplementary Table 2).

To further identify the functional and meaningful modules most associated with IPAH, we utilized WGCNA in these four datasets (Figure 5A). We selected genes with p < 0.05 and logFCs > 0.25 from the ranked gene list to cover enough genes in WGCNA. Finally, 14 modules were identified by setting the soft-thresholding power as 6 (scale-free R² = 0.9) and cut height as 0.25 (Figures 5B–D; non-clustering DEGs shown in gray). The correlations between module and IPAH were illustrated in a heatmap, and the midnight-blue module was found to be the most highly correlated module with IPAH (Figures 5E,F). The midnight-blue module contained 42 genes shown in Figure 5G (correlation coefficient = 0.49, p = 1E−08G). We selected top 20 hub genes from the midnight-blue module: CXCL10, IFNG, CCL3, CXCL9, VCAM1, ANKRD22, IL13, SAA1, ITGAMT, TNFSF8, NCR1, EOMES, KLRC3, SLC14A1, ADORA3, ZAP70, CD96, CCL5, LCK, and BTLA. Venn diagram displayed the overlap of genes between PPI analysis and WGCNA. The results showed that seven hub genes (INFG, IL13, CXCL10, SAA1, CCL3, CXCL9, and ADORA3) were found in the top 20 hub genes of PPI and top 20 genes of the midnight-blue module (Figure 5H).

It is widely acknowledged that miRNA exerts an essential effect on PAH by negatively regulating target mRNA. Therefore, we utilized the GEO2R tool to identify the DE-miRNAs in IPAH and predicted the downstream target genes of DE-miRNAs by using the miRNet database. Eventually, we predicted 1,816 and 773 target genes for the 13 upregulated DE-miRNAs and nine downregulated DE-miRNAs, respectively. Furthermore, the upregulated DE-miRNA–target gene network and downregulated DE-miRNA–target gene network were established and presented (Figures 6A,D). The number of target genes for each DE-miRNA is listed in Figures 6B,E. The common genes in DEGs and predicted target genes were identified and shown by the Venn diagram (Figures 6C,F). As a consequence, we found that four target genes (LILRA2, RPS4Y1, SFN, and VENTX) were negatively regulated by increased miRNA-500a-3p, miRNA-31-5p, or miRNA-6074; and eight upregulated genes (SLFN5, ANKRD50, XIAP, SYTL3, PDE4D, VCAM1, FGD4, and SECISBP2L) were the target of miRNA-1178-3p, miRNA-302f, or miRNA-495-5p (Figure 7).

Chronic hypoxia exposure with SU5416 injection is a well-studied and published method to make animal models of PAH. However, there are currently no perfect animal models that replicate entirely human PAH. To investigate whether the pathological changes of IPAH in human also occur in CH+SU mice and to validate the expressions of hub genes found in our analysis, we established PAH model. After a 28-day exposure to CH with the administration of SU5416, mice had a significant elevation in RVSP (Figure 8A) and RV/LV+S compared with those in control mice (Figures 8B,C). Meanwhile, CH+SU treatment resulted in pulmonary vascular remodeling, evidenced by the increased media (Figures 8D–F). Then, we detected the expression of DEGs that interest us in a successfully established PAH mice model. As shown in Figure 9, the expressions of PKP2 (Figure 9A) were significantly upregulated. At the same time, ADORA3, PROK2, and IL-13 (Figures 9B–D) were downregulated in the lungs from CH+SU mice, which were consistent with the results in the datasets of IPAH patients. Meanwhile, we found inconsistent gene expression between the animal model and IPAH patients; for example, CXCL10 (Figure 9E) was downregulated, while SFN (Figure 9F) was upregulated in CH+SU mice, which were contrary to the findings in IPAH patients. The expression of SFRP2 showed no difference between normoxic and CH+SU mice (Figure 9G) but significantly upregulated in IPAH patients based on RRA analysis.