Introduction The Clinical Challenge of High-Altitude Pulmonary Hypertension

Ascent to elevations exceeding 3,000 meters exposes the pulmonary vasculature to sustained hypoxic stress.1,2 The physiological response—hypoxic pulmonary vasoconstriction (HPV)—serves an adaptive purpose: redirecting blood flow from poorly ventilated lung regions to optimize gas exchange. This protective reflex becomes maladaptive under chronic hypoxia. Sustained vasoconstriction triggers vascular remodeling, medial hypertrophy, and ultimately fixed elevation of pulmonary artery pressure.3,4 The clinical syndrome that emerges—hypoxia-induced pulmonary hypertension (HPH)—carries substantial morbidity and mortality when left untreated.

The pathophysiological cascade has been extensively characterized. Hypoxia-inducible factors (HIFs), particularly HIF-1α and HIF-2α, orchestrate transcriptional responses to low oxygen.5,6 These master regulators drive expression of genes promoting glycolysis, erythropoiesis, and angiogenesis while suppressing oxidative metabolism. The signaling network extends to involve endothelin-1, serotonin, bone morphogenetic protein receptor type 2 (BMPR2), and numerous other mediators.7,8 Despite this mechanistic understanding, translation to effective therapies has proven challenging. Phosphodiesterase-5 inhibitors, endothelin receptor antagonists, and prostacyclin analogs address vascular tone but fail to reverse established remodeling.9,10 Epidemiological data underscore the problem’s scope. Conservative estimates place 10–15% of individuals ascending above 4,500 meters at risk for clinically significant HPH.11 The global population living permanently above 2,500 meters approximates 140 million—concentrated on the Tibetan Plateau, Andean highlands, and Ethiopian uplands.12 Military operations, mining activities, and tourism add transient exposures. Chronic mountain sickness, of which pulmonary hypertension constitutes a major component, affects up to 18% of Han Chinese migrants to Tibet.13 Economic consequences—healthcare expenditure, reduced productivity, disability—amplify the public health imperative for effective interventions.

Metabolic Reprogramming in Pulmonary Vascular Disease

A dimension of HPH pathophysiology receiving increased attention involves cellular metabolism.14,15 Hypoxic cells face a fundamental bioenergetic challenge: generating ATP when the terminal electron acceptor for oxidative phosphorylation becomes limiting. The canonical response involves metabolic reprogramming toward glycolysis—what has been termed “Warburg-like” rewiring—a metabolic shift toward glycolysis, analogous to the phenotype observed in cancer cells—to address the bioenergetic challenge of hypoxic conditions.16,17

This metabolic switch carries several implications. Glycolysis generates ATP rapidly but inefficiently (2 ATP per glucose versus 36 from complete oxidation). The accompanying lactate production acidifies the microenvironment. Perhaps most significantly, glycolytic intermediates feed biosynthetic pathways supporting cell proliferation—the pentose phosphate pathway for nucleotide synthesis, serine-glycine metabolism for one-carbon units.18,19 Proliferating pulmonary artery smooth muscle cells (PASMCs) and adventitial fibroblasts drive the medial hypertrophy and neointima formation characteristic of vascular remodeling.20 Mitochondrial dysfunction extends beyond simple substrate limitation. Hypoxia paradoxically increases mitochondrial reactive oxygen species (ROS) production, particularly from complex III.21,22 These ROS serve signaling functions—stabilizing HIF-1α, activating calcium channels—but also impose oxidative damage. The antioxidant systems tasked with ROS neutralization, particularly glutathione and taurine-dependent pathways, face depletion under chronic oxidative stress.23,24

Amino acid metabolism intersects critically with vascular function. Arginine serves as substrate for nitric oxide synthase (NOS); its depletion compromises endothelium-dependent vasodilation.25,26 Arginase upregulation—documented in HPH—diverts arginine toward ornithine and polyamine synthesis, further reducing NO bioavailability.27 Branched-chain amino acid catabolism fuels the TCA cycle during metabolic stress but may impair protein synthesis capacity.28

Importantly, while metabolic reprogramming represents a central feature of hypoxia-induced pulmonary hypertension, the disease is multifactorial in nature. In addition to metabolic alterations, hypoxic pulmonary vasoconstriction, endothelial dysfunction, aberrant vascular remodeling, inflammatory activation, oxidative stress, and canonical hypoxia signaling pathways such as HIF-mediated transcriptional programs all contribute critically to disease initiation and progression.4–11 These interconnected mechanisms collectively shape pulmonary vascular remodeling under chronic hypoxic exposure. Accordingly, metabolic dysregulation should be interpreted as one key but non-exclusive component within a broader pathological network rather than an isolated driver of disease.

Traditional Chinese Medicine and Systems Pharmacology

Traditional Chinese medicine (TCM) approaches disease from a systems perspective fundamentally distinct from Western reductionist pharmacology.29,30 Rather than targeting single molecular entities, TCM formulas employ multi-component mixtures hypothesized to engage multiple pathways synergistically. This therapeutic philosophy aligns conceptually with network pharmacology—the emerging recognition that complex diseases may require multi-target interventions.31,32

Shengxian Decoction (SXT) exemplifies this paradigm. Documented in Zhang Xichun’s “Yi Xue Zhong Zhong Can Xi Lu” (Medical Records Combining Chinese and Western Medicine), the formula addresses “qi deficiency leading to sinking” (qi xu xia xian)—a syndrome characterized by fatigue, dyspnea, and cardiovascular collapse.33 The classical presentation maps remarkably well to hypoxia-related cardiopulmonary dysfunction, suggesting potential therapeutic relevance despite the pre-modern conceptual framework.

The formula comprises five herbs with distinct pharmacological profiles. Huangqi (Astragalus membranaceus), the principal component at 18 g, yields astragaloside IV—a cycloartane-type saponin with documented anti-inflammatory, antioxidant, and cardioprotective properties.34,35 Mechanistic studies have implicated PI3K/Akt signaling, Nrf2-mediated antioxidant response, and AMPK activation in its effects.36 Zhimu (Anemarrhena asphodeloides) contributes steroidalsaponins (timosaponins) that moderate Astragalus’s warming properties while exerting anti-inflammatory effects via NF-κB inhibition.37 Jiegeng (Platycodon grandiflorus) provides platycodin D and related triterpenoid saponins traditionally believed to “guide” other constituents to the lung; modern studies confirm pulmonary tropism and mucolytic activity.38 Chaihu (Bupleurum chinense) yields saikosaponins with hepatoprotective and immunomodulatory properties.39 Shengma (Cimicifuga foetida) contributes triterpene glycosides implicated in anti-inflammatory and vasodilatory effects.40 Importantly, accumulating experimental evidence suggests that SXT exerts protective effects in hypoxia-related cardiovascular disorders. Previous studies have shown that SXT can improve myocardial energy metabolism through activation of the AMPK/PGC-1α pathway, increase ATP production, enhance mitochondrial respiratory chain activity, and alleviate chronic heart failure–associated cardiac dysfunction.41,42 In addition, SXT has been reported to protect the microvascular endothelium by upregulating VE-cadherin and tight-junction proteins, promoting angiogenesis, and preserving vascular basement membrane integrity.41,43 These findings suggest that SXT may influence several biological processes directly relevant to pulmonary hypertension, including energy metabolism, vascular remodeling, endothelial function, and hypoxia adaptation. However, whether these effects occur in HPH and the molecular networks underlying such actions remain largely unknown. While network pharmacology predicts SXT engages pathways such as HIF-1 and PI3K-Akt signaling, unbiased multi-omics profiling is needed to validate these predictions and identify unrecognized convergent nodes.44

Multi-Omics Integration for Mechanism Elucidation

Multi-omics strategies offer powerful tools for dissecting complex therapeutic mechanisms.45,46 Metabolomics provides a functional readout of cellular state—the downstream product of gene expression, protein activity, and environmental interaction.47,48 Unlike genomics or transcriptomics, which capture regulatory potential, metabolomics reflects what is actually happening biochemically. This functional proximity to phenotype makes metabolomics particularly informative for pharmacological studies.49

Transcriptomics complements metabolomics by revealing the regulatory logic underlying observed changes.50 Gene expression profiling identifies which pathways are being transcriptionally activated or repressed. Integration of the two layers enables correlation analysis linking metabolite alterations to their putative enzymatic and regulatory determinants.51

Despite increasing evidence supporting the pharmacological activities of SXT and its constituent compounds, several critical knowledge gaps remain. First, the molecular basis of SXT’s protective effects in hypoxia-induced pulmonary hypertension has not been systematically characterized. Second, it remains unclear how hypoxia-induced metabolic remodeling is coordinated with transcriptional alterations and whether SXT can reverse these changes in an integrated manner. Third, potential metabolite-gene signatures that may serve as biomarkers for monitoring therapeutic responses have not been explored. These gaps provide a strong rationale for applying an integrated metabolomic-transcriptomic strategy to investigate the systems-level effects of SXT in HPH. In particular, we hypothesize that pyruvate-centered metabolic regulation, involving key enzymes such as pyruvate kinase M2 (PKM2) and lactate dehydrogenase A (LDHA), may represent a central convergence point linking glycolytic reprogramming and mitochondrial dysfunction in hypoxia-induced pulmonary vascular remodeling, and may represent one potential pathway modulated within this broader multi-mechanistic disease context affected by SXT.

Study Objectives

The present study employs integrated metabolomic-transcriptomic analysis to interrogate SXT’s therapeutic mechanism in a rat model of HPH. Three specific aims guide our investigation:

(i) To delineate the metabolic and transcriptional signatures of hypobaric hypoxia–induced pulmonary injury, thereby establishing a pathophysiological baseline for evaluating therapeutic effects.

(ii) Delineate the dose-dependent effects of SXT on metabolomic and transcriptomic profiles, testing whether molecular restoration scales with dosing and which pathways respond most robustly.

(iii) Construct a cross-platform network model integrating metabolite alterations with gene expression changes, identifying hub nodes and pathway convergence points that may represent key therapeutic targets.

Through this approach, we aim not only to characterize the molecular alterations associated with HPH and their modulation by SXT, but also to identify potential metabolite-gene networks and candidate biomarkers associated with treatment response. By integrating transcriptomic and metabolomic data, this study seeks to provide systems-level evidence for the biological effects of SXT and generate testable hypotheses regarding its potential mechanisms of action in HPH.

Materials and Methods Reagents and Materials

Methanol, acetonitrile (LC-MS grade), and formic acid were purchased from Fisher Scientific (Pittsburgh, PA, USA). Leucine-enkephalin was obtained from Waters Corporation (Milford, MA, USA). Authentic metabolite standards for taurine, pyruvate, lactate, glutathione, and branched-chain amino acids were purchased from Sigma-Aldrich (St. Louis, MO, USA). TRIzol reagent was from Invitrogen (Carlsbad, CA, USA). Sildenafil was purchased from Pfizer Inc. (New York, NY, USA). Urethane was purchased from Sigma-Aldrich (St. Louis, MO, USA). Hematoxylin and eosin (H&E) staining kit was purchased from Solarbio Life Sciences (Beijing, China). BL-420N biological signal acquisition system was purchased from Chengdu Taimeng Technology Co., Ltd. (Chengdu, China). All other chemicals were of analytical grade. Crude herbs for Shengxian Decoction were procured from Tongrentang Pharmaceutical (Beijing, China).

Animals and Experimental Design

Thirty male Sprague-Dawley rats (200–250 g, 8–10 weeks old) were obtained from the Experimental Animal Center of Qinghai University Clinical Medical College. Animals were housed under controlled conditions (22±2°C, 50±10% relative humidity, 12-h light/dark cycle) with ad libitum access to standard laboratory chow and water. All animal experiments were conducted between January 2024 and March 2025. All experimental procedures were approved by the Ethics Committee of Qinghai University Affiliated Hospital (School of Clinical Medicine) (Approval No. P-SL-2023-453) and were conducted in accordance with the Guide for the Care and Use of Laboratory Animals (NIH Publication No. 85–23, revised 2011). All procedures complied with the guidelines for the care and use of laboratory animals issued by the American Veterinary Medical Association (AVMA).

Following one-week acclimatization, rats were randomly assigned to six groups (n=6 per group) using a computer-generated randomization sequence: (1) Control—normoxic housing with vehicle gavage; (2) Model—hypobaric hypoxia with vehicle gavage; (3) XDNF group—hypobaric hypoxia with sildenafil treatment (30 mg/kg/day); (4) SXT-L—hypobaric hypoxia with low-dose SXT (1.8 g/kg/day); (5) SXT-M—hypobaric hypoxia with medium-dose SXT (3.6 g/kg/day); (6) SXT-H—hypobaric hypoxia with high-dose SXT (7.2 g/kg/day).

After randomization, each animal was assigned a unique identification code for tracking throughout the study. Allocation concealment and blinded procedures were strictly implemented during the entire experimental process. Drug administration was performed independently by designated personnel who were not involved in subsequent analyses. Personnel responsible for sample preparation, histopathological evaluation, metabolomic profiling, transcriptomic sequencing, and statistical analyses were blinded to group allocation throughout the study. Group information was disclosed only after completion of all experimental procedures, data acquisition, and statistical analyses.

Sildenafil (XDNF) was included as a positive control because it is a clinically established phosphodiesterase-5 (PDE5) inhibitor widely used for pulmonary hypertension treatment. Inclusion of the XDNF group enabled validation of model responsiveness and comparative evaluation of the therapeutic efficacy of SXT.

SXT was originally described by Zhang Xichun in the classical medical text Yixue Zhongzhong Canxilu. The original prescription consists of 6 qian of Astragali Radix, 3 qian of Anemarrhenae Rhizoma, 1.5 qian of Bupleuri Radix, 1.5 qian of Platycodonis Radix, and 1 qian of Cimicifugae Rhizoma. According to the traditional Qing Dynasty weight system, 1 qian is equivalent to approximately 3.125 g. Based on this conversion, the total daily clinical dose for an adult is approximately 40.64 g. The rat dose was determined according to the body surface area (BSA)-based interspecies dose conversion method, in which the equivalent dose for rats is 6.25 times the human dose on a body-weight basis. Assuming an average adult body weight of 70 kg, the clinical human dose was calculated as 40.64 g/70 kg/day. The corresponding rat-equivalent dose was calculated as follows: Rat dose (g/kg/day) = 6.25 × Human dose (g/kg/day) = 6.25×40.64 / 70 = 3.6 g/kg/day. Accordingly, 3.6 g/kg/day was designated as the medium-dose SXT group. The low-dose and high-dose groups were set at 0.5-fold (1.8 g/kg/day) and 2-fold (7.2 g/kg/day) of the equivalent dose, respectively.

Doses were calculated based on BSA conversion from clinical human doses, representing 0.5×, 1×, and 2× the human-equivalent dose, respectively. Sample size was determined by power analysis based on preliminary data. Assuming a 30% difference in key metabolites between groups, with a significance level of α = 0.05 and a statistical power of 0.80, a minimum of five animals per group was required. To account for potential attrition, six rats were included in each group. Although the sample size of n = 6 per group is relatively small, it is consistent with the commonly accepted scale used in small-animal omics studies of pulmonary hypertension. Group sizes of 4–6 animals are generally considered statistically appropriate for paired transcriptomic and metabolomic analyses and have been widely adopted in previous studies.52 No animals were excluded during the study; therefore, all six rats in each group were included in the final analyses. Importantly, all transcriptomic, metabolomic, histological, and biochemical assessments were performed using tissues collected from the same cohort of rats, ensuring direct correspondence among datasets generated from different analytical platforms. In addition, all experimental procedures, including disease modeling, drug administration, and tissue collection, were conducted according to a standardized protocol throughout the study, thereby minimizing inter-individual variability and enhancing data comparability.

Preparation of Shengxian Decoction

SXT was prepared according to the classical formula proportions: Astragalus membranaceus (18 g; batch No. 200501; Huayun Chinese Medicinal Decoction Pieces Co., Ltd., Bozhou, China), Anemarrhena asphodeloides (9 g; batch No. 201001; Huayun Chinese Medicinal Decoction Pieces Co., Ltd)., Platycodon grandiflorus (4.5 g; batch No. 190501; Huayun Chinese Medicinal Decoction Pieces Co., Ltd)., Bupleurum chinense (4.5 g; batch No. 200501; Huayun Chinese Medicinal Decoction Pieces Co., Ltd)., and Cimicifuga foetida (3 g; batch No. 160101; Huayun Chinese Medicinal Decoction Pieces Co., Ltd). These quantities correspond to one dose for a 70 kg adult human.

Authenticated herbs were combined and decocted twice in 10 volumes (w/v) of distilled water for 1 hour per decoction. Combined extracts were filtered through 200-mesh sieves, concentrated under reduced pressure (50°C, 100 mbar) to relative density 1.10–1.15, and lyophilized to obtain dry powder.

Quality control employed high-performance liquid chromatography (HPLC) with photodiode array detection. Marker compound content in the final extract was: astragaloside IV (0.82±0.05 mg/g), timosaponin AIII (1.23±0.08 mg/g), platycodin D (0.56±0.04 mg/g), saikosaponin A (0.34±0.03 mg/g).

Batch-to-batch variation remained below 10% for all markers across three independent preparations. Before administration, lyophilized powder was reconstituted in distilled water at appropriate concentrations.

Hypobaric Hypoxia Model

HPH was induced using a hypobaric hypoxia chamber (DYC-300, Guizhou Fenglei Aviation Equipment Co., Ltd., Guizhou, China). The chamber simulated 5,000 m altitude conditions: barometric pressure 404 mmHg (approximately 54 kPa), ambient oxygen concentration ~10.8%.

Temperature and humidity within the chamber were maintained at 22±2°C and 50±10%, respectively. Carbon dioxide was scrubbed using soda lime.

Rats in hypoxia groups (Model, SXT-L, SXT-M, SXT-H) underwent daily 8-hour exposures (9:00–17:00) for 28 consecutive days. Ascent and descent rates were controlled at 10 m/s equivalent to minimize barotrauma. Control animals remained at ambient pressure (sea level equivalent, ~760 mmHg) in identical housing conditions within the same facility. Drug or vehicle (distilled water, 10 mL/kg) was administered by oral gavage daily at 8:00, one hour before hypoxia exposure commenced.

Hemodynamic Measurements and Right Ventricular Hypertrophy Assessment

After 28 consecutive days of intervention, rats were anesthetized by intraperitoneal injection of urethane (1 g/kg). Excess hair on the neck was removed, and a polyethylene catheter prefilled with heparinized saline was inserted into the pulmonary artery through the right jugular vein via a cervical incision. After stabilization of the pressure waveform, mean pulmonary artery pressure (mPAP) was recorded using a BL-420N biological signal acquisition system. Following hemodynamic measurements, the right ventricle (RV) was rapidly separated on ice from the left ventricle plus interventricular septum (LV + S). The right ventricular hypertrophy index (RVHI) and right ventricular weight index (RVWI) were calculated using the following formulas: RVHI = RV / (LV + S); RVWI= RV / body weight.

Histopathological Examination and Morphometric Analysis

After 28 days of treatment, rats were anesthetized with urethane (1 g/kg, intraperitoneal injection), and lung tissues were rapidly harvested on ice. Lung sections were stained with hematoxylin and eosin (H&E) and scanned using a digital slide scanner. Pulmonary arterioles with external diameters ranging from 30 to 80 μm were selected under microscopy, and vascular morphology was evaluated under a 400× optical microscope. Morphometric parameters included vascular wall thickness (WT), total vascular area bounded by the external elastic lamina (TA), external diameter (ED), and lumen area (LA). Pulmonary vascular remodeling indices were calculated according to the following equations: WT% = (2 × WT / ED) × 100%; WA% = (TA − LA) / TA × 100%.

Sample Collection and Processing

Twenty-four hours after the final hypoxia exposure, animals were fasted overnight for 12 h with free access to water. Rats were deeply anesthetized by isoflurane inhalation for euthanasia. Isoflurane (veterinary drug approval No. 153717015; Qingdao Obofang Pharmaceutical Technology Co., Ltd., Qingdao, China) was administered via inhalation until loss of righting reflex and absence of pedal withdrawal reflex were confirmed, ensuring adequate depth of anesthesia prior to exsanguination. Whole blood was then collected via abdominal aorta puncture using heparinized syringes.

Following blood collection, the pulmonary circulation was perfused via the right ventricle with 20–30 mL of ice-cold phosphate-buffered saline (PBS) at a constant flow rate of approximately 5 mL/min to remove residual blood before lung tissue harvest. Whole blood was allowed to clot at room temperature for 30 minutes, then centrifuged (3,000 rpm, 10 min, 4°C). Serum was aliquoted (200 µL per tube) and stored at −80°C until analysis. All serum samples were collected between 9:00 and 11:00 to minimize circadian variation.

After exsanguination under deep anesthesia, euthanasia was confirmed by cessation of heartbeat and respiration. Lungs were perfused with ice-cold PBS via the right ventricle to remove residual blood. Right lung lobes were dissected, flash-frozen in liquid nitrogen within 5 minutes of harvest, and stored at −80°C for RNA extraction. Left lung lobes were fixed in 4% paraformaldehyde for histopathological examination.

Untargeted Metabolomics Sample Preparation

Serum samples (100 µL) were thawed on ice and mixed with 400 µL cold methanol:acetonitrile (1:1, v/v) containing internal standards (2-chlorophenylalanine, 10 µg/mL). Mixtures were vortexed for 30 seconds, incubated at −20°C for 1 hour to precipitate proteins, then centrifuged (14,000 rpm, 15 min, 4°C). Supernatants were transferred to clean tubes, dried under nitrogen at 37°C, and reconstituted in 100 µL 50% acetonitrile. After centrifugation (14,000 rpm, 10 min, 4°C), supernatants were transferred to autosampler vials for analysis.

Quality control (QC) samples were prepared by pooling equal volumes (10 µL) from each study sample. QC samples were processed identically to study samples and injected at the beginning of the analytical run, after every 10 study samples, and at the end to monitor system stability and analytical drift.

UHPLC-Q-TOF-MS Analysis

Chromatographic separation employed a Waters ACQUITY UPLC system (Waters Corporation, Milford, MA, USA) equipped with a BEH C18 column (2.1 × 100 mm, 1.7 µm particle size). Mobile phases consisted of water containing 0.1% formic acid (A) and acetonitrile containing 0.1% formic acid (B). Gradient elution proceeded as follows: 0–2 min, 5% B; 2–11 min, 5–100% B; 11–13 min, 100% B; 13–13.5 min, 100–5% B; 13.5–16 min, 5% B for column re-equilibration. Flow rate was 0.4 mL/min, column temperature 40°C, and injection volume 5 µL.

Mass spectrometry employed a Waters Xevo G2-XS Q-TOF mass spectrometer with electrospray ionization (ESI) in both positive and negative modes. Source parameters: capillary voltage 2.5 kV (ESI+) or 2.0 kV (ESI-); cone voltage 40 V; source temperature 120°C; desolvation temperature 450°C; desolvation gas flow 800 L/h. Data were acquired in MSE mode (continuum) over m/z 50–1,000 with scan time 0.2 s. Leucine-enkephalin ([M+H]+ = 556.2771) served as lock mass for real-time mass accuracy correction.

Data Processing and Metabolite Identification

Raw data were processed using XCMS Online (Scripps Research Institute) for peak detection, alignment, and integration. Parameters: ppm=15, peakwidth=c(5,20), snthresh=6, mzwid=0.015. Features present in <80% of samples within any group were removed. Remaining features underwent probabilistic quotient normalization and log2 transformation. Metabolite annotation followed Metabolomics Standards Initiative (MSI) guidelines. Level 1 identification required matching retention time (±0.2 min) and MS/MS spectrum to authentic standards.

Level 2 identification required accurate mass (<10 ppm), isotope pattern matching, and MS/MS spectral similarity (>70%) to database entries in HMDB (version 5.0), KEGG, or METLIN. Features not meeting these criteria were assigned Level 3 annotations (putative compound classes based on accurate mass).

Because metabolomic and transcriptomic data were generated from the same animals, metabolite–gene associations were evaluated using paired sample matching based on individual animal identifiers.

For cross-omics integration, metabolite–gene associations were identified using a two-step strategy combining statistical correlation analysis and biological knowledge-based filtering. First, Spearman correlation analysis was performed between differential metabolites and differentially expressed genes using paired samples. Associations were considered significant when |r| > 0.6 and the Bonferroni-adjusted p value was < 0.001. Robust regression analysis was subsequently conducted as a sensitivity analysis to verify association stability. Second, statistically significant metabolite–gene pairs were retained only when biological support was available, including (i) mapping to the same KEGG pathway or (ii) documented enzyme–substrate, regulator–target, or other functional relationships curated in the KEGG or Reactome databases. Therefore, metabolite–gene links represented statistically validated and biologically supported associations and were interpreted as covariation relationships rather than direct causal regulatory interactions.

RNA Sequencing and Transcriptomics RNA Extraction and Quality Control

Total RNA was extracted from lung tissue (~50 mg) using TRIzol reagent according to the manufacturer’s protocol. RNA pellets were dissolved in RNase-free water and treated with DNase I (Qiagen; 1 U DNase I per 1 μg RNA) at 37 °C for 30 min to remove genomic DNA contamination, followed by enzyme inactivation according to the manufacturer’s instructions. RNA concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific), and samples with an A260/A280 ratio outside the range of 1.8–2.1 were excluded. RNA integrity was assessed using an Agilent Bioanalyzer 2100; only samples with a RNA Integrity Number (RIN) ≥ 7.0 were used for subsequent library construction. cDNA libraries were amplified by PCR for 12 cycles.

Library Preparation and Sequencing

Sequencing libraries were constructed using NEBNext Ultra II RNA Library Prep Kit (New England Biolabs) following manufacturer’s instructions. Briefly, mRNA was enriched using poly(A) selection with oligo(dT) magnetic beads, fragmented to ~200 bp, and reverse transcribed to cDNA. After end repair, A-tailing, and adapter ligation, libraries were amplified by PCR (12 cycles). Library quality was verified by Bioanalyzer (expected peak ~300 bp) and quantified by qPCR. Sequencing was performed on Illumina NovaSeq 6000 platform (paired-end 150 bp) targeting ~40–50 million reads per sample. Raw data quality was assessed using FastQC (version 0.11.9). Low-quality reads (Phred score <20) and adapter sequences were trimmed using Trimmomatic (version 0.39).

Read Alignment and Quantification

Clean reads were aligned to the rat reference genome (Rattus norvegicus Rnor_6.0, Ensembl release 104) using STAR aligner (version 2.7.9a) with default parameters. Alignment rates averaged 92.3±1.8% across samples. Gene-level read counts were generated using featureCounts (Subread package, version 2.0.1) with default parameters for paired-end data.

Differential Expression Analysis

Differential expression analysis employed DESeq2 (version 1.34.0) in R. Genes with total counts <10 across all samples were filtered. Normalization used DESeq2’s median-of-ratios method. Differentially expressed genes (DEGs) were defined by |log2 fold change| > 1 and Benjamini- Hochberg adjusted p-value (FDR) < 0.05. Volcano plots were generated using EnhancedVolcano package.

Bioinformatics Analysis Pathway Enrichment

Metabolic pathway enrichment analysis was performed in MetaboAnalyst 5.0 against the KEGG Rattus norvegicus pathway library. Enrichment significance was assessed by hypergeometric test with Holm-Bonferroni correction. Pathway impact scores were calculated based on pathway topology (relative betweenness centrality).

Transcriptomic pathway analysis employed clusterProfiler (version 4.2.2) for Gene Ontology (GO) and KEGG enrichment. Gene set enrichment analysis (GSEA) was performed using fgsea package (version 1.20.0) with KEGG and Hallmark gene sets from MSigDB. Significance thresholds: FDR < 0.05 for enrichment, |NES| > 1.5 for GSEA.

Weighted Gene Co-Expression Network Analysis

Weighted gene co-expression network analysis (WGCNA) (version 1.71) identified co-expressed gene modules. Soft-thresholding power (β=12) was selected based on scale-free topology criterion (R2 > 0.85). Minimum module size was set to 30 genes. Module eigengenes were correlated with treatment status using Pearson correlation. Hub genes within modules were identified by module membership (kME > 0.8) and gene significance.

Multi-Omics Integration

Metabolite-gene correlations were computed using Spearman’s rank correlation across all samples after z-score normalization of both datasets. To minimize false discoveries, pairs required |r| > 0.6 and Bonferroni-corrected p < 0.001. Robust regression using M-estimators served as sensitivity analysis against outlier-driven correlations. Integrated networks were visualized in Cytoscape (version 3.9.1). Network topology metrics (degree, betweenness centrality, closeness centrality) identified hub nodes. Network robustness was assessed by sequential node removal analysis.

Statistical Analysis

Continuous data are presented as mean ± standard deviation (SD). Normality was assessed by Shapiro–Wilk test. For normally distributed data, two-group comparisons used Student’s t-test; multi- group comparisons used one-way ANOVA with Tukey’s post-hoc test. For non-normally distributed data, Mann–Whitney U-test (two-group) or Kruskal–Wallis test with Dunn’s correction (multi-group) was applied. Multivariate analyses included principal component analysis (PCA) for unsupervised pattern recognition and orthogonal partial least squares-discriminant analysis (OPLS-DA) for supervised classification. OPLS-DA model validity was assessed by R2Y, Q2Y, and permutation testing (1,000 permutations). To statistically support PCA-based group separation, permutational multivariate analysis of variance (PERMANOVA) was performed on Euclidean distance matrices using the adonis2 function in the vegan package with 9,999 permutations. Pairwise comparisons were adjusted using the Benjamini–Hochberg method. Receiver operating characteristic (ROC) curves evaluated biomarker discriminatory performance. Area under the curve (AUC) and 95% confidence intervals were calculated. To reduce potential performance overestimation associated with the limited sample size, leave-one-out cross-validation (LOOCV) and five-fold cross-validation were performed to assess the robustness of the biomarker panel. All statistical analyses were performed in R (version 4.1.2) and Python (version 3.9.7). Significance was defined as p < 0.05 or FDR < 0.05 unless otherwise specified.

Results Validation of HPH Model and Therapeutic Effects of SXT on Hemodynamic and Right Ventricular Remodeling Parameters

To verify the successful establishment of the HPH model, mPAP, RVHI, and RVWI were measured. Compared with the Control group, the Model group exhibited a significant increase in mPAP (134.62%, p < 0.05), confirming successful induction of pulmonary hypertension. Following intervention, both XDNF and SXT treatment significantly reduced mPAP compared with the Model group. Specifically, mPAP was decreased by 27.74% (SXT-L), 26.71% (SXT-M), and 37.53% (SXT-H), while the XDNF group also showed a significant reduction in mPAP (p < 0.05) (Supplementary Figure 1A).

In addition, the RVHI in the Model group was significantly increased by 101.63% compared with the Control group (p < 0.05). High-dose SXT treatment (SXT-H) significantly reduced this index by 44.76% compared with the Model group (p < 0.05) (Supplementary Figure 1B).

Similarly, the RVWI was increased by 21.65% in the Model group relative to the Control group (p < 0.05), whereas SXT-H administration reduced this index by 17.10% compared with the Model group (p < 0.05) (Supplementary Figure 1C).

Histopathological examination of lung tissues (H&E staining) revealed normal alveolar architecture without obvious pathological injury in the Control group. In contrast, the Model group exhibited marked pulmonary structural damage, including necrosis and inflammatory cell infiltration. Morphometric analysis further showed significant increases in pulmonary arterial wall thickness percentage (WT%) and wall area percentage (WA%) compared with the Control group. SXT treatment alleviated these histopathological alterations in a dose-dependent manner, accompanied by significant reductions in WT% and WA%. Specifically, WT% was reduced by 40.16%, 50.12%, and 53.62% in the SXT-L, SXT-M, and SXT-H groups, respectively (p < 0.05). WA% was significantly decreased in the SXT-M and SXT-H groups by 26.43% and 26.04%, respectively (p < 0.05) (Supplementary Figure 2).

Collectively, these results confirm successful establishment of the HPH model and demonstrate that SXT effectively attenuates HPH and right ventricular remodeling, with overall trends suggesting dose-dependent histopathological improvement.

Metabolomics Data Quality and Global Profiling

UHPLC-Q-TOF-MS analysis detected 11,087 molecular features across all samples in combined positive and negative ionization modes. Analytical performance met quality benchmarks: >90% of features showed coefficient of variation (CV) <15% in pooled QC samples, internal standard CVs ranged from 5.2% to 8.7%, and retention time drift remained <0.1 min across the analytical batch. Missing values constituted <5% of the dataset and were imputed using k-nearest neighbor algorithm. Following QC-based filtering (removal of features with CV >30% in QC samples), 8,234 features proceeded to downstream analysis. Metabolite annotation achieved the following distribution: 2,156 features (26.2%) reached MSI Level 1 identification against authentic standards; 4,823 (58.6%) achieved Level 2 based on accurate mass, isotope pattern, and MS/MS library matching; the remainder fell into Level 3 category. Critically, the metabolites driving key findings—taurine, pyruvate, lactate, glutathione, and branched-chain amino acids—all received Level 1 confirmation, ensuring robust pathway interpretations. The annotated metabolome showed broad chemical diversity (Table 1). Lipid species predominated numerically (n=1,845), with phosphatidylcholines (687), sphingolipids (312), and lysophospholipids (198) representing major subclasses. Amino acids and derivatives (n=312) covered proteinogenic and non-proteinogenic species. Carbohydrate-related features (n=234) included hexoses, pentoses, and sugar phosphates integral to central carbon metabolism. Nucleotide derivatives (n=178) spanned both purine and pyrimidine classes.

Table 1 Top 20 Differential Metabolites Between Model and Control Groups

Principal component analysis revealed distinct metabolic phenotypes across experimental groups (Figure 1A). PC1 and PC2 together explained 41.2% of total variance (33.4% and 7.8%, respectively). Control and Model groups showed clear separation along PC1, confirming substantial metabolic perturbation induced by hypobaric hypoxia. SXT treatment groups occupied intermediate positions with dose-graded distribution: SXT-L clustered closer to the Model group, SXT-M occupied an intermediate position, and SXT-H showed a tendency toward the Control group rather than complete overlap, suggesting a partial reversal of disease-associated metabolic alterations rather than full normalization. This pattern is consistent with dose-associated metabolic shifts.

Figure 1 Principal component analysis (PCA) of the serum metabolome (A) and lung transcriptome (B). Samples are colored by group (Ctrl, green; Model, red; SXT-L/M/H, graded blues) and drawn as filled circles; open rings mark group centroids. Axes give the variance explained ((A) PC1 33.4%, PC2 7.8%; (B) PC1 31.5%, PC2 17.0%). Group structure was tested by PERMANOVA (adonis2, Euclidean distance, 9999 permutations): metabolome F = 2.19, R2 = 0.38; transcriptome F = 1.88, R2 = 0.41 (both p < 0.001), with significant pairwise contrasts after Benjamini–Hochberg correction. Model is displaced from Ctrl along PC1, with SXT groups at intermediate positions.

PERMANOVA analysis further confirmed significant differences in global metabolomic profiles among groups (F = 2.19, R2 = 0.38, p < 0.001), with significant separation between the Control and Model groups after Benjamini–Hochberg correction (adjusted p < 0.01). PERMANOVA further confirmed significant transcriptomic differences among groups (F = 1.88, R2 = 0.41, p < 0.001). OPLS-DA confirmed robust discrimination between Control and Model (R2Y=0.96, Q2Y=0.89). Permutation testing (1,000 iterations) yielded Q2 intercept of −0.21, validating model integrity against overfitting.

Transcriptomic Profiling Reveals Hypoxia-Induced Transcriptional Suppression

RNA sequencing generated 42.3±3.2 million reads per sample, with alignment rates averaging 92.3±1.8%. After filtering low-count genes, 14,582 genes remained for differential expression analysis. Comparison of Model versus Control identified 609 differentially expressed genes (DEGs) meeting criteria of |log2FC| > 1 and FDR < 0.05. Notably, downregulated genes (n=390, 64%) substantially outnumbered upregulated genes (n=219, 36%). This asymmetry suggests that chronic hypoxia operates predominantly through transcriptional suppression rather than activation—a finding with implications for understanding disease mechanism and therapeutic strategy.

Principal component analysis of transcriptomic data mirrored metabolomic patterns (Figure 1B). PC1 and PC2 captured 48.5% of variance (31.5% and 17.0%, respectively). Model samples formed a distinct cluster displaced from Control, while SXT groups showed graded intermediate positioning. Volcano plot visualization (Figure 2) highlighted the transcriptional landscape. Among the upregulated genes, canonical hypoxia-responsive transcripts predominated, including glycolytic enzymes such as hexokinase 2 (HK2), 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase 3 (PFKFB3), LDHA, PKM2, as well as hypoxia-inducible factor (HIF) target genes, including vascular endothelial growth factor A (VEGFA), BCL2 interacting protein 3 (BNIP3), and pyruvate dehydrogenase kinase 1 (PDK1), together with other angiogenesis-related factors. In contrast, downregulated genes were primarily associated with oxidative phosphorylation, including NADH oxidoreductase subunit A4 (NDUFA4), NADH oxidoreductase subunit B5 (NDUFB5), cytochrome c oxidase subunit 7A2 (COX7A2), and ATP synthase peripheral stalk subunit OSCP (ATP5O); fatty acid β-oxidation, including carnitine palmitoyltransferase 1A (CPT1A), acyl-CoA dehydrogenase medium chain (ACADM), and hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit alpha (HADHA); and mitochondrial biogenesis, including peroxisome proliferator-activated receptor gamma coactivator-1 alpha (PGC-1α), mitochondrial transcription factor A (TFAM), and nuclear respiratory factor 1 (NRF1).

Figure 2 Volcano plot of differentially expressed genes (Model vs Ctrl). Red, up-regulated in Model (n = 219); blue, down-regulated (n = 390); gray, not significant (thresholds: |log2FC| > 1.0, FDR < 0.05; dashed lines). Representative genes are annotated: up-regulated glycolytic/HIF-target genes (HK2, PFKFB3, LDHA, PKM2, VEGFA, PDK1, ADM) and down-regulated oxidative phosphorylation, fatty-acid β-oxidation and mitochondrial-biogenesis genes (NDUFA4, NDUFB5, COX7A2, ATP5O, CPT1A, ACADVL, HADHB, PGC-1α, TFAM, NRF1).

Coherent Metabolic Pathway Alterations

Differential metabolites revealed coherent pathway-level patterns. Branched-chain amino acids (BCAAs)—leucine, isoleucine, valine—declined 38–55% in Model serum (Table 1). This depletion may reflect either impaired hepatic synthesis, increased peripheral catabolism for energy production, or both. BCAA catabolism generates acetyl-CoA and succinyl-CoA that fuel the TCA cycle, suggesting mobilization to meet energy demands under hypoxic stress. The downstream implications extend to protein synthesis capacity, as BCAAs serve as both substrates and signaling molecules for mTOR-dependent translation.

Sulfur-containing amino acids showed modest declines: methionine (0.78-fold) and cysteine (0.72- fold). These species serve as precursors for glutathione synthesis, and their depletion alongside reduced GSH levels (0.52-fold) and elevated GSSG/GSH ratio points to glutathione system stress. Taurine, another sulfur-containing compound with antioxidant properties, showed the most dramatic reduction (0.45-fold)—consistent with consumption under oxidative stress. Energy metabolism intermediates painted a coherent picture of mitochondrial dysfunction. Pyruvate accumulated 2.31-fold while TCA cycle intermediates—citrate (0.56-fold), succinate (0.62-fold), α-ketoglutarate (0.74-fold), malate (0.71-fold)—were uniformly depleted. This pattern indicates a bottleneck at mitochondrial pyruvate entry: glycolysis operates but oxidative metabolism stalls. The elevated lactate (1.87-fold) and lactate-to-pyruvate ratio corroborate this interpretation, reflecting increased cytosolic NADH/NAD⁺ ratio characteristic of anaerobic glycolysis.

Arginine metabolism alterations carry vascular implications. Arginine declined 35% (0.65-fold), with downstream metabolites showing divergent changes: citrulline decreased (0.76-fold) while ornithine increased (1.35-fold). This pattern suggests arginase pathway activation—diverting arginine from nitric oxide synthesis toward polyamine production. Reduced NO bioavailability would impair endothelium-dependent vasodilation, contributing to sustained vasoconstriction.

SXT restored these metabolic alterations in dose-dependent fashion. At highest dose, pyruvate (0.92-fold versus Model), lactate (1.08-fold), and TCA intermediates approached Control levels.

Taurine and GSH showed partial restoration (67–78% of Control). The pattern suggests SXT addresses the metabolic bottleneck at mitochondrial pyruvate utilization while bolstering antioxidant capacity.

Pathway Enrichment Analysis

Metabolomic pathway enrichment (Figure 3) identified several significantly affected pathways. Top hits included: aminoacyl-tRNA biosynthesis (p=0.014, impact=0.42), which links amino acid availability to protein synthesis capacity; glycolysis/gluconeogenesis (p=0.018, impact=0.52), reflecting the glycolytic shift; taurine/hypotaurine metabolism (p=0.032, impact=0.38), indicating antioxidant pathway stress; thiamine metabolism (p=0.028, impact=0.31); and glutathione metabolism (p=0.039, impact=0.45).

Figure 3 KEGG pathway enrichment of the differential metabolites. (A) Lollipop plot ranked by significance (x-axis, −log10 p; dashed line, p = 0.05); dot size denotes the number of mapped metabolites (hits), and pathways with fewer than three hits are shaded and interpreted with caution. (B) Bubble plot of pathway impact (x-axis) versus significance (y-axis, −log10 p), with color encoding −log10 p. Key pathways are labeled with p-value and impact: aminoacyl-tRNA biosynthesis (p = 0.014, impact = 0.42), glycolysis/gluconeogenesis (0.018, 0.52), thiamine metabolism (0.028, 0.31), taurine and hypotaurine metabolism (0.032, 0.38), and glutathione metabolism (0.039, 0.45).

Transcriptomic GSEA revealed concordant pathway alterations. Positively enriched gene sets (upregulated in Model) included: hypoxia response (NES=2.34, FDR<0.001), glycolysis (NES=1.98, FDR=0.002), HIF-1 signaling (NES=1.85, FDR=0.004), and angiogenesis (NES=1.72, FDR=0.008).

Negatively enriched gene sets (downregulated in Model) included: oxidative phosphorylation (NES=- 2.56, FDR<0.001), fatty acid metabolism (NES=−1.92, FDR=0.003), TCA cycle (NES=−1.78, FDR=0.006), and mitochondrial biogenesis (NES=−1.65, FDR=0.012).

The concordance between metabolomic and transcriptomic pathway enrichments strengthens confidence that observed changes represent genuine pathway-level rewiring rather than isolated fluctuations. Eight pathways showed dual-layer enrichment: glycolysis/gluconeogenesis, oxidative phosphorylation, glutathione metabolism, arginine/proline metabolism, taurine metabolism, TCA cycle, pyruvate metabolism, and BCAA degradation. These cross-platform convergence points represent high-confidence targets of both hypoxia pathophysiology and SXT therapeutic action.

Multi-Omics Network Integration

Integration of metabolomic and transcriptomic data generated a correlation network comprising 47 nodes (23 metabolites, 24 genes) and 86 edges (Figure 4). Each edge represents Spearman correlation |r| > 0.6 with Bonferroni-corrected p < 0.001. Network visualization employed force- directed layout to reveal modular structure.

Figure 4 Multi-omics correlation network (force-directed layout; 47 nodes, 86 edges). Pink nodes, metabolites (n = 23); blue nodes, genes (n = 24); node size reflects degree centrality. Edges represent Spearman |r| > 0.6 (green, positive; red, negative) supported by shared pathway membership; edge direction does not imply causal regulation. The network comprises an energy-metabolism core (metabolite hub Pyruvate; gene hubs PKM2 and LDHA) and an amino-acid module, bridged by α-ketoglutarate and glutamate. Only hub and bridge nodes are labeled.

Topology analysis revealed densely connected energy metabolism core. Pyruvate emerged as the most connected metabolite hub (degree=12), with edges linking to glycolytic enzymes (HK2, PFKFB3, LDHA, PKM2), TCA cycle genes (CS, IDH1, SDHA), and gluconeogenic regulators (PCK1, FBP1).

Amino acid metabolism formed a secondary module with BCAAs clustering with their catabolic enzymes (BCAT1, BCKDHA, BCKDHB). Bridge metabolites—α-ketoglutarate and glutamate — connected the two modules, reflecting their dual roles in carbon and nitrogen trafficking.

The network exhibited small-world properties: clustering coefficient 0.42 (significantly higher than random networks of equal size, p<0.001) and average path length 2.3 edges. These characteristics indicate efficient information flow alongside vulnerability to hub perturbation.

In silico robustness analysis confirmed hub criticality. Sequential removal of pyruvate fragmented the network into three disconnected components. LDHA removal produced similar fragmentation. By contrast, random node removal required deleting >40% of nodes to achieve comparable fragmentation. This differential vulnerability underscores the structural importance of identified hubs.

On the transcriptome side, PKM2 (degree=9) and LDHA (degree=8) occupied analogous hub positions. Both enzymes sit at critical control points of glycolytic flux: PKM2 controls phosphoenolpyruvate-to-pyruvate conversion and exists in regulatable tetrameric/dimeric forms; LDHA catalyzes pyruvate-to-lactate interconversion. Their network centrality, combined with extensive metabolite correlations, positions them as prime therapeutic targets.

Dose-Response Relationships and Expression Patterns

Heatmap visualization of top 50 DEGs (Figure 5) illustrates transcriptional patterns across groups. Hierarchical clustering segregated samples by treatment condition. Model animals exhibited characteristic signature: upregulation of hypoxia-responsive/glycolytic genes (Cluster 1), downregulation of oxidative metabolism genes (Cluster 2). SXT treatment produced graded normalization—low-dose animals clustered nearer Model, medium-dose occupied intermediate positions, high-dose approached Control.

Figure 5 Heatmap of the top 50 differentially expressed genes (Model vs Ctrl; n = 50), ranked by adjusted p-value. Color shows z-score-normalized expression (blue, low; red, high). Genes (rows) and samples (columns) were ordered by hierarchical clustering (dendrograms shown), and the top bar indicates treatment group (Ctrl, Model, SXT-L/M/H). Two gene clusters are resolved: C1, hypoxia-responsive/glycolytic genes up-regulated in Model, and C2, oxidative-metabolism genes (oxidative phosphorylation, TCA cycle and fatty-acid oxidation) down-regulated in Model; SXT shifts both clusters toward the Ctrl pattern.

Dose Response Index (DRI) quantified dose-dependence: DRI = |Effect_High - Effect_Model| / |Effect_Control - Effect_Model|. Values range from 0 (no restoration) to 1 (complete restoration to Control level). ANOVA confirmed significant dose-dependence (p<0.001) for 78% of differential metabolites and 65% of DEGs.

Taurine exemplified strong dose-dependency (DRI=0.89): levels progressed from 45% of Control in Model to 67% (SXT-L), 78% (SXT-M), 92% (SXT-H). The dose-response curve fit a sigmoidal model (R2=0.94) with ED50 ≈ 4.2 g/kg. Pyruvate (DRI=0.76) and lactate (DRI=0.58) showed similar patterns.

Among genes, glycolytic enzymes (HK2, PFKFB3, LDHA) showed strong dose-dependent downregulation (DRI=0.72–0.85), while mitochondrial genes (NDUFB5, COX7A2) showed dose- dependent upregulation (DRI=0.65–0.78).

At high dose, 89% of differential metabolites shifted toward Control. A plateau effect emerged: SXT-M achieved ~60–70% of maximal response seen at SXT-H. This observation carries practical implications—moderate dosing may deliver substantial benefit with potentially improved tolerability.

Candidate Biomarker Identification

ROC analysis evaluated discriminatory performance of candidate biomarkers (Table 2). Taurine achieved highest single-metabolite AUC (0.94, 95% CI: 0.87–0.99) for distinguishing Model from Control. Pyruvate (AUC=0.91) and glutamine (AUC=0.89) also demonstrated strong discrimination. Each captures a distinct pathophysiological facet: taurine reflects antioxidant status, pyruvate indicates glycolytic flux, glutamine reports anaplerotic capacity.

Table 2 Biomarker Performance Characteristics

Combination of these three metabolites into a multi-marker panel further improved performance, yielding an apparent AUC of 0.98 (95% CI: 0.95–1.00). To