The human blood transcriptome exhibits time-of-day-dependent response to hypoxia: Lessons from the highest city in the world.

High altitude exposes humans to hypobaric hypoxia, which induces various physiological and molecular changes. Recent studies point toward interaction between circadian rhythms and the hypoxic response, yet their human relevance is lacking. Here, we examine the effect of different high altitudes in conjunction with time of day on human whole-blood transcriptome upon an expedition to the highest city in the world, La Rinconada, Peru, which is 5,100 m above sea level. We find that high altitude vastly affects the blood transcriptome and, unexpectedly, does not necessarily follow a monotonic response to altitude elevation. Importantly, we observe daily variance in gene expression, especially immune-related genes, which is largely altitude dependent. Moreover, using a digital cytometry approach, we estimate relative changes in abundance of different cell types and find that the response of several immune cell types is time- and altitude dependent. Taken together, our data provide evidence for interaction between the transcriptional response to hypoxia and the time of day in humans.

Introduction

High-altitude exposure (i.e., over 2,500 m) leads to hypobaric hypoxia, which gradually increases with elevation (Luks, 2015; West, 2012). The primary and acute response aims to preserve sufficient tissue oxygenation through activation of the sympathetic nervous system and increase in ventilation, heart rate, and blood flow (Bartsch and Gibbs, 2007; Luks, 2015). Subsequently, the hypoxic response involves the activation of gene-expression program that is mainly regulated by the hypoxia-inducible factor (HIF) pathway. One of its hallmarks is an increase in hemoglobin mass due to renal hypoxemia, HIF activation, and erythropoietin release (Luks, 2015; West, 2012). Despite the wide knowledge on the HIF gene program activation (Schodel and Ratcliffe, 2019; Semenza, 2012), the overall transcriptional response induced by hypobaric hypoxia secondary to high-altitude exposure in humans is relatively uncharted (Yasukochi et al., 2020).

Himalayan and Andean highlanders, living for thousands of years at high altitude, have acquired progressive genomic adaptations that allow them to cope with a hypoxic environment (Azad et al., 2017; Bigham and Lee, 2014; Pamenter et al., 2020). In contrast, lowlanders acutely exposed to high altitude are susceptible to develop acute mountain sickness with severe complications such as high-altitude pulmonary or cerebral edema, which reflect poor acclimatization to hypoxia (Hackett and Roach, 2001). Hence, the mechanisms that are implicated in proper acclimatization of lowlanders to hypoxia are of interest in respect to high-altitude exposure and are potentially relevant for hypoxic conditions associated with various pathologies (Grocott et al., 2007).

Recent studies, mostly with cultured cells and animal models, point toward inter-connections between hypoxia signaling pathways and the circadian clock (Adamovich et al., 2022; O'Connell et al., 2020; Peek, 2020). Specifically, we showed that physiological oxygen rhythms can synchronize circadian clocks and that acute hypoxia can induce circadian misalignment between peripheral clocks in mice (Adamovich et al., 2017; Manella et al., 2020). Concurrently, the circadian clock was found to regulate oxygen rhythmicity in mice and temporally gate the response to hypoxia (Adamovich et al., 2019; Manella et al., 2020; Peek et al., 2017; Wu et al., 2017). In humans, several studies indicated that high altitude can influence physiological daily rhythms, including body temperature and heart rate (Bosco et al., 2003; Mortola, 2007; Ponchia et al., 1994; Vargas et al., 2001). Yet, little is known at the molecular level (i.e., gene expression) on the time-of-day effects of low oxygen availability in humans.

To examine the effect of exposure to hypoxia on gene expression, as well as its potential interaction with the time of day, we collected data during the 5300 Expedition (Champigneulle et al., 2021; Hancco et al., 2020; Oberholzer et al., 2020). This expedition mostly studies the health consequences of the extreme environmental conditions on the indigenous habitants of the highest city in the world, La Rinconada, Peru, located 5,100 m above sea level in the Peruvian Andes. In the present field study, which was conducted under an extremely challenging experimental setting in one the most remote places on earth, we examined the effect of high altitude and its interaction with time of day on the expedition members, who were all lowlanders. The subjects’ blood was sampled at three different altitudes (i.e., 200, 3,800, and 5,100 m) around the clock. Whole-blood RNA transcriptome analyses showed that exposure to high altitude carries substantial effects on gene expression. Moreover, our data revealed that the daily variance in expression of a subset of genes is influenced by altitude. These genes were enriched for different innate immunity and erythropoietic processes and, alongside our digital cytometry analysis, suggested that several immune cell types respond in a time- and altitude-dependent manner.

Results

The effects of high altitude on physiological parameters

In order to test the effects of high altitude on various physiological and molecular parameters in humans, 8 male subjects were tested at three different altitudes: ∼3,800 m at the city of Puno, Peru (PU), ∼5,100 m at the town of La Rinconada, Peru (RI), and, as a baseline reference, at ∼200 m above sea level (SL) in Grenoble, France. Whole-blood samples were collected from every subject at 4 h intervals throughout 24 h in each altitude. Sampling and measurements were performed on the seventh day in either PU or RI (Figure 1A). All subjects spent a pre-acclimatization period of 10 days at high altitude before the week in PU to minimize effects related to acute response to high altitude as well as other potential confounding effects that might stem from the nature of this field study (e.g., jet lag). All subjects were healthy adult males (35.3 ± 10.8 years old) with normal body mass indices (BMIs; 21.56 ± 0.78), (for anthropometric details, see Table S1).

Figure 1

Effects of high altitude on physiological parameters

(A) The study design. Data were collected from 8 male subjects at three different altitudes (SL, sea level, 200 m; PU, Puno, 3,800 m; RI, La Rinconada, 5,100 m) according to the scheme. In each 24-h sampling round, blood samples were collected from each subject every 4 h for subsequent transcriptomic analysis. In addition, subjects were continuously monitored for activity and heart rate (HR). Blood saturation was measured at 4-h intervals alongside the blood collections. All other physiological parameters were monitored in each altitude once a day.

(B) Blood saturation (SpO2), PetO2, and PetCO2 measurements.

(C) Hemoglobin mass per body weight (BW), hemoglobin concentration ([Hb]), and hematocrit measurements. The measurements in (B) and (C) were performed during the sampling round in each altitude, between 10:00 a.m. and 2 p.m.

(D) Daily mean HR and light-intensity activity levels during the rest phase. Longitudinal HR and activity data are available in Figure S1.

(E) Rest phase (clock time 10 p.m. to 6 a.m.) mean HR and light-intensity activity levels.(∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 ∗∗∗∗p < 0.0001, ns, non-significant, two-sided paired t test, n = 8 subjects). Each line represents measurements from a single subject (B–E). Note, for technical reasons, SL Hb mass, Hb measurements, and hematocrit for one individual are lacking.

We first measured the effects of the different altitudes on various physiological parameters. These include blood O2 saturation (SpO2) and the end-tidal partial pressure of O2 and CO2 (PetO2 and PetCO2, respectively), which are considered as reliable proxies for arterial blood O2 and CO2 concentrations, respectively (Ramos et al., 2016). As expected, we observed a decrease in blood O2 with altitude as measured by both SpO2 and PetO2 (Figure 1B). PetCO2, however, only showed a significant decrease from SL to PU (Figure 1B). Of note, SpO2 levels were monitored at 4 h intervals during the blood sampling and were relatively constant throughout the day (JTK_CYCLE, non-significant for 24 h rhythms) (Figure S1A). Hemoglobin mass and concentration, as well as hematocrit, were elevated in the higher altitudes, namely PU and RI (Figure 1C). Nevertheless, we did not observe a further increase in these parameters from PU to RI. Hence, it appears that upon elevation from PU to RI, the reduction in oxygen availability was not accompanied by a further increase in hemoglobin content.

The participants were equipped with heart-rate monitors and accelerometers during the sampling days, which recorded their heart rate and physical activity continuously for 24 h (Figures S1B–S1E). As expected, the subjects were more active and exhibited higher heart rate during the light compared with the dark phase (Figure S1F). The activity recorded during the dark phase is mostly attributed to the blood sampling for which the subjects were required to wake up. The participants’ activity intensity was mostly light, and a major portion of it was classified as sedentary even during the light phase. In addition, the participants were slightly more active in the high altitudes compared with SL (Figures 1D and S1C–S1F), probably because they were engaged in the expedition’s tasks, which were not required at SL. We observed an increase in heart rate with altitude (Figure 1D). Heart-rate differences during the light phase might stem from differences in locomotor activity levels. However, significant differences in heart rate strictly during the rest phase were observed as well (Figure 1E). These differences were likely due to the effect of altitude, as activity levels in this time frame were comparable (Figure 1E).

Overall, our physiological measurements confirm the hypoxemic state of the participants at higher altitudes. Notably, some parameters showed graded dependence on the altitude (e.g., saturation and heart rate at rest), while others did not change significantly between PU and RI (e.g., hemoglobin).

High altitude affects whole-blood gene expression

To study the effects of high altitude and time of day on whole-blood transcriptome, we performed RNA sequencing of whole-blood samples obtained at different times throughout the day in each of the altitudes, as described above. Principal-component analysis revealed that altitude explains a major part of the variability between the whole-blood transcriptome samples (Figure S2A). Accordingly, we initially focused on genes that were significantly affected by altitude (i.e., showed a significant overall change between SL, PU, and RI). Out of 12,726 genes detected, 3,775 genes (∼30%) were significantly affected by different altitudes with varied magnitudes of response (likelihood-ratio test [LRT] test, q < 0.01, |log2fold change(FC)| > 0.5; Figures S2B–S2D; see Table S2). We used k-means unsupervised clustering, as detailed in the STAR Methods, to group these genes based on their expression profiles into 6 clusters (Figures 2A and 2B). Clusters A1, A2, and A3 contain genes that were up-regulated upon altitude elevation, while cluster A6 contains down-regulated genes. Unexpectedly, clusters A4 and A5 show a more complex behavior, in which the response to altitude is non-monotonic: in cluster A4, there is up-regulation in PU followed by down-regulation in RI, and in cluster A5, there is down-regulation in PU followed by up-regulation in RI. To identify related biological functions for these different expression profiles, we performed pathway enrichment analysis using the GO Biological Processes (BP) database and upstream transcriptional regulators enrichment analysis using the ChEA 2016 (ChEA) database (Figures 3 and S3; Table S3).

Figure 2

High altitude affects whole-blood gene expression

(A) Heatmap representation of all genes significantly affected by altitude (Q < 0.01, |log2FC| > 0.5, LRT, n = 8 subjects), clustered according to their different expression profile (cluster serial numbering appears on the left). The data presented as a subject-wise Z score. (SL, Sea Level, 200m; PU, Puno, 3800m; RI, La Rinconada, 5100m).

(B) Boxplot representation of the mean expression per condition of all genes in each cluster from (A) and across all subjects, (Q < 0.01, |log2FC| > 0.5, LRT, n = 8 subjects). In boxplots; center indicates the median, the bottom and top edges indicate the 25th and 75th percentiles, each cluster size is indicated in parenthesis.

Figure 3

Functional signature of different expression clusters in response to high altitude

Pathway enrichment analysis on clusters of significantly altitude-affected genes as in Figure 2 (Q < 0.01, |log2FC| > 0.5, LRT), based on GO Biological Processes (BP) database. Terms were filtered and sorted based on adjusted p-value (Adj. P) and the number of hits from this term (N).

(A) Enrichment analysis for the monotonically upregulated cluster A3 (Adj. p < 0.01, N>=10).

(B) Enrichment analysis for the monotonic down-regulated cluster A6 (Adj. p < 0.01, N>=10).

(C) Enrichment analysis for the non-monotonic clusters A4-A5 (Adj. p < 0.05, N ≥ 5).

We first examined the clusters that showed monotonic behavior, namely clusters A1–A3 and A6. The up-regulated cluster A3 was mainly enriched for protein-biosynthesis-related processes (Figure 3A). The other up-regulated clusters (A1 and A2) showed no significant enrichment within the selected cutoffs. Subsequent analysis on the up-regulated clusters based on ChEA revealed enrichment for target genes of transcription factors involved in cell proliferation and RNA processing such as TTF2, FOXM1, E2F1, and MYC (Figure S3). On the other hand, the down-regulated genes (cluster A6) were enriched for innate immunity related terms, covering neutrophil activation and degranulation and ERK signaling (Figure 3B). Consistently, it was also enriched for targets of the TAL1, TCF7, and ELK3, which play a role in hematopoietic cell development (Figure S3).

To further corroborate these findings, we analyzed the distribution of innate-immunity-related genes and ribosomal proteins across the clusters and found a distinct behavior (Figure 4A). While genes encoding ribosomal proteins were mostly abundant in the up-regulated clusters, the innate-immunity-related genes were exclusively present in the down-regulated cluster (Figure 4A). Along this line, we found that the expression level of several ribosomal proteins (RPLs) is elevated with altitude (Figure 4B), whereas the expression of genes implicated in innate immunity (e.g., calcium-dependent nucleotidase [CANT1], the cytidine deaminase [CDA], and complement factor D [CFD], also known as adipsin) is down-regulated.

Figure 4

Distinct response of genes related to innate immunity and protein production to altitude

(A) Histogram showing the distribution across clusters of genes encoding ribosomal proteins or involved in innate immunity (GO:0002446).

(B) Representative ribosomal proteins encoding genes with significant altitude response (Q < 0.01, |log2FC| > 0.5, LRT, n = 8 subjects).

(C) Representative genes involved in innate immunity with significant altitude response (Q < 0.01, |log2FC| > 0.5, LRT, n = 8 subjects) (SL, 200 m; PU, 3,800 m; RI, 5,100 m). In boxplots, center indicates the median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively.

HIF1A is known as a master regulator of cellular response to hypoxia. We therefore examined HIF1Α transcript levels in whole-blood samples at different altitudes (Figure S4) and observed a decline in HIF1Α transcript levels with high altitude. While HIF-1A protein is transiently up-regulated in hypoxia, mRNA is frequently found to be repressed in cell culture assays (Bruning et al., 2011; Cavadas et al., 2015; Chamboredon et al., 2011; Uchida et al., 2004). We next examined whether HIF1A targets are enriched among our altitude-responsive genes. Our BP and ChEA analyses did not show enrichment for the HIF pathway. Furthermore, we specifically examined changes in expression levels of several well-established target genes of HIF1A pathway (Figure S4) and did not observe a coherent response that is consistent with activation of HIF1A pathway in an altitude-dependent manner. This is in line with the reduced involvement of HIF1A in chronic hypoxia compared with acute exposure (Saxena and Jolly, 2019). Taken together, based on our analyses and previous evidence from animal models, it is likely that the observed transcriptional response to high altitude is mostly HIF1A independent.

Collectively, we concluded that altitude induces the expression of genes involved in protein translation and processing. This may suggest higher demand for protein production (e.g., to increase hemoglobin levels) to support adaptation to high altitude and low oxygen availability. On the other hand, we observed down-regulation of genes implicated in innate immunity such as myeloblast differentiation and neutrophil activation, which might lead to alterations in innate immunity functions. It is possible that these two processes, namely the protein synthesis and innate immunity response, act together to shape the blood immune cell composition upon high altitude, as suggested below.

Different altitudes above SL differentially affect gene expression

Almost 1,500 vertical meters separate PU and RI, generating a further 15 mmHg drop in oxygen partial pressure, and, as shown above, affect both physiological (Figure 1) and molecular parameters (Figure 2). Interestingly though, clusters A4 and A5 exhibited a response to altitude that differed in its directionality between PU and RI (Figure 2). To examine the potential functional significance of these differentially expressed genes between PU and RI, we performed an enrichment analysis focusing exclusively on clusters A4 and A5 (Figure 3C). Cluster A4 includes genes that were up-regulated in PU compared with SL but down-regulated in RI compared with PU. This cluster was enriched for interferon type 1-related genes. Cluster A5, the mirror image of cluster A4, contained genes that were down-regulated in PU but up-regulated in RI. This cluster was enriched with terms related to megakaryocyte, myeloid cell, and erythrocyte differentiation, as well as platelet degranulation (Figure 3C). In addition, these clusters showed enrichment for targets of TAL1/SCL, EKLF, and GATA1/2, transcription regulators involved in hematopoietic processes (Figure S3). Hence, although altitude monotonically down-regulated genes involved in myeloblast differentiation and neutrophil activation (cluster A6), genes related to differentiation of erythrocytes and megakaryocytes exhibited a non-monotonic effect and were down-regulated in PU and elevated in RI. Here, again, these specific changes in expression of hematopoietic genes are likely to shape the immune cell composition or response upon changes in altitude.

We concluded that high altitude differently affects the transcriptional response in PU compared with RI and does not necessarily unidirectionally follow the drop in oxygen availability. Nevertheless, it is important to note that due to limitations in the study design as discussed below, we cannot distinguish between effects that stem from differences in altitude versus differences in duration at high altitude.

Whole-blood gene expression shows altitude-dependent daily variance

In recent years, a connection between hypoxia and circadian rhythms was elucidated, mainly in animal models (Adamovich et al., 2017, 2019, 2022; Liu et al., 2022; Manella et al., 2020) and to some extent in humans (Bosco et al., 2003). This prompted us to ask whether daily variance in gene expression changes between the different altitudes, characterized by varying oxygen availability. As described above, we collected the blood samples at 4-h intervals over 24 h, which enabled assessment of the time-of-day effect in addition to the overall effect of altitude. We first analyzed the temporal transcriptomic data using several well-established algorithms for analyzing circadian rhythms, namely JTK, Lomb-Scargle periodogram, ARSER (Wu et al., 2016), and harmonic regression (Luck et al., 2014). These algorithms were designed to identify periodicity of a specific frequency, with an overall profile similar to a cosinor wave. However, analyzing the data based on this approach resulted in non-significant results across all algorithms we have tried (false discovery rate [FDR] <0.05; Table S4). It is likely that given the nature of this field study, the rhythms previously described in whole-blood RNA (Archer et al., 2014; Moller-Levet et al., 2013) were disrupted, for example, by the night awakenings for blood sampling.

In an attempt to obtain insights regarding potential time-of-day effects, we took a different approach and explored the differences between the active and rest phases. To this aim, we separately pooled the samples that were obtained in the light (10 a.m., 2 p.m., and 6 p.m.) or dark (10 p.m., 2 a.m., and 6 a.m.) phases, which corresponded to the activity and rest phases, respectively, in line with our activity monitoring profiles (Figure S1). We presumed that this procedure would result in attenuation of the experimental noise and provide higher statistical power.

Overall, 807 genes (∼20% of the detected genes) showed significant difference with varied magnitudes between the light and dark phases in any of the altitudes (LRT test, q < 0.05) (Figures 5A and S5A–S5C), out of which ∼30% overlapped with previously reported genes that exhibited diurnal expression pattern (Figure S5D) (Moller-Levet et al., 2013) and ∼40% overlapped with genes that were significantly affected by altitude (Figure S5E; Table S2). The significant genes were grouped based on unsupervised clustering into 6 distinct clusters (Figures 5A and 5B). Clusters T1–T4 are characterized by higher expression levels in the light versus the dark phase, while clusters T5 and T6 have higher expression in the dark phase. We also observed differences in the effect of altitude on the daily variance: in clusters T1 and T5, the daily variance was overall unaffected by altitude, while in the other clusters, an interaction with altitude was apparent (one-way ANOVA and Tuckey’s post-hoc test, |ΔZ score| ≥ 0.2, p < 0.01) (Figure 5B). The altitude-dependent clusters T4 and T6 exhibited increased and decrease magnitude of daily variance with altitude, respectively. Notably, the T2 and T3 clusters presented a non-monotonous interaction with altitude. Overall, the majority of the time-affected clusters were enriched for the terms related to innate immunity and, in particular, to neutrophil activation, as well as targets of the transcription factors RELA and IRF8, mediators of inflammation and immune responses (Figures S5F and S6). Representative genes from these clusters include and IL1B (cluster T1) and CCR7 (cluster T5), as well as NR1H2 (cluster T6) (Figure 5C). In cluster T2, the daily variance in gene expression was less pronounced in PU compared with in RI and SL. This cluster was enriched for a distinct subset of terms, many of which related to megakaryocyte, erythrocyte, and myeloid differentiation (Figure S6B), consistent with the expression of GATA1 (Figure 5C), and enrichment for GATA1, SCL/TAL1, and KLF1 targets (Figure S5F). Cluster T2 highly resembled the non-monotonic altitude-affected cluster A5 (Figure 2), as was also apparent from the enrichment results (Figures 3C and S3).

Figure 5

Daily variance in whole-blood transcriptome is altitude dependent

(A) Heatmap representation of all genes that showed significant daily variance between the light and dark phases (in any of the altitudes; Q < 0.05, LRT, n = 8 subjects), clustered according to their different expression profiles (cluster serial numbering appears on the left). The data presented as subject-wise Z score. (SL, 200 m; PU, 3,800 m; RI, 5,100 m).

(B) Boxplot representation of the difference between the light- and dark-phase expression levels per altitude of all genes in each cluster from (A) averaged across all subjects (cluster size in parenthesis).

(C) Representative gene profiles for each cluster in (A). In boxplots, center indicates the median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively.

Next, we analyzed the response of clock genes to altitude and time of day (Figure 6). The circadian repressor and immunomodulatory factor NR1D1 showed a significant response to both altitude and time of day. Its expression was elevated in the dark compared with in the light phase, and overall, it was down-regulated with altitude. The expression profile of DBP resembled that of NR1D1 and was altitude dependent. ARNTL, a target of NR1D1 repression, exhibited an opposite phase compared with NR1D1, namely its expression levels were elevated in the light compared with in the dark phase. However, it was non-monotonically affected by altitude. RORA, an activator of ARNTL expression, was up-regulated by altitude with no time-of-day effect. Notably, changes in expression level of several clock genes followed a trend in response to altitude; CRY2 and PER1 were down-regulated, while NR1D2, PER3, and RORC were up-regulated, upon altitude.

Figure 6

Time and altitude effects on clock genes

Boxplot representation of clock genes’ expression across different altitudes and time of day. (n = 8 subjects, significant altitude effect marked as A, significant time effect as T, according to LRT as in Figures 2 and 5, respectively). In boxplots, center indicates the median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively.

Altogether, our analysis identified daily variance in the whole-blood transcriptome, which was altitude- and likely oxygen-availability dependent. As in the case of altitude effect, some of the altitude-dependent time effects also showed clear non-monotonic behaviors. All clusters contained genes largely related to functions of different cell types of the myeloid lineage. However, while the non-monotonic altitude-dependent daily variance was associated with megakaryocyte and erythrocyte differentiation, genes related to neutrophil activation showed diurnal variations in expression that were altitude independent. Furthermore, we detected both altitude and time-of-day effects on the expression of several clock genes, which might suggest an effect of altitude and low oxygen availability on the clock’s function.

Digital cytometry reveals altitude- and time-dependent variance in blood cell-type composition

Hitherto, our analysis revealed that many of the altitude- and time-dependent differences in gene expression are immune related. Whole-blood samples contain a mixture of transcripts from multiple cell types. Hence, it is conceivable that changes in the transcriptome are partially due to differences in the cytological composition of the blood. This is especially true for the immune response, which involves cell infiltration, proliferation, and death of specific cell types. To gain insights regarding the differences in cytological composition, we employed CIBERSORTx, a digital cytometry approach that assesses the cell-type proportions in mixed transcriptomes using expression signatures unique to each cell type in the sample (Newman et al., 2019). As a signature dataset, we used the LM22, which contains signatures of 22 different leukocyte cell types and was previously used for assessing cell-type proportions in whole-blood samples (Newman et al., 2015).

We found that, indeed, several cell types showed a significant difference in abundance following altitude changes. For instance, monocyte and M2 macrophage proportions increased, while resting natural killer (NK) cells and regulatory T cells decreased with altitude (Figure 7). Moreover, some cell types, predominantly neutrophils, showed a daily difference in their abundance, in accord with previous reports and in line with our enrichment analyses (Figure S6) (Aroca-Crevillen et al., 2020; Crespo et al., 2020; Ella et al., 2016). Interestingly, we also found cell types that had an altitude-dependent daily variance, e.g., monocytes had an increased abundance in the dark phase, and this daily difference increased with altitude.

Figure 7

Altitude and time effects on blood cellular composition based on digital cytometry

Boxplot representation of digital cytometry analysis on the RNA sequencing (RNA-seq) data using CIBERSORTx with the LM22 signature matrix. Each panel represents the relative abundance of a single cell type across the experimental conditions. ∗p < 0.05, ∗∗p < 0.01, ns, non-significant, two-sided paired t test; #p < 0.05, ##p < 0.01, ###p < 0.001, repeated-measures two-way ANOVA, for the denoted effect (time, time of day effect; Alt, altitude effect, or interaction); n = 8 subjects. (SL, 200 m; PU, 3,800 m; RI, 5,100 m). In boxplots, center indicates the median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively.

These results suggested that altitude alters the blood cytological composition and, in some cases, affects the daily variance in cell-type abundance. Together with the results above, our findings highlight the intricate effects of high altitude and low oxygen availability on the immune system.

Discussion

The extreme living conditions at high altitude confer multi-level responses and adaptations, including myriad behavioral, physiological, and molecular changes (Grocott et al., 2007). The physiological and molecular response mainly depends on the time spent in the high altitude (i.e., acute versus chronic exposure) and on the intensity of the environmental changes (i.e., the altitude above SL) (Grocott et al., 2007). The blood plays a pivotal role in this response, primarily given its principal role as an oxygen carrier but also due to its other roles (e.g., coagulation, immune, and endocrine signaling). Few studies in the past investigated the blood transcriptional response to high altitude in humans (Chen et al., 2012; Gaur et al., 2020; Liu et al., 2017). The current study has several important advantages compared with previous works. First, it examines the effect of high altitude on lowlanders specifically following an initial acclimatization step (mid-term exposure). Second, it includes two different high-altitude levels (3,800 and 5,100 m), where the latter is considered a “very high altitude” and thus enables comparison between two distinct high-altitude regimens. Third, our design includes sample collection from multiple time points around the clock, which addresses the question of time-of-day effects.

Our analysis revealed two distinct responses to high altitude, namely monotonic and non-monotonic responses. Most of the genes belonged to clusters that were monotonically affected by altitude. Gene enrichment among these clusters spotted genes related to innate immunity, most notably neutrophil functions, as well as protein biosynthesis. Innate-immunity-related genes and genes involved in protein biosynthesis showed an opposite response to altitude: the former were down-regulated, while the latter were up-regulated, with altitude. The less-trivial and counter-intuitive non-monotonic response consisted of genes that respond to the change between PU and RI in the opposite direction compared with that between SL and PU. Interestingly, clusters that behaved in this manner were enriched for example with genes related to erythrocyte, myeloid cell, and megakaryocyte differentiation. We speculate that this effect can be partly seen as a decreased activation due to sufficient pre-acclimatization and adaptation in PU, which is followed by re-activation at the higher altitude of RI. It is also possible that above a certain altitude and hypoxia severity, between PU and RI, a molecular turn point occurs, and the response to altitude changes. Together, these results highlight the complex molecular adaptation to altitude, which is likely affected by different environmental variables (e.g., oxygen, temperature) to regulate different biological processes at different rates.

The second unique behavior was of genes that presented an interaction between their time-of-day effect and their altitude effect. These genes showed a daily variance in their expression levels, but this variance was dependent on the altitude. Intriguingly, there was a major overlap between the non-monotonic altitude clusters (clusters A4 and A5) and the temporal non-monotonic interaction clusters (clusters T2 and T3). This was also reflected in the enrichment analysis, which draws attention toward the branch of the hematopoietic lineage that gives rise to erythrocytes and megakaryocytes (which are the progenitors of the platelets). The effect of high altitude on erythropoiesis is well known (Haase, 2013). In addition, previous studies have shown that high altitude increases the rate of blood coagulation (Hancco et al., 2020) and affects platelets transcriptome and proteome (Shang et al., 2020). Our findings suggest that these effects of high altitude and low oxygen availability might be time- and altitude dependent.

The effect of high altitude on the immune system was not thoroughly investigated in humans (Mishra and Ganju, 2010), and the studies conducted thus far raise a complex and not always coherent picture. While some studies indicate immune suppression at high altitude, others point toward increase in the pro-inflammatory response (Pham et al., 2021; Wang et al., 2021). Some of the discrepancies might stem from differences in experimental protocols and might relate to time-of-day effect as we report herein. Sustained hypoxia, a major factor in the pathophysiology of mountain sickness, was shown to modulate the immune response in animal models, including exacerbation of inflammatory responses. The crosstalk between hypoxia and inflammation is implicated in various morbidities such as chronic obstructive pulmonary disease and obstructive sleep apnea, as well as cardiovascular diseases (Kammerer et al., 2020; Liu et al., 2017; Mishra and Ganju, 2010). Our digital cytometry analysis shows an altitude- and time-dependent increase in M2 macrophages and monocytes but an altitude-dependent decrease in NK and regulatory T cells. Unfortunately, due to the lack of appropriate lab equipment, we collected and froze only whole-blood samples, which are inappropriate for further cytological or molecular analyses.

It is well established that many elements of the immune system are regulated in a daily or circadian manner (Scheiermann et al., 2013). However, the combined effect of time and altitude on immunity remains understudied. Interestingly, we found that the expression of some clock genes is affected either by altitude or time of day. NR1D1 showed a day-time difference in expression, and its overall expression was significantly down-regulated by altitude. It is noteworthy that NR1D1 has been implicated in regulation of the immune response (Gibbs et al., 2012; Poole and Kitchen, 2022). NR1D1−/− mice loss their circadian gating of endotoxin response, and pharmacologic or genetic modulation of NR1D1 alters the production and release of pro-inflammatory cytokines such as interleukin-6 (IL-6). Hence, it is tempting to speculate that NR1D1 might play a role as the molecular link between the circadian clock and the immune function in response to low oxygen availability.

Our gene enrichment analysis as well as the digital cytometry analysis indicate that indeed both altitude and time of day affect immune cells, and moreover, on some occasions, the two interact with each other. It is thus possible that some processes of altitude adaptation are optimized by daily rhythmicity.

In summary, our data provide a unique perspective on the effect of high altitude and low oxygen availability on the blood transcriptome and can serve as a good basis for further investigation.

Limitations of the study

The current study includes some inherent limitations. First, the SL condition was measured in a different continent quite a long time (4 months) after the expedition ended. While this ensures that we measure the true baseline without residual long-lasting effects from exposure to high altitude, it can also introduce confounding effects. Along the same line, we cannot completely separate between the effect of cumulative time in high altitude and of the relative altitude. In other words, we cannot tell which differences between PU and RI are due to longer exposure and which are due to the difference in altitude. Another caveat was the use of acetazolamide by 6 out of the 8 participants, which might increase the inter-subject variability. This sub-optimal design is the result of the obvious difficulties of conducting such a field study: the study was annexed to an existing expedition, faced schedule and equipment constrains, and was conducted in harsh conditions. Additionally, it is important to mention that high altitude carries many additional effects beyond reduced oxygen availability. These include changes in barometric pressure, humidity, temperature, radiation, and others (Coppel et al., 2015; Faiss et al., 2013), and we cannot distinguish between the effects stemming from hypoxia from those stemming from hypoxia-independent parameters.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Gad Asher (gad.asher@weizmann.ac.il).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Participants

Eight healthy male lowlanders living near sea level (<400 m) participated in the present study (Table S1). All of them are non-smokers, free of any known health problems and chronic medications. All participants were scientists involved in a large research project focusing on high-altitude adaptation (Expedition 5300 research program) in Peru. All participants were fully informed about the study and signed a written consent before inclusion. The study was approved by the ethics committee of Universidad National Mayor de San Marcos (Lima, Peru, CIEI-2019-002) and was conducted in accordance with the Declaration of Helsinki.

Study protocol overview

To examine the around the clock transcriptomic changes induced by exposure to high altitude, blood samples were taken at 4-h intervals throughout 24 h (i.e., 6:00, 10:00, 14:00, 18:00, 22:00 and 02:00) from each participant, at 3 different altitudes, during and after the scientific expedition conducted in the Altiplano area in southeast Peru. Participants were first pre-acclimatized (3 days in 3,800 m, and 7 days in 5,100 m), and then spent an additional week at 3,800 m. The first high-altitude 24-h sampling round was conducted on the 7th day at 3,800 m (barometric pressure 483 mmHg, Puno, Peru). Afterwards, the participants acclimatized for one week at an altitude of 5,100 m (barometric pressure 418 mmHg, La Rinconada, Peru), and a second high-altitude 24-h sampling round was conducted on the 7th day. The 24-h control sea level sampling round was conducted 4 months after returning from high altitude. Travelling to high altitude was performed by airplane and by car. During the three sampling rounds, participants had a standardized life style and only performed their scientific research activities. No coffee, alcohol and stimulants were taken. Lifestyle conditions were controlled and standardized for all participants, especially meal times (breakfast between 7:30 and 8:00, lunch between 13:00 and 13:30, dinner between 20:00 and 20:30), bedtime (22:30) and wake-up time (5:30). No medications were taken by the participants during the study periods, except acetazolamide for 6 of the participants to prevent acute mountain sickness at high altitude, which was similarly administrated at all three altitudes (125mg, twice daily at 8:00 and 12:00, the day before and during the 24-h sampling day) (Luks et al., 2019). None of the participants suffered from acute mountain sickness (as defined by the 2018 Lake Louis Score (Roach et al., 2018)) during the study period. Sea level measurements of hemoglobin concentration, hemoglobin mass and hematocrit are lacking for one of the participants due to technical reasons.

Method details

Accelerometer and heart rate data collection and analysis

During each 24-h sampling round, physical activity and heart rate (HR) were monitored using a wrist-worn triaxial accelerometer (Actigraph GT9X, ActiGraph Corp, Pensacola, FL, USA) associated with a Bluetooth HR transmitter chest belt (Polar H10 Heart Rate Monitor, Polar Electro Oy, Kempele, Finland). From activity counts registered per minute, data were aggregated per hour to determine fraction of time spent in sedentary, light or moderate-to-vigorous activity.

Hemoglobin and gas measurements

Measurements of blood parameters (hematocrit, hemoglobin, oxygen saturation) and gas exchange were performed during the sampling round in each altitude, between 10:00 and 14:00.

Hematocrit was measured in duplicate from a venous blood sample, using microcentrifuge method (HemataSTAT II, Separation Technology Inc., Sandford, USA). Blood hemoglobin concentration ([Hb]) was measured using an ABL80 blood gas analyzer (Radiometer, Copenhagen, Denmark).

Hypoxemia levels induced by high-altitude exposure were assessed at each altitude location by pulsed oximetry (SpO2), done after a 5 min sitting rest period through a finger sensor (NELLCOR OxiMaxN-65, TycoHealthcare, CA). In addition, at each altitude location, a daily measurement of end-tidal partial oxygen (PetO2) and carbon dioxide (PetCO2) pressures were also performed (AGM100, MediPines, Yorba Linda, CA, USA), as a surrogate for arterial pressure values.

Total hemoglobin mass was determined based on the carbon monoxide (CO) rebreathing method, using an automated system (OpCO, Detalo Health, Copenhagen, Denmark) with a dose of 1.2 mL CO·kg−1, as previously described (Oberholzer et al., 2020). Just prior and at 10 min of CO rebreathing, a venous sample was obtained to measure the hematocrit (microcentrifuge method), the percentage of carboxyhemoglobin and the hemoglobin concentration (ABL80, Radiometer, Copenhagen, Denmark).

RNA extraction

At each altitude level (SL, PU, and RI) and time around the clock (06:00, 10:00, 14:00, 18:00, 22:00 and 02:00), a whole blood sample (2.5 mL) was collected in a dedicated blood RNA tube (PAXgene, PreAnalytiX, Qiagen, Hombrechtikon, Swizerland) from a venipuncture of the arm for each participant (144 samples in total). Blood RNA tubes were immediately frozen and stored at −20°C after sampling for further analysis.

RNA from the samples collected on the PAXgene tubes was isolated according to the manufacturer protocol (PAXgene Blood RNA kit #762174, PreAnalytiX, Qiagen, Hombrechtikon, Swizerland). Briefly, PAXgene tubes were thawed at 20°C for 2 h 30 min. Then, RNA was isolated, eluted in water (ultrapure DEPC-treated water #750024, Thermo Fisher Scientific), flash-frozen in liquid nitrogen and stored at −80°C. RNA concentration and purity was measured using NanoDrop ND100 spectrophotometer (Thermo Fisher Scientific).

MARS-seq library preparation, sequencing, and processing

We used a derivation of MARS-seq as described (Jaitin et al., 2014), originally developed for single-cell RNA-seq to produce expression libraries, and exclusively sequencing the 3′-end of the transcripts. The prepared bulk MARS-seq libraries were sequenced on NovaSeq 6000 (Illumina), using the NovaSeq 6000 SP Reagent Kit v1.5 (100cycles) (Illumina). The entire library was run twice, in order to increase depth. The fastq files from both runs were concatenated for downstream processing.

Processing of raw sequencing data into read counts was performed via the User-friendly Transcriptome Analysis Pipeline (UTAP) (Kohen et al., 2019). In short, Reads were trimmed using cutadapt (Martin, 2011) and mapped to genome (/shareDB/iGenomes/Mus_musculus/UCSC/mm10/Sequence/STAR_index) using STAR (Dobin et al., 2013) (default parameters). The pipeline quantifies the genes annotated in RefSeq (that have expanded with 1000 bases toward 5′ edge and 100 bases toward 3′ bases). Counting (UMI counts) was done using HTSeq-count in union mode (Anders et al., 2015). Normalization of the counts was performed using DESeq2 (version 1.34.0) (Love et al., 2014) with default parameters.

The raw fastq files and raw counts summary table are available through the GEO database (accession number GSE196728).

RNA-seq statistical analysis

Principal component analysis was performed on subject-wise normalized values, with “stats” package, version 4.1.2 in R (Core R Team).

Differential gene expression was evaluated by DESeq2 (Love et al., 2014), using the Likelihood Ratio Test (LRT) test with two designs. Both designs were based on the full model Y ∼ Subject + Altitude + Time + Altitude:Time, which is a repeated-measure model due to the Subject term. The first design tested for any effect for altitude by using the reduced model Y ∼ Subject + Time. The second design tested for any time effect within each altitude, so the reduced model was Y ∼ Subject + Altitude. In both cases, the Time variable had two levels: Light (time points 10, 14, 18) and Dark (time points 22, 2, 6). The normalized data together with the statistical tests results are available in Table S2.

Clustering and gene enrichment testing

Clustering was performed on subject-wise Z-score normalized data. Significant genes in each comparison (see cutoffs in the main text and figure legends) where clustered according to their expression profiles, using the k-means method with the with Hartigan and Wong algorithm (Hartigan and Wong, 1979) as implemented in “stats” package, version 4.1.2 in R. We used 6 centers with 50 random starts; default settings were used for the rest of the parameter. The number of centers was selected to balance between the clusters’ size and the separation quality.

Gene enrichment in each cluster was tested using the “EnrichR” package, version 3.0 (Kuleshov et al., 2016), using GO (Ashburner et al., 2000; The Gene Ontology, 2019) and ChEA2016 (Lachmann et al., 2010) databases. Term size and p value cutoffs are stated in the figure legends. The full enrichment tests results can be found in Table S3.

Rhythmicity analysis

To assess rhythmicity in activity and heart rate measurements, we used JTK_CYCLE in repeated measures design as implemented in MetaCycle R package, version 1.2.0 (Wu et al., 2016).

To assess rhythmicity in gene expression, we used several algorithms: JTK_CYCLE, Lomb-Scargle Periodogram, ARSER (all from MetaCycle R package, version 1.2.0 (Wu et al., 2016)), and Harmonic Regression (HarmonicRegression R package, version 1.2.0 (Luck et al., 2014)). We used them both in repeated measures mode and in bulk mode. In all cases, no more than two genes passed FDR <0.05. The results are available in Table S4.

Digital cytometry

The CIBERSORTx online tool (Newman et al., 2019) was used to assess the proportion of different cell types in the blood samples based on expression signatures. For the expression signatures we utilized the “LM22” datasets which was generated by the CIBERSORTx authors (Newman et al., 2019) and includes 22 leukocyte cell types. Settings: run in absolute mode, number of permutations set to 500, verbose output was enabled, B-mode batch correction was applied. The full CIBERSORTx output is provided in Table S5.

Quantification and statistical analysis

All the statistical analyses and data visualization were performed with R (v4.1.2). Data visualization with R was performed using the ggplot2 (Wickham, 2016) and ComplexHeatmap (Gu et al., 2016) packages. In all boxplots: middle line = median, box = 25th to 75th Inter Quantile Range (IQR), whiskers = the largest/smallest value no greater/smaller than 1.5∗IQR, outlier points = measurements outside this range. All the statistical details of experiments (the statistical tests used, the value of n, and what it represents alongside the definition of center, dispersion and precision measures) are provided in the figure legends.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Human blood samples	This study	N/A
Critical commercial assays
PAXgene Blood RNA kit	PreAnalytiX, Qiagen, Hombrechtikon, Switzerland	#762174
Ultrapure DEPC-treated water	Thermo Fisher Scientific	#750024
Deposited data
RNA-seq data	This study	GEO accession number: GSE196728
CIBERSORTx data	This study	Table S5
Individual Subject Data	This study	Table S1
Software and algorithms
R (v4.1.2)	CRAN	N/A
CIBERSORTx	Newman et al., 2019	https://cibersortx.stanford.edu/
Other
Micro-centrifuge method: HemataSTAT II	Separation Technology Inc., Sandford, USA	#100-100
ABL80 blood gas analyzer	Radiometer, Copenhagen, Denmark	MAPSSS-000120
AGM100	MediPines, Yorba Linda, CA, USA	RS-0004-21
OpCO	Detalo Health, Copenhagen, Denmark	https://detalo-health.com/
Wrist-worn triaxial accelerometer: Actigraph GT9X	ActiGraph Corp, Pensacola, FL, USA	GT9X
Bluetooth HR transmitter chest belt: Polar H10 Heart Rate Monitor	Polar Electro Oy, Kempele, Finland	Polar H10
Oxygen finger sensor: NELLCOR OxiMaxN-65	Medtronic, Dublin, Ireland	OxiMaxN-65