Genome-Wide DNA Methylation Changes Associated With High-Altitude Acclimatization During an Everest Base Camp Trek

The individual physiological response to high-altitude hypoxia involves both genetic and non-genetic factors, including epigenetic modifications. Epigenetic changes in hypoxia factor pathway (HIF) genes are associated with high-altitude acclimatization. However, genome-wide epigenetic changes that are associated with short-term hypoxia exposure remain largely unknown. We collected a series of DNA samples from 15 participants of European ancestry trekking to Everest Base Camp to identify DNA methylation changes associated with incremental altitude ascent. We determined genome-wide DNA methylation levels using the Illumina MethylationEPIC chip comparing two altitudes: baseline 1,400 m (day 0) and elevation 4,240 m (day 7). The results of our epigenome-wide association study revealed 2,873 significant differentially methylated positions (DMPs) and 361 significant differentially methylated regions (DMRs), including significant positions and regions in hypoxia inducible factor (HIF) and the renin–angiotensin system (RAS) pathways. Our pathway enrichment analysis identified 95 significant pathways including regulation of glycolytic process (GO:0006110), regulation of hematopoietic stem cell differentiation (GO:1902036), and regulation of angiogenesis (GO:0045765). Lastly, we identified an association between the ACE gene insertion/deletion (I/D) polymorphism and oxygen saturation, as well as average ACE methylation. These findings shed light on the genes and pathways experiencing the most epigenetic change associated with short-term exposure to hypoxia.


INTRODUCTION
Altitude acclimatization in humans is characterized by complex physiological responses, which include the cardiovascular, hemopoietic, respiratory, and metabolic systems [for review, see Palmer (2010)]. Each system responds uniquely to low oxygen environments. For example, cardiovascular output increases (i.e., increased heart rate and stroke volume) upon initial altitude exposure and returns to pre-altitude baseline after several days of acclimatization [for review, see Naeije (2010)]. The respiratory system's response is to initiate hyperventilation. The hypoxic ventilatory response (HVR) is elicited shortly upon exposure to high altitude, with ventilatory acclimatization emerging following 5-7 days of sustained exposure to hypoxia (Powell et al., 1998). Lastly, the hemopoietic response in the form of increased erythrocyte production is evident after several days to weeks of exposure (Rodriguez et al., 2000). Each of these responses facilitates acute acclimatization to the low ambient oxygen tension present at high altitudes, allowing humans to acclimatize to hypoxic conditions. Epigenetic change is one mechanism through which physiological acclimatization may occur. Epigenetic modifications can affect gene expression and include DNA methylation, histone tail modifications, and short RNA regulation. The most well-studied epigenetic mark is DNA methylation, the addition of a methyl group primarily to cytosine bases. DNA methylation patterns can change upon exposure to various environmental conditions, including exposure to different diets, stress, and toxicants (Dolinoy et al., 2007;Baccarelli et al., 2009;Colacino et al., 2012;Childebayeva et al., 2019b). Previous studies have demonstrated that changes in DNA methylation are associated with exposure to the low oxygen environment of high altitude (Alkorta-Aranburu et al., 2012;Childebayeva et al., 2019aChildebayeva et al., ,b, 2020. These studies show that genes in the hypoxia inducible factor (HIF) pathway exhibit changes in DNA methylation associated with high-altitude exposure.
The HIF pathway is the main oxygen sensing pathway that regulates cellular homeostasis in metazoans (Bigham and Lee, 2014). The pathway takes its name after the master transcriptional regulator HIF, a heterodimeric transcription factor that is formed by one of three α-subunits (HIF-1α, HIF-2α, or HIF-3α) and a β-subunit (also known as ARNT). In normoxia, HIF1A is hydroxylated and subsequently degraded by the ubiquitin-proteosome pathway. Under hypoxia, this hydroxylation is inhibited by the lack of oxygen availability, leading to the dimerization of HIF and activation of target genes. HIF is responsible for transducing changes in oxygen tension to changes in gene expression through hypoxia response elements (HREs) (Wang and Semenza, 1995;Kaelin and Ratcliffe, 2008;Semenza, 2012). The renin-angiotensin system or RAS is a second pathway that is involved in the response to hypoxia. It is one of the body's most important regulators of blood pressure and inflammation (Muller et al., 1997;Rupert et al., 2003). The RAS protein, angiotensin converting enzyme (ACE), is a central peptide in bloodpressure regulation responsible for converting angiotensin-I to the vasoconstrictor, angiotensin-II. An insertion/deletion (I/D) polymorphism in ACE is associated with physical performance at high altitude (Woods et al., 2002;Tsianos et al., 2005). The I-allele has been associated with higher levels of submaximal oxygen saturation (SaO 2 ) among Andean Quechua (Bigham et al., 2008), and in trekkers of European ancestry (Woods et al., 2002).
Previous research by our group has shown that acclimatization to hypoxia is associated with DNA methylation changes in HIF pathway genes including EPAS1, EPO, PPARa, and RXRA (Childebayeva et al., 2019a). However, it is not well understood what other genes and pathways display DNA methylation changes upon exposure to hypoxia. To understand how acclimatization to hypoxia affects genomewide DNA methylation patterns, we performed an epigenomewide association study in individuals trekking to Everest Base Camp. Our analysis compared baseline methylation measured in Kathmandu, Nepal at 1,400 m (day 0) with methylation measured at a high-altitude location, Pheriche, Nepal at 4,240 m (day 7 of the trek).

Ethics Statement
Ethical approval was received from the Syracuse University Institutional Review Board (Protocol 18-006) and the University of Michigan Institutional Review Board (HUM00141118). The study abided by the Canadian Government Tri-Council policy on research ethics with human participants (T2) and the Declaration of Helsinki, except for registration in a database. Ethical approval was received also from the Mount Royal University Human Research Ethics Board (Protocol 100012 and 101361) and harmonized with the Nepal Health Research Council (Protocol 109-2017).

Study Design and Sample Collection
Thirty-two samples (16 samples at 1,400 m and 16 samples at 4,240 m) corresponding to 16 unique individuals were selected from a larger participant cohort from the research expedition to Everest Base Camp in the Nepal Himalaya (Childebayeva et al., 2019a). Briefly, study participants and researchers flew from Kathmandu (baseline) to Lukla from where the research group trekked for 10 days from 2,800 to 5,160 m (Figure 1). In the morning between 06:00 and 08:00 local time at 1,400 m (Kathmandu; day 0) and 4,240 m (Pheriche; day 7), saliva samples for DNA and physiological measures were taken following one night of sleep at each altitude. Physiological measurements included hemoglobin concentration [Hb], patient end-tidal carbon dioxide (P ET CO 2 ), a measure of CO 2 partial pressure in expired air, which reflects the CO 2 level in the arterial blood, and peripheral oxygen saturation (SaO 2 ). Detailed information on phenotype collection and sampling is provided in Childebayeva et al. (2019a). All participants were healthy, non-pregnant, non-lactating, non-smokers between 19 and 41 years of age. All participants were of self-reported European ancestry and had at least 1 year since their last altitude experience. Participant characteristics can be found in Table 1.

Phenotype Testing
We performed linear mixed models using the R package lmerTest to test for significant differences in phenotypes between Kathmandu and Pheriche ( Table 1). The following model was tested: Phenotype ∼ Altitude + Sex + Age + (1 | ID).  Age, year 23.6 (6.0) Data are means (SD). Significance symbols denote the difference between Kathmandu baseline and each altitude. # p-value < 0.10; *p-value < 0.01; and **p-value < 0.001.

DNA Methylation
We generated DNA methylation data for ∼850,000 CpG sites using the Illumina Infinium R MethylationEPIC BeadChip assay for 32 samples in our study. We used the EZ-96 DNA Methylation TM Kit (Zymo Research, Irvine, CA, United States) to bisulfite convert each DNA sample following the standard protocol with alternative incubation conditions optimized for the Illumina Infinium R MethylationEPIC BeadChip assay. We used R for data processing and analysis implementing the packages minfi, ChAMP, and SmartSVA (Aryee et al., 2014;Morris et al., 2014;Chen et al., 2017). Based on QC metrics, two samples from the same participant failed and were excluded from all analyses; thus, the final sample size was n = 30 (15 at 1,400 m and 15 at 4,240 m). Data normalization was performed using the funnorm normalization function in minfi (Aryee et al., 2014). We removed all probes that were above the 10e5 detection p-value threshold (N = 8,126) in more than 5% of the samples, all crossreactive probes, probes associated with sex chromosomes, probes containing SNPs with MAF > 5% at target CpG sites, single base extension sites of type I probes, and in the body of the probe (Chen et al., 2013). All analyses were performed with N = 657,569 sites after normalization and probe removal. Samples were tested for batch effects using singular value decomposition (SVD) analysis in champ. SmartSVA (Chen et al., 2017) was used to perform a surrogate variable test, and the surrogate variable 1 was used for correcting for any saliva cell type differences associated with altitude. SmartSVA is a surrogate variable analysis method that can be used for reference-free adjustment for cell mixtures (Chen et al., 2017).

Differential Methylation Testing
Fully processed M-values were tested for differential methylation using the package limma (Ritchie et al., 2015). The following model was used to test for the differentially methylated positions (DMPs): DNA methylation ∼ Sample ID + Altitude + Surrogate Variable 1 (from smartSVA). P-values were adjusted for multiple testing using the false-discovery rate (FDR) following the Benjamini-Hochberg procedure (Hochberg and Benjamini, 1990) to produce FDR-corrected q-values. Differentially methylated regions (DMRs) were determined using DMRcate with default parameters (lambda = 1,000, C = 2, min.CpG sites = 2) (Peters et al., 2015). Pathway enrichment was performed using the package methylGSA (Ren and Kuan, 2019).
Angiotensin converting enzyme genotyping was performed using the same protocol as in Bigham et al. (2008). We extracted ACE CpG sites from the MethylationEPIC array to assess its methylation status independent from the epigenome-wide association analysis. We tested the relationship between ACE I/D status and SaO 2 separately for Kathmandu and Pheriche using linear modeling and adjusting for age and sex. The relationship between ACE genotype and phenotypes, as well as ACE genotype and ACE DNA methylation, was tested using linear mixed modeling in R using the lmerTest package (Kuznetsova et al., 2017). The linear mixed model was adjusted for altitude, age, sex, and individual IDs. Plotting was performed using the ggplot2 package (Wickham, 2009).

Participant Demographics
Our study group included n = 15 participants of self-reported European ancestry, with 53% females and the average BMI of 22.60 (SD 2.36) at baseline. Participant characteristics can be found in Table 1.

Physiological Changes With Altitude Exposure
We detected significant physiological changes between altitude 1,400 m (day 0) and 4,240 m (day 7) (henceforth physiological variables are referred to as phenotypes in this manuscript) in arterial oxygen saturation (SaO 2 ), hemoglobin concentration [(Hb)], and end-tidal carbon dioxide partial pressure (P ET CO 2 ) ( Table 1). Briefly, we observed a significant increase in [Hb] and a significant decrease in SaO 2 and P ET CO 2 with increasing altitude. The physiological responses we have reported are expected at high altitude, i.e., lower arterial oxygen saturation due to decreased ambient PO 2 , a decrease in P ET CO 2 indicating an increase in alveolar ventilation, and higher [Hb], reflecting the body's physiological response to low-oxygen conditions by increasing hemoglobin production.
In order to detect biological pathways overrepresented among the significant CpG sites from the analysis of differential methylation, we performed a pathway enrichment analysis using the methylgometh function in the R package methylGSA (Ren and Kuan, 2019). Ninety-five significant pathways were identified by methylgometh including the GO pathways regulation of glycolytic process (GO:0006110), regulation of hematopoietic stem cell differentiation ( GO:1902036), and regulation of angiogenesis (GO:0045765) (Supplementary Table 2). Other pathways of interest included brain development (GO:0007420), negative regulation of neuron differentiation (GO:0045665), and interleukin-1-mediated signaling pathway (GO:0070498).
We then tested for DMRs, i.e., contiguous regions in the genome that show differential methylation between phenotypes or groups. We used DMRcate (Peters et al., 2015) to find DMRs between low-(1,400 m) and high-altitude (4,240 m) samples. Using this approach, we identified 361 significant DMRs out of 657,408 possible DMRs (Supplementary Table 3). These included DMRs near/in genes associated with the HIF pathway: HIF1A and ENO1 (glycolytic enzyme), and the RAS pathway: ABL1, FGFR3, KRAS, RASA3, and RGL2.

Phenotype Associations
To determine if changes in DNA methylation could be driving acclimatization, we performed association testing between significant genome-wide methylation positions and phenotypes associated with high-altitude acclimatization. To do so, we focused our analysis on significant CpG sites identified in the DMP analysis (N = 2,873) and phenotypes that were significantly different between the groups ( Table 1) including SaO 2 , [Hb], and P ET CO 2 . Two CpG sites, cg16546681 (chr1:155244518, q-value = 0.01, β regression coefficient = 6.46) in the gene CLK2 and cg14548038 (chr9:140178418, q-value = 0.03, β regression coefficient = 4.73) upstream of the gene TOR4A, were significantly positively associated with SaO 2 (%). No significant associations were identified for [Hb] or P ET CO 2 after correcting for multiple comparisons.

ACE I/D, Oxygen Saturation, and DNA Methylation
We tested the relationship between ACE, a gene associated with high-altitude performance, and high-altitude phenotypes [SaO 2 , P ET CO 2 , (Hb)]. Individuals in this study were genotyped for the ACE I/D (rs4646994) polymorphism. We performed a genotypic test, wherein I/I and I/D genotypes were compared to D/D genotype, and identified a significant association between ACE genotype and SaO 2 . Individuals with genotypes I/D (β regression coefficient = 1.69, p-value < 0.01) and I/I (β regression coefficient = 1.85, p-value < 0.05) had significantly higher SaO 2 than individuals with the D/D genotype at 1,400 m (Kathmandu); the relationship was not significant for 4,240 m (Pheriche) (Figure 2). In an additive model, the I-allele was associated with increased SaO 2 (β regression coefficient = 1.03, p-value < 0.05) at 1,400 m; the relationship was also not significant for 4,240 m. In a dominant model, individuals who were either heterozygotes or homozygotes for the I-allele (grouped together) displayed higher SaO 2 (β regression coefficient = 1.71, p-value < 0.01) at 1,400 m; the relationship was approaching significance (β regression coefficient = 2.79, p-value = 0.09) for 4,240 m. Our results suggest that the dominant model, wherein individuals carrying either the I/D or I/I alleles have higher oxygen saturation than individuals carrying the D/D allele, is best suited to explain the relationship between SaO 2 ACE I/D in our study.
Average ACE methylation was positively associated with the I-allele when we tested the relationship using an additive model (β regression coefficient = 0.31, p-value = 0.03) ( Figure 2B). We also tested the relationship between individual ACE CpG sites and high-altitude phenotypes. ACE CpG sites, cg02040921 and cg09920557, were associated with SaO 2 (cg02040921: p-value = 0.04; cg09920557: p-value = 0.06). Increased methylation of ACE CpGs was associated with increased SaO 2 (Figures 2C,D). No significant associations were identified for [Hb] or P ET CO 2 .

DISCUSSION
The role of epigenetic change, including DNA methylation, in acclimatization to short-term hypoxia exposure is not well characterized. We aimed to fill this gap using genome-wide DNA methylation data from the same individuals measured at different altitudes during a trek to Everest Base Camp. We identified significant associations between genome-wide DNA methylation and short-term altitude exposure, among which were CpG sites and regions associated with HIF pathway, including HIF1A, and RAS pathway genes. We identified both a significant CpG position (DMP) and a DMR associated with hypoxia inducible factor 1A or HIF1A, which is a central gene in the body's hypoxic response (Slemc and Kunej, 2016). In normoxia, HIF1A is degraded via ubiquitination but is allowed to accumulate in hypoxic conditions. This allows its protein product to bind to a constitutively expressed HIF1B forming a heterodimer that activates downstream genes (Wenger, 2002). HIF1A activity is under epigenetic control in human cancer cells and hematopoietic cell lines (Walczak-Drzewiecka et al., 2010;Nguyen et al., 2013;Cimmino et al., 2019). Importantly, the HIF1A-associated DMR identified here overlaps with the promoter region of the gene, suggesting that methylation at this locus may be associated with changes in gene expression.
We found significant CpG sites associated with the RAS pathway, including ones in the genes ANGPT1 and RASA3 (RAS P21 protein activator 3). Angiopoietins 1 and 2 are regulated by HIF1, and ANGPT1 expression is associated with increased number of vessels without excessive permeability (Kelly et al., 2003). ANGPT1 can be activated and repressed by HIF1 in a cell-specific manner (Kelly et al., 2003). RASA3 (RAS P21 protein activator 3) is a Ras-GTPase activating protein that causes anemia and thrombocytopenia in mice when mutated  (Blanc et al., 2012). RAS pathway is another canonical hypoxiainduced pathway. RAS has been linked to blood pressure (Fontes et al., 1994), cardiovascular disease (Lee et al., 1993), and primary hypertension (Frossard et al., 1998). The role of RAS in hypoxia has been explored in association with high-altitude pulmonary edema via the regulation of the pulmonary vascular tone (Stobdan et al., 2011).
We also found significant DNA methylation changes in genes outside of canonical pathways implicated in highaltitude acclimatization (i.e., HIF and RAS), including significant DNA methylation changes in genes associated with cognitive impairment [ASH1L (de Ligt et al., 2012;Crawley et al., 2016;Xi et al., 2020) and TNIK (Coba et al., 2012;Anazi et al., 2016)]. Cognitive decline is a common side effect of highaltitude hypoxia (Regard et al., 1989;Yan, 2014;Gao et al., 2015) that becomes apparent 1-2 weeks after initial exposure (Bolmont et al., 2000) and may improve to some degree upon acclimatization (Heinrich et al., 2019). This timing of the cognitive decline is consistent with our study design wherein we identified methylation changes in ASH1L and TNIK after 1 week of high-altitude exposure. In addition to methylation differences in genes associated with cognitive function, we also identified changes in several genes associated with inflammation. These included CpG sites in the genes IL12B (Glas et al., 2012;Liu et al., 2012) and TRIM31 (Song et al., 2016;Wang et al., 2018).
We specifically focused on the RAS pathway gene ACE as it is centrally involved in circulatory homeostasis, and the ACE I/D polymorphism has been linked to endurance performance (Myerson et al., 1999), adaptation of highland resident/native populations (Qadar Pasha et al., 2001;Bigham et al., 2008), and performance at altitude (all those other citations). The ACE I allele is associated with lower ACE activity (Costerousse et al., 1993) and higher SaO 2 (Woods et al., 2002), potentially as a result of an increased HVR (Patel et al., 2003). We identified an association between ACE genotypes I/D and I/I with higher SaO 2 , which is consistent with previous research showing a significant relationship between ACE and SaO 2 (Woods et al., 2002;Bigham et al., 2008).
We found that the ACE I-allele was associated with higher average ACE methylation, which has been shown before in a study of birth weight and ACE (Rangel et al., 2014). Notably, the ACE I-allele is associated with lower serum and tissue ACE activity (Rigat et al., 1990;Costerousse et al., 1993;Woods et al., 2000). Since methylation is commonly associated with gene silencing, the association between ACE I-allele and higher DNA methylation suggests that ACE methylation may be involved in mediating decreased ACE expression in individuals with the I-allele.
Individuals at high altitude displayed increased [Hb] and decreased SaO 2 and P ET CO 2 compared to low altitude. We found two CpG sites, in the gene CLK2 and near the gene TOR4, that were associated with SaO2. CDC like kinase 2 or CLK2 suppresses PPARGC1A transcriptional activity on gluconeogenic genes (Sahu et al., 2019) and thus downregulates hepatic gluconeogenesis and glucose output. We found CLK2 methylation to be positively associated with SaO 2 , suggesting that CLK2 expression is potentially decreased in hypoxic conditions, given methylation is linked to gene repression. Interestingly, the CpG site in CLK2 is upstream of the gene PKLR that is significantly differentially methylated in high-compared to lowaltitude Quechua (Childebayeva et al., 2020). We also found a CpG site upstream of TOR4A (Torsin family 4 member A), which is associated with dystonia (Cascalho et al., 2017). Dystonia is linked to hypoxic exposure, more specifically cerebral anoxia/hypoxia (Kuoppamaki et al., 2002;Kern et al., 2016), and our finding might indicate a potential epigenetic mechanism playing role in the development of this condition.
Tissue types can show different methylation profiles across the body, and the degree to which they correlate varies by study design, type of sample, or age (Langie et al., 2017). For example, there is evidence of a low correlation between salivary and blood global DNA methylation (Godderis et al., 2015). Here, we analyzed saliva. Saliva is an attractive tissue for the analysis of DNA methylation in field studies given its relative ease of collection compared to blood or other tissues (Langie et al., 2017). By focusing on a singular tissue type, our results may be restricted to salivary tissue alone. However, salivary DNA methylation patterns have been shown to correlate with DNA methylation from blood (Thompson et al., 2013;Langie et al., 2016), intestinal mucosa (Hearn et al., 2019), and the brain (Smith et al., 2015). Furthermore, saliva panels have shown proteomic changes upon hypoxic exposure in cell cultures (Jain et al., 2020), suggesting the relevance of this tissue for analyzing the overall hypoxic response. Therefore, we suggest that our analysis of saliva is an important first step in identifying DNA methylation changes to acute hypoxia that may be relevant to other bodily tissues.
Overall, our data demonstrate that various pathways and systems are affected by exposure to high altitude, including the HIF pathway, RAS pathway, cognitive performance, and inflammatory systems. Moreover, we identified a significant association between SaO 2 and ACE I/D, and associations between ACE I/D and ACE methylation, further highlighting the connection between ACE and SaO 2 as well as the role of ACE in altitude acclimatization.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author AC, upon request.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Syracuse University Institutional Review Board (Protocol 18-006), University of Michigan Institutional Review Board (HUM00141118), Mount Royal University Human Research Ethics Board (Protocol 100012 and 101361), and Nepal Health Research Council (Protocol 109-2017). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
AC, TB, and AB: conceptualization. AC, TD, and JW: data curation. AC: formal analysis. TD and AC: funding acquisition. AC and TH: investigation. AC and AB: writing-original draft preparation. AB, AC, TB, TD, and TH: writing-review and editing. All authors contributed to the article and approved the submitted version.