A total of 113 Chinese Han individuals (104 males, 9 females), aged between 18 and 63, were enrolled in this cross-sectional study. Among them, 56 had experienced or were experiencing acute mountain sickness (eAMS+), while 57 did not. DAHA ranged from 1 to 23 days (Supplementary Figure S1; Table 1). Participants were categorized into four groups based on DAHA: T1 (≤3 days, n = 22), T2 (4–7 days, n = 18), T3 (8–14 days, n = 29), and T4 (15–23 days, n = 44). The Kruskal–Wallis test revealed no significant demographic differences (sex, age, BMI, eAMS) among the four DAHA groups (Table 2). Chi-square test or t-test suggested no demographic differences among eAMS + group and eAMS- group (Table 3).

We analyzed 44 physiological parameters (Supplementary Table S1) to define acclimatization traits. After excluding sample S187 for excessive missing data (17 out of 44 items, 38.64%), principal component analysis (PCA) was performed on 112 individuals. Four additional outliers (S381, S472, S15, S13) out of 95% confidence intervals were removed (Figures 1a,b). Multiple linear regression model (adjusted for age, sex, and BMI) identified four parameters exhibiting a significant linear relationship with DAHA after adjusting for multiple comparisons via the Benjamini–Hochberg method (FDR <0.05). These parameters included saturation of peripheral oxygen (SpO₂, r = 0.331), hemoglobin (HGB, r = 0.602), hematocrit (HCT, r = 0.553), and the standard deviation of red blood cell distribution width (RDW-SD, r = 0.46) (Figures 1c–f; Supplementary Table S2). However, only RDW-SD remained significant in the subsequent Kruskal‒Wallis test (FDR = 0.0066, Supplementary Table S3).

To explore the relationship between GEPs and DAHA, we performed transcriptomics on a subset of 48 individuals with complete RNA sequencing data (including 35 eAMS+ and 13 eAMS- participants). Their DAHA ranged from 2 to 19 days and were divided into T1 (≤3 days, n = 4), T2 (4–7 days, n = 11), T3 (8–14 days, n = 17), and T4 (15–19 days, n = 16) (Supplementary Figures S2a,b). Using Fuzzy C-Means Clustering (FCM) clustering algorithm, genes were grouped into 12 clusters with distinct temporal dynamics. Among these, Cluster 4 (C4, n = 1,389) and Cluster 9 (C9, n = 1,469) showed consistently increasing and decreasing gene expression over time, respectively (Figure 2a; Supplementary Figure S2c). C4 was enriched in adaptive immune response, nonsense-mediated decay independent of the exon junction complex, steroid biosynthesis, and the regulation of protein transport. And C9 was linked to the regulation of the response to biotic stimuli, the innate immune response, neutrophil degradation, interferon gamma signaling, the regulation of immune effector processes, and other immune or inflammatory pathways (Supplementary Figures S2d,e).

Moreover, differentially expressed gene (DEG) analysis across DAHA groups revealed distinct gene changes in each period: P1 (T2 vs. T1, 681 up, 169 down), P2 (T3 vs. T2, 205 up, 108 down), and P3 (T4 vs. T3, 217 up, 486 down) (p < 0.05, |logFC| > 0.3785) (Figure 2b). Hallmark Gene Set Enrichment Analysis (HGSEA) revealed that heme metabolism was upregulated in P1 then downregulated in P2 and P3, aligned with, the “acute increase followed by decrease” expression trend of Cluster 8 (C8). Interferon alpha/gamma response were not different in P1 but significantly downregulated in P2 and P3 (p < 0.05, q < 0.25, and |Normalized Enrichment Score (NES)| > 1) (Figure 2c). C8 was specially enriched in heme metabolism (p = 1e-10, q = 4e-9, NES = 2.13), while the interferon alpha/gamma response was enriched in C9. In contrast, no significant HGSEA results were identified for C4 (Supplementary Table S4). Based on the expression trends and enrichment results, the functional characteristics of C4, C8, and C9 were further compared.

We performed WGCNA on the transcriptomic data to investigate gene modules associated with DAHA. After excluding 3 outlier samples identified by hierarchical clustering identified, a soft threshold power of 6 was selected for network construction (Supplementary Figures S3a,b). The gene set was divided into modules using an adjacency matrix, which quantifies the pairwise correlations with genes. To ensure robust module definitions, modules with over 75% similarity were merged, resulting in 26 distinct modules (Supplementary Figure S3c). Among these, the darkorange module (DO) presented a significant negative correlation with DAHA (r = −0.3, p = 0.048). Conversely, the darkmagenta module (DM) were positively correlated (r = 0.39, p = 0.0086). Significant correlations with eAMS were identified: positive for the blue module (BL; r = 0.35, p = 0.019) and negative for DO (r = −0.34, p = 0.021) (Supplementary Figure S3d).

Enrichment analysis showed that DO-associated genes (n = 114) were linked to several key immune processes such as the antiviral defense response and interferon response, sharing similar expression trends and functions with C9. In contrast, the DM (n = 37) were involved in memory and ion transport, while BL (n = 1911) was significantly enriched in protein catabolic, cellular stress responses, metabolic processes, and autophagy regulation (Supplementary Figure S4).

Using integrated analysis of time series and WGCNA, we identified 18 high-confidence genes shared between the DM module and C4, both positively associated with the DAHA (Figure 3a). These genes were associated primarily with the negative regulation of hydrolase activity (3 genes), inorganic ion transmembrane transport (4 genes), and intracellular chemical homeostasis (3 genes) (Supplementary Figure S5a). Although 15 of these genes were found in the STRING database, protein‒protein interaction (PPI) network analysis revealed no significant associations among them (Supplementary Figure S5b).

The gene set with a decreasing expression trend with DAHA, the DO module and C9 shared 90 immunity-related genes involved in antiviral defense, interferon responses, and tuberculosis immunity. To pinpoint key genes, we imported the PPI results into Cytoscape. The top 5% (n = 5) of genes with the highest degree centrality were selected as hub genes: STAT1, OAS1, OAS2, RSAD2, and IFI44. Considering the consistent expression trend and functional profile of this intersection gene set, we designated it as Pattern 1 (Figures 3a,b,d).

The BL module shared 85 genes with C4, enriched in pathways including the Nop56p-associated pre-rRNA complex (4 genes), insulin receptor signaling pathway (3 genes), cholesterol metabolic process (3 genes), positive regulation of the MAPK cascade (5 genes) and positive regulation of translation (3 genes) (Supplementary Figure S5c). It also shared 86 genes with C9, which were linked to organic acid catabolic processes (4 genes), proteinogenic amino acid metabolic processes (3 genes), and herpes simplex virus 1 infection (3 genes) (Supplementary Figure S5d).

For C8, the BL module presented the highest gene overlap. These genes were enriched in processes like protein ubiquitination, erythrocyte differentiation, autophagy, and mitophagy. The top 20 hub genes were identified via PPI network analysis. These genes comprised UBA52, UBB, FECH, SNCA, SLC4A1, GATA1, RBX1, ALAS2, SKP1, EPB42, GYPB, KLF1, PSMD4, ANK1, CDC34, PINK1, BCL2L1, DMTN, DCAF12, and GYPA. The “mountain-shaped” expression pattern of these 818 genes was defined as Pattern 2 (Figures 3a,c,e).

Furthermore, we investigated the correlations between the hub genes and the 44 physiological parameters mentioned above. Hub genes in Pattern 1 were negatively correlated with plateletcrit (PCT) and platelets (PLTs). Key correlations included STAT1 negative with the percentage of monocytes (MONOP) and mean corpuscular hemoglobin concentration (MCHC). OAS1 positive with creatinine (CREA) and aspartate aminotransferase (AST), OAS2 positive with HCT and HGB, and IFI44 weakly positive with the white blood cell (WBC) count.

In Pattern 2, nearly all hub genes correlated negatively with immune cells like WBC, monocyte count (MONO), MONOP, neutrophil count (NEUT), lymphocyte count (LYMPH) and hematological parameters like HCT, red blood cell count (RBC) and HGB, but were positively correlated with CREA. UBA52, the hub gene with the highest degree in Pattern 2, exhibited positive correlation with the percentage of eosinophils (EOP), aspartate aminotransferase/alanine aminotransferase (ASL/ALT), and direct bilirubin (DBIL), and negative correlation with the platelet distribution width (PDW), low-density lipoprotein cholesterol (LDLC), and total cholesterol (TC). Notably, PSMD4 presented the strongest positive correlation with CREA among all Pattern 2 hub genes (Supplementary Figure S6).

WGCNA revealed distinct gene network expression in DO and BL between the eAMS+ and eAMS- groups (Supplementary Figure S3d). Using a general linear model with AMS experience as the independent variable and DAHA, age, sex and BMI as the covariates, 583 upregulated genes and 104 downregulated genes were identified (Supplementary Table S5). Enrichment analysis showed that the upregulated genes were involved in adaptive immune response activity, erythrocyte differentiation and the CAMKK2 pathway (Supplementary Figure S7a). Conversely, the downregulated genes were linked to leukocyte-mediated cytotoxicity and antibacterial humoral response (Supplementary Figure S7b).

We identified 398 upregulated genes and 10 pathways shared with Pattern 2, such as gas transport, erythrocyte differentiation, autophagy regulation and so on. In addition, only one downregulated gene overlapped with Pattern 1, though both showed enrichment in immune response and interferon gamma signaling (Figures 4a,b). The HGSEA results further revealed upregulation of heme metabolism and xenobiotic metabolism pathway, and downregulation of interferon response and inflammatory pathways in the eAMS + group (Figure 4c). Heme metabolism was also found in Pattern 2, while the interferon responses were common in Pattern 1 (Figure 4d). These findings highlight molecular distinctions between eAMS+ and eAMS- groups and their strong association with both Pattern 1 and Pattern 2.

To contextualize our findings, we integrated our Pattern 2 data with the public dataset GSE75665, which tracks transcriptomic changes in human blood before and after 3-day high-altitude exposure. After high-altitude exposure, 679 genes were upregulated and 1,020 downregulated compared with their baseline conditions. Of these, 211 upregulated genes overlapped with Pattern 2 in our study. Enrichment analysis revealed that the gene set shared was primarily associated with proteolysis, erythrocyte differentiation, porphyrin metabolism, and the regulation of autophagy and protein catabolism (Figures 5a,b).

At low altitudes, 64 of the 182 DEGs between AMS and non-AMS groups overlapped with Pattern 2. These included five Pattern 2 hub genes (BCL2L1, DCAF12, CDC34, PINK1 and UBB), highlighting their potential as predictive biomarkers for AMS (Figure 5c). After high-altitude exposure, the overlapped genes with Pattern 2 increased to 128 DEGs with eleven hub genes (ALAS2, BCL2L1, DCAF12, CDC34, EPB42, GATA1, PINK1, PSMD4, SNCA, UBA52, and UBB) including the five initial biomarkers (Figure 5d). Conversely, in the GSE75665 data, Pattern 1 was nearly absent, found in none of the 182 low-altitude DEGs and only 2 high-altitude DEGs (EIF2AK2 and TRIM22) (Figure 5c; Supplementary Tables S6,S7).

The study was designed as a cross-sectional investigation of physiological and molecular characteristics associated with high-altitude exposure. A total of 113 Chinese Han individuals who recently arrived at a high altitude (LiTang, 4,104 m) were enrolled between February 29th and 6 March 2024. And those with diabetes, respiratory disease, cardiovascular or neurological diseases were excluded. All participants underwent a single standardized physical examination and blood sampling during their stay at high altitude. Meanwhile, participants completed structured case report form (CRF) questionnaires to document demographic characteristics and AMS symptoms occurring within the first 3 days after ascent. Throughout the high-altitude exposure period, they received standardized meals with comparable nutritional contents and were living and working in the same environment. The study was approved by the Ethics Committee of West China Hospital (No. 2021–716). All participants provided written informed consent for the use of their samples and data for research purposes.

The symptoms of AMS include headache, gastrointestinal issues, fatigue and/or weakness, and dizziness/light-headedness. AMS was diagnosed via the 2018 Lake Louise Acute Mountain Sickness Score (LLAMS) (Roach et al., 2018), which identifies individuals with a headache score of at least one point and a cumulative score of at least 3 points. For participants who examined within the first 3 days after ascent, AMS symptoms were recorded prospectively. For participants who had already resided at high altitude for more than 3 days at the time of assessment, AMS symptoms occurring within the first 3 days after ascent were recorded retrospectively based on participant recall, using a structured questionnaire administered by trained investigators. Participants who experienced or were experiencing AMS during the early post-ascent period were defined as the eAMS + group, and individuals who did not meet these criteria were classified as eAMS negative (eAMS-) group, which served as the internal control group.

Blood was collected from a peripheral vein after 12 h of fasting at high altitude. As this was a cross-sectional study, blood sampling was performed once at a single time point for each individual. The duration of high-altitude exposure at the time of blood collection was retrospectively determined based on each participant’s documented arrival date at high altitude, resulting in exposure durations ranging from 1 to 23 days.

Blood samples collected for complete blood cells (CBC) count and blood biochemical tests were immediately sent to LiTang County People’s Hospital for testing. A subset of 48 participants, selected from the overall cohort, provided samples for RNA sequencing. All blood samples of the RNA-sequencing set were stored at −80 °C with PAXgene Blood RNA Tubes (PAXgene, PreAnalytix, Hombrechtikon, Switzerland, distributed by Qiagen, catalog no. 762165) until subsequent RNA extraction and sequencing.

Total RNA was extracted from blood samples collected from a subset of 48 participants selected from the original cohort of 113 individuals, using lysis, PAXgene spin columns, and DNase treatment. After quality assessment, mRNA was enriched, globin mRNA removed, and fragmented for cDNA synthesis. Following library preparation with dUTP-based second-strand synthesis, end repair, A-tailing, and adapter ligation, DNA nanoballs (DNBs) were generated and sequenced on the BGI-T7 platform (BGI-Shenzhen, China) using combinatorial probe-anchor synthesis (cPAS) to achieve high-throughput pairwise sequencing.

The raw sequencing data were filtered via SOAPnuke (Chen et al., 2018) (version 1.6.5) to remove the following reads: 1) reads containing adapters (adapter contamination), 2) reads with an unknown base N content exceeding 1%, and 3) low-quality reads (where the proportion of bases with a quality score below 15 exceeded 40% of the total read length). Clean reads were aligned to the Homo sapiens reference genome (GRCh38.p14) via STAR (Dobin et al., 2013) (version 2.7.10b) and gene expression counts were quantified via featureCounts (Liao et al., 2014) (version 2.0.4). Normalization was calculated via the trimmed mean of M values (TMM) method in the R package edgeR (Robinson et al., 2010) (version 4.0.16). Differentially expressed genes (DEGs) between the eAMS+ and eAMS- groups were identified using Generalized Linear Model (GLM) (adjusted for DAHA, age, sex and BMI) with significance defined as p < 0.05 and |logFC| > 0.3785.

Based on the duration of high-altitude acclimatization (Lucas et al., 2011; Gaur et al., 2020), the subset of 48 participants for RNA sequencing analysis had exposure duration ranging from 2 to 19 days (Supplementary Figure S2a). These participants were divided into 4 time-point groups (T1: ≤3 days; T2: 4–7 days; T3: 8–14 days; T4: 15–19 days) (Supplementary Figure S2b). Time series gene expression analysis was performed using FCM algorithm from the R package Mfuzz (Kumar and E Futschik, 2007) (version 2.62.0) to categorize genes into 12 clusters on the 4 time-point groups. Furthermore, Weighted Co-Expression networks were constructed using the R package WGCNA (Langfelder and Horvath, 2008) (version 1.73) to identify the relationships among genes and investigate their associations with specific traits.

Gene enrichment analysis was conducted via the Metascape (Zhou et al., 2019) (version 3.5) (https://metascape.org/) with default settings by inputting unique gene symbols from each gene set of interest. Hallmark Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) was utilized to identify significant pathways (p < 0.05, q < 0.25, |NES| > 1) between the AMS and non-AMS groups. PPI network from STRING (Szklarczyk et al., 2023) (version 12.0) (https://string-db.org/) was visualized in Cytoscape (Shannon et al., 2003) (version 3.10.3). Hub genes were defined as those in the top 2.5% or 5% of node degree based on the size of the gene set.

For compensating the lack of baseline data and prediction biomarker identification, we conducted a systematic literature search for studies associated with AMS. Inclusion criteria were: 1). The transcriptomics data is available; 2). The participants should have baseline data and AMS outcomes following high-altitude exposure; 3). As genetic and environmental heterogeneity significantly influences gene expression patterns (Martin et al., 2014), we restricted the baseline data to the Chinese population. Following our inclusion criteria, only one public transcriptome data GSE75665 (Liu et al., 2017) was included. This dataset comprised 20 samples divided into two groups (AMS and non-AMS) and two time points (baseline and 72 h after high-altitude exposure). DEGs based on longitudinal groups were analyzed via a density-based pruning algorithm. For the baseline cross-sectional comparison (AMS vs. non-AMS = 5 vs. 5), we employed the same analytical approach, excluding the lowest 10% expressed genes and considering |logFC| > 0.3785 as statistically significant.