Rat PASMCs isolation and cell culture
Rat PASMCs were isolated as previously described.12 Briefly, six-week old Sprague-Dawley rat purchased from Charles River (Beijing, China) was sacrificed, which was approved by the Ethics and Animal Care and Use Committee of Henan University. The lungs were excised and rinsed with phosphate-buffered saline (PBS). The pulmonary arteries were isolated followed by removal of the adventitia under the microscope. Minced arteries were attached to bottom of petri dish and then immersed by Dulbecco’s modified eagle medium/F12 containing 20% fetal bovine serum, 100 U/mL penicillin and 100 μg/mL streptomycin in a 37 °C, 5% CO2 humidified incubator. Three days later, non-adherent cells were removed, and the adherent cells that had grown to 90% confluence were considered as passage 0 PASMCs. Passages 3 were used for the subsequent hypoxia experiments.
Rat PASMCs at a density of 1x105/mL were seeded into 6-well plates and then placed in a 37 °C, 5% CO2 humidified incubator. 24 h later cells were starved with Dulbecco’s modified eagle medium/F12 containing 0.5% fetal bovine serum. Medium were refreshed and cells were then either cultured in normoxia condition (21% O2) or hypoxia condition (1%O2) in the incubator 24 h post-starvation. 24 h later cells under normaxia or hypoxia (n = 6 per group) were collected respectively. Cells were washed and centrifuged at 300 g at 4°C for 10 minutes and resuspend with 1 mL PBS. Cells were then counted with erythrocytometers and diluted to a final density of 1x105/mL with PBS. 1 mL cell suspension per sample was transferred to 1.5 mL clean tube and centrifugated at 300 g at 4°C for 5 minutes. Supernatants were discarded and cell pellets were stored at -80°C before further use. No repeated freezing and thawing occurred before sample processing to avoid potential degradation risks of metabolites.
1mL of MeOH:ACN:H2O (2:2:1, v/v) solvent mixture was added into the cell samples. The cell mixtures were vortexed for 30 sec followed by sonification at 4℃ in the water bath for 10 min. The mixtures were transferred to liquid nitrogen for 1 min and thawed at room temperature followed by sonification at 4℃ in the water bath for 10 min. After the repeat for 3 times, mixtures were put at -20°C for 1 h to facilitate protein precipitation. Mixtures were centrifuged at 13,000rpm at 4°C for 15 min, supernatants were discarded, and pellets were vacuum dried and reconstituted with 100 μL of ACN:H2O (1:1, v/v). Mixtures were vortexed for 30 sec and sonicated as previously stated followed by centrifugation at 13,000rpm at 4°C for 15 min. Supernatants were stored at -80 °C prior to Liquid chromatography-mass spectrometry (LC/MS) analysis.
Liquid chromatographic separation for processed plasma was achieved on a ZORBAX Eclipse Plus C18 column (2.1 × 100 mm, 3.5 μm, Agilent, USA) maintained at 45 °C, whereas mass spectrometry was performed on a Nexera X2 system (Shimadzu, Japan) coupled with a Triple TOF 5600 quadrupole-time-of-flight mass spectrometer (AB SCIEX, USA). The temperature of the sample chamber was maintained at 7°C. The gradient elution steps were shown in Supplementary Table 1. The injected sample volume was 10 μL for each run in the full loop injection mode, and the flow rate of the mobile phase was 0.5 mL/min. The mobile phase A was mainly composed of water and contains 0.1% formic acid. The mobile phase B was mainly composed of acetonitrile and contains 0.1% formic acid. Water (LC-MS grade) and acetonitrile (LC-MS grade) were purchased from Fisher Scientific. The purity of formic acid purchased from Acros was greater than 98%.
Data availability and process
In terms of metabolomics study, data preprocessing was performed before pattern recognition. The original data was processed by the instrument's own metabolomics processing software Progenesis QI (Waters Corporation, Milford, USA) for baseline filtering, peak identification, integration, retention time correction, peak alignment and normalization. Finally, a data matrix of retention time, mass-to-charge ratio and peak intensity were obtained. The integrated data matrix was imported into the SIMCA-P+ (v 14.1) software package (Umetrics, Umeå, Sweden), and partial least squares discriminant analysis (PLS-DA) was used to distinguish the overall difference in metabolic profile between groups. In PLS-DA analysis, variables with a variable weight value (Variable Important in Projection, VIP) > 1 were considered to be distinguishing among groups.
With regard to transcriptomic analysis, we accessed the data set GSE190913 based on GPL1674 including mouse lung tissues during development of hypoxia-induced PH (days 1-21) and resolution of PH after return to normoxia (days 22-35). Mice were sampled during nine time-points and there were 4 replicates at each time-point. We selected samples at day 1, day 14, day 21 and day 35 after hypoxia for further analysis. Differential expression analysis was performed using the limma package in R (v3.6.3.).14 P value < 0.05 and |fold change (FC)| > 1.5 were considered statistically significant for identifying differentially expressed genes (DEGs). To validate the hub gene in lung tissues from PH patients, we selected two data sets GSE11726115 consisting of 25 control subjects and 58 PH patients, and GSE2498816 including 22 controls and 62 PH patients. Both were based on GPL6244. The ComBat function in the sva package was used to merge these two data sets and remove the batch effect.17 In addition, single cell RNA sequencing data from lung tissues of monocrotaline, Sugen 5416/hypxoixa (SuHx) or control rats were downloaded from an open-access online platform (http://mergeomics.research.idre.ucla.edu/PVDSingleCell/).18
Connectivity Map analysis
The Connectivity Map (CMap) (https://portals.broadinstitute.org/cmap) is an open resource that connect genes, small molecules, and disease by virtue of common gene-expression signatures.19 Drugs with Mean < −0.4 and p < 0.05 in CMap analysis were considered to be potential small molecular compounds that can reverse altered expression of user-defined gene list in cell lines.
Identification of hub genes and network interaction visualization
Top 10 enriched metabolite pathways were identified by enrichment analysis of metabolite sets with online open-access platform Metaboanalyst20 (v5.0, https://www.metaboanalyst.ca). Metabolism associated genes (MAGs) were acquired from Genecards21 (https://www.genecards.org) database. The terms of top 10 enriched metabolite pathways were used as key word for analysis, and 1259 genes with relevance score > 8 were regarded as MAGs. The intersection of MAGs and DEGs in mouse lung tissues of hypoxia induced PH was listed as gene set 1. The intersection of MAGs and DEGs in smooth muscle cells of SuHx versus control rats at single cell level was listed as gene set 2. The union set of these two gene sets were regarded as hypoxia induced metabolism associated genes. Protein protein interaction were analyzed and visualized by online tool STRING (v11.0, https://string-db.org).22 Hub genes were identified by cytoHubba plugin with MCC algorithm or by MCODE plugin in Cytoscape (v 3.8.2).23
Construction of LASSO model and receiver operating characteristic (ROC) curve analysis
To distinguish PH patients from control subjects, Least Absolute Shrinkage and Selection Operator (LASSO) model was constructed with the expression profiles of 6 metabolism associated hub genes by glmnet package in R. A model index for each sample was generated using the regression coefficients from the LASSO analysis to weight the expression value of the selected genes. GSE117261 data set were randomly assigned to training set (70%) and test set (30%). ROC curves were generated to evaluate the ability of LASSO model to identify PH by ROCR package.
Data visualization and statistics
The sample distribution of hypoxia exposed rat PASMCs (Hx) and control PASMCs (Nor), or mouse lung tissues under hypoxia at the indicated days were visualized in scatter plot with R package ggplot2. The distribution of data sets GSE117261 and GSE24988 before and after removal of batch effect was also visualized in scatter plot with ggplot2 package. Distinct metabolites/genes between groups were visualized in volcano plot. Heatmap was generated to visualize the expression of all the distinct metabolites or indicated genes in each individual sample using pheatmap package in R. Venn diagram was generated to visualize the overlap of the indicated gene sets using R package VennDiagram. Pathway enrichment of indicated genes were analyzed in functional annotation bioinformatics tool DAVID24 (v6.8, https://david.ncifcrf.gov/tools.jsp) and were then plotted with barplot function in R. Uniform manifold approximation and projection (UMAP) was generated based on Seurat_umap.coords data to visualize each clusters under different conditions. Hub gene netwok was visulized by Cytoscape.
Data are presented as the mean ± standard error of them mean (SEM). When two groups were compared, statistical differences were assessed with unpaired 2-tailed unpaired t-test if normally distributed. Otherwise, Mann-Whitney U test was utilized. Comparisons of more than three groups were performed by analysis of variance (ANOVA) and Tukey’s post-hoc test or Kruskal-Wallis test, as appropriate (GraphPad Prism 8). A two-sided P value of ≤ 0.05 was used to determine significant differences among groups.