Genome-Wide Epigenetic Signatures of Adaptive Developmental Plasticity in the Andes.

High-altitude adaptation is a classic example of natural selection operating on the human genome. Physiological and genetic adaptations have been documented in populations with a history of living at high altitude. However, the role of epigenetic gene regulation, including DNA methylation, in high-altitude adaptation is not well understood. We performed an epigenome-wide DNA methylation association study based on whole blood from 113 Peruvian Quechua with differential lifetime exposures to high altitude (>2,500) and recruited based on a migrant study design. We identified two significant differentially methylated positions (DMPs) and 62 differentially methylated regions (DMRs) associated with high-altitude developmental and lifelong exposure statuses. DMPs and DMRs were found in genes associated with hypoxia-inducible factor pathway, red blood cell production, blood pressure, and others. DMPs and DMRs associated with fractional exhaled nitric oxide also were identified. We found a significant association between EPAS1 methylation and EPAS1 SNP genotypes, suggesting that local genetic variation influences patterns of methylation. Our findings demonstrate that DNA methylation is associated with early developmental and lifelong high-altitude exposures among Peruvian Quechua as well as altitude-adaptive phenotypes. Together these findings suggest that epigenetic mechanisms might be involved in adaptive developmental plasticity to high altitude. Moreover, we show that local genetic variation is associated with DNA methylation levels, suggesting that methylation associated SNPs could be a potential avenue for research on genetic adaptation to hypoxia in Andeans.

Significance

High-altitude adaptation is a fascinating example of the human ability to adapt to marginal environments. The role and contribution of epigenetic modifications are not well understood in the context of this adaptation. We studied the potential role of epigenetic modifications in high-altitude adaptation in a group on indigenous Andean Quechua who are genetically adapted to life at high altitude. We found changes in DNA methylation, an epigenetic modification, associated with differential exposure to high altitude in genes involved in body’s response to hypoxia. Our findings suggest that epigenetic modifications might be involved in high-altitude adaptation.

Introduction

High-altitude regions exceed 2,500 m above sea level, and are characterized by decreased oxygen availability, high levels of UV radiation, and low ambient temperatures, making them a challenging environment for human life (Baker and Little 1976). High-altitude natives including Andeans, Tibetans, and Ethiopians display unique evolutionary adaptations allowing them to overcome hypoxic stress including population-specific genetic adaptations (Beall et al. 2010; Bigham et al. 2010; Simonson et al. 2010; Alkorta-Aranburu et al. 2012; Brutsaert et al. 2019). Among high-altitude adapted populations, Andeans are one of the most well characterized. They have inhabited the Andean Altiplano (average elevation of 12,000 ft) for up to 11,000 years, displaying unique physiological adaptations including increased chest circumference, elevated hemoglobin concentration, elevated arterial oxygen saturation, and a blunted (low) hypoxic ventilatory response (Bigham and Lee 2014). Genetic signatures of high-altitude adaptation have been reported in hypoxia-inducible factor (HIF) pathway genes that are central to oxygen sensing and regulation including EGLN1 and EPAS1 (Bigham et al. 2010; Foll et al. 2014; Eichstaedt et al. 2017; Brutsaert et al. 2019). Moreover, polymorphisms in EGLN1, the HIF gene, PRKAA1, and the renin–angiotensin system gene, ACE, have been associated with altitude-adaptive phenotypes, including the maximal oxygen uptake (VO2Max) during exercise, birth weight, and arterial oxygen saturation (SaO2), respectively (Bigham et al. 2008, 2014; Brutsaert et al. 2019). Together, these studies support evolutionary adaptation to high altitude among Andeans.

In addition to classic evolutionary adaptation, growth and development at high altitude contribute to altitude-adaptive traits among Andeans through the process of developmental adaptation. The term “developmental adaptation,” which also can be understood as adaptive developmental plasticity or irreversible adaptive plasticity (Storz and Scott 2019), was coined by Dr Roberto Frisancho to describe how exposure to high altitude during growth and development affected lung volume in Andeans (Frisancho 1970, 1977). Here, we use the term developmental adaptation throughout to refer to the process of adaptive developmental plasticity in Andeans. It is based on the hypothesis that adult phenotypes are affected by the environment during the developmental phase of early life (Frisancho et al. 1973, 1997; Brutsaert 2001). Importantly, adaptive developmental plasticity (Lafuente and Beldade 2019), or developmental adaptation, is distinct from immediate physiological plasticity, or acclimatization, to high altitude wherein sojourners to altitude experience physiological changes over days to weeks to overcome transient low-oxygen exposure. Acclimatization is reversible upon return to low altitude whereas developmental adaptation is not.

There are several examples of Andean high-altitude adaptive phenotypes where developmental adaptation plays a key role in the manifestation of the adult phenotype. Lifelong Andean residents of high altitude display greater lung volumes than lifelong sea-level residents of Andean ancestry (Frisancho et al. 1997; Kiyamu et al. 2012). Peruvian Quechua born and raised at altitude display higher SaO2 during exercise compared with Peruvian Quechua born and raised at sea level, but tested in altitude conditions (Bigham et al. 2008; Kiyamu, León-Velarde, et al. 2015). Peruvian Quechua residents of sea level with full developmental exposure to high altitude display significantly higher submaximal oxygen uptake (VO2peak) in simulated hypoxic conditions compared with low-altitude-born Peruvian Quechua (Kiyamu, Elías, et al. 2015). These studies support the idea that growth and developmental at high altitude play a pivotal role in Andean altitude-adaptive phenotypes specifically through the process of developmental adaptation. It is important to note that developmental adaptation among Andeans has a genetic component, but needs the environmental stimulus for the manifestation of the trait. This is illustrated by SaO2. Peruvian Quechua born and raised at altitude display higher SaO2 during exercise than Europeans who are born and raised at high altitude, and those same Europeans have similar SaO2 to European newcomers who are acclimatized to high altitude (Brutsaert et al. 2000).

The molecular mechanism underlying the environmental aspect of developmental adaptation to high altitude has not been identified. One possible means by which the environment can influence individual traits is through epigenetic modifications that change the expression of DNA (Frisancho 2009; Szyf and Bick 2013; Maccari et al. 2014; Julian 2017; Childebayeva et al. 2019). Epigenetics is the study of mitotically heritable and stable modifications to the DNA that do not involve changes to the nucleotide sequence. The most well-studied epigenetic mark is DNA methylation, which results from the addition of a methyl group to the cytosine in cytosine-guanine dinucleotides.

Early development is considered a critical window for epigenetic plasticity. Shortly after fertilization and before preimplantation, DNA methylation undergoes reprogramming wherein the maternal and paternal genomes are demethylated and a fetal methylome is established (Deng et al. 2009; Messerschmidt et al. 2014). During this critical window, DNA methylation is liable and sensitive to environmental exposures (Weaver et al. 2004; Dolinoy, Huang, et al. 2007; Nagy and Turecki 2012). DNA methylation patterns established prenatally as well as early in the postnatal period are associated with methylation patterns later in life (Dolinoy, Weidman, et al. 2007; Heijmans et al. 2008; Non et al. 2016), and thus may be involved in the organism’s phenotypic response to their environment. This phenomenon is known as the predictive-adaptive response (Gluckman and Hanson 2004; Jablonka and Raz 2009; Kuzawa and Thayer 2011), and it occurs when the parental environment during prenatal development as well as early peri/postnatal exposures can be used by the individual to “prepare” for future environmental conditions. Furthermore, epigenetic variation has been thought to contribute to the ability of an organism to adapt to a new environment when genetic variation is limited. In this case, epigenetic changes can potentially act as mechanisms of phenotypic plasticity (Jablonka 2013).

Several studies have shown that DNA methylation differences are associated with hypoxia exposure in various populations. First, our previous work (Childebayeva et al. 2019) demonstrated that Quechua individuals born and raised at low versus high altitude have significant differences in DNA methylation of the repetitive element LINE1 and the gene EPAS1, a well-known high-altitude adaptive gene that has been shown to be under genetic selection in Tibetans (Beall et al. 2010; Yi et al. 2010) and Andeans (Foll et al. 2014; Eichstaedt et al. 2017). Second, genome-wide DNA methylation differences in peripheral blood mononuclear cells have been identified between offspring of mothers with and without hypertensive pregnancy at altitude in La Paz, Bolivia (Julian et al. 2015). Third, males presenting with excessive erythrocytosis, a preclinical form of chronic mountain sickness, exhibit epigenetic differences compared with healthy matched controls (Julian 2017). Despite these advances in our understanding of epigenetic changes associated with high-altitude exposure, genome-wide DNA methylation patterns have not been studied in healthy Andeans to characterize the role of epigenetic change in altitude adaptation for this population.

Here, we hypothesized that early developmental exposure to high-altitude hypoxia in Andeans will be associated with differences in DNA methylation that could contribute to their pattern of adaptive developmental plasticity. We analyzed genome-wide DNA methylation signatures among 113 Peruvian Quechua who were recruited based on a migrant study design. This migrant study design, wherein individuals from the same population differ based on their early developmental and lifelong exposure to high altitude, allowed us to explore the potential influence of differential exposure to high-altitude hypoxia on the epigenome of Andean Quechua while controlling for underlying genetic effects on DNA methylation patterns.

Results

Study Design and Participant Characteristics

Andean Quechua are genetically and physiologically adapted to high altitude (Bigham et al. 2010; Kiyamu, León-Velarde, et al. 2015), including the individuals in our study group (Brutsaert et al. 2019). We used a migrant study design wherein three groups of Andean Quechua participants with different lifetime and developmental exposures to high altitude were recruited. They included: 1) high-altitude Quechua (HAQ)—highland individuals with lifetime exposure to high altitude who were born, raised, and resided at high altitude (n = 301). HAQ were recruited in the city of Cerro de Pasco, Peru (4,300 m); 2) migrant Quechua (MQ)—individuals with variable lifetime exposure to high altitude who were born in highland Peru (>3,000 m) and moved to low-altitude Lima, Peru (150 m) during their lifetime (n = 152). MQ were recruited in the city of Lima, Peru (150 m); and 3) low-altitude Quechua (LAQ)—individuals with no lifetime exposure to high altitude who were born and raised in Lima, Peru, but their parents and both sets of grandparents were of highland Quechua ancestry (n = 150). LAQ also were recruited in Lima, Peru (fig. 1). All participants were healthy males and females between the ages of 18 and 35 years whose parents and grandparents were of high-altitude origin, that is, they were born and raised at high altitude. We conducted a genome-wide principal component analysis (PCA) of genetic polymorphisms in the three groups comparing them with three world-wide populations from the 1000 Genomes Project including 90 East Asians, 60 Centre de Polymorphism Humain (CEPH) Europeans, and 60 Yorubans. We used microarray genotype data generated using the Affymetrix (Santa Clara, CA) Axiom Biobanking array that includes 600,000 markers. The PCA was performed on 328,260 autosomal SNPs with genotyping rates >95%, MAF > 0.0001, that passed Affy QC, and that were not in strong linkage disequilibrium (R2 < 0.8). The PCA demonstrated that three groups form a population cluster based on their genetic variation (supplementary fig. 1, Supplementary Material online).

. 1

Study design. (A) Map of Peru depicting high-altitude, Cerro de Pasco (4,300 m), and low-altitude, Lima (150 m), participant recruitment locations. (B) Genome-wide DNA methylation was determined in 113 Peruvian Quechua with differential exposures to high-altitude hypoxia over their lifetime. High-altitude Quechua (HAQ) were born, raised, and reside in or near Cerro de Pasco, Peru, Migrant Quechua (MQ) were born and raised at high altitude (>3,000 m), but moved to low altitude Lima (0 m) during lifetime, and low-altitude Quechua (LAQ) were born, raised, and reside at sea level in Lima, Peru. Number of participants analyzed with the EPIC chip is denoted in the parentheses.

The three study groups differ by: 1) prenatal and early postnatal exposure to hypoxia, what we call here developmental exposure, 2) lifelong high-altitude exposure, and 3) high-altitude acclimatization status. HAQ had full developmental exposure and were fully acclimatized to high altitude. MQ were developmentally exposed to hypoxia, but were not acclimatized to high altitude as residents of low altitude. MQ study participants varied in their degree of postnatal exposure to altitude depending on the age of down-migration to low altitude (on an average, 14.30 ± 6.62 years). LAQ had no developmental exposure to hypoxia and were not acclimatized to high altitude. These differences allowed us to identify DNA methylation differences associated with exposure to high-altitude hypoxia at different timepoints: early developmental exposure versus no exposure (prenatal and perinatal exposure) (LAQ vs. MQ comparison), lifelong exposure versus no exposure (LAQ vs. HAQ comparison), and lifelong high-exposure versus developmental exposure (HAQ vs. MQ). (See fig. 1 for a summary of the high-altitude exposure differences among the three groups.) It is important to note that the participants in this study are genetically undifferentiated (supplementary fig. 1, Supplementary Material online), and thus methylation differences observed between the groups are not expected to reflect underlying genetic variation.

We randomly selected a subset of participants (n = 113 from 603 total) from each group (HAQ = 38, MQ = 39, and LAQ = 36) for genome-wide DNA methylation analysis. Participant characteristics are provided in table 1. Consistent with previous research, the three study groups were significantly different (P < 0.05) by height, weight, fractional concentration of exhaled nitric oxide (FeNO), and hemoglobin (Hb) concentration. FeNO is a marker of respiratory nitric oxide (NO) synthesis. NO is a potent vasodilator whose levels decrease in hypoxic conditions given oxygen is essential for its production (Le Cras and McMurtry 2001). Body mass index (BMI), sex, age, forced vital capacity (FVC), and arterial oxygen saturation (SaO2) at rest measured at altitude (HAQ) or in simulated altitude conditions (LAQ and MQ) were not significantly different between the groups. Average age of down-migration for MQ was 14.30 ± 6.62 years of age. HAQ were recruited in Cerro de Pasco, a major mining center of Peru with high levels of inorganic lead in the soil as well as high whole-blood lead levels (WBLL) reported among women and children living there (Conklin et al. 2008; van Geen et al. 2012). Individuals recruited in Cerro de Pasco (HAQ) had significantly higher WBLLs (4.91 ± 1.72 µg/dl) compared with the individuals recruited in Lima (LAQ [2.11 ± 0.67 µg/dl] and MQ [1.99 ± 0.61 µg/dl]) (β = 3.05, P value <2e-16) (table 1). We corrected for WBLLs in all statistical models to adjust for the effect of inorganic lead on DNA methylation, given evidence for the effect of this heavy metal on DNA methylation (Goodrich et al. 2016).

Table 1

Participant Characteristics

Variable	All (n = 113)	Low-altitude Quechua (LAQ, n = 36)	Migrant Quechua (MQ, n = 39)	High-altitude Quechua (HAQ, n = 38)
% Female	50.44	52.78	46.15	52.63
Age (years)	24.38 ± 4.76	23.78 ± 4.25	25.23 ± 5.44	24.08 ± 4.46
Age migration (years)	NA	NA	14.30 ± 6.62	NA
Height (cm)*	158.91 ± 7.84	162.38 ± 7.16	158.08 ± 7.29	156.47 ± 8.02
Weight (kg)**	63.51 ± 10.28	67.73 ± 11.51	63.25 ± 10.40	59.79 ± 7.23
BMI	25.13 ± 3.52	25.61 ± 3.60	25.28 ± 3.73	24.51 ± 3.21
WBLL (μg/dl)***	3.08 ± 1.83	2.11 ± 0.67	1.99 ± 0.61	4.91 ± 1.72
FeNO (ppb)***	18.79 ± 12.55	20.83 ± 10.09	22.84 ± 15.27	17.60 ± 9.25
Hemoglobin (g/dl)***	15.04 ± 2.35	13.71 ± 1.49	13.88 ± 1.50	17.49 ± 1.72
FVC (l)a	4.44 (4.31, 4.57)	4.08 (3.88, 4.28)	4.37 (4.25, 4.49)	4.66 (4.45, 4.87)
SaO2 (%)	89.42 ± 3.12	89.62 ± 3.13	89.88 ± 3.21	89.08 ± 2.97

Note.—Data are means ± SD, 95% confidence intervals (CI) are in brackets. NA, not applicable; WBLL, whole blood lead level; BMI, body mass index; FeNO, fractional exhaled nitric oxide; FVC, forced vital capacity; SaO2, arterial oxygen saturation measured at rest at altitude (HAQ) or simulated altitude conditions (LAQ and MQ).

FVC mean adjusted for height using the ggeffects package in R (Lüdecke 2018).

P < 0.05 for LAQ versus MQ and LAQ versus HAQ based on t-test;

P < 0.05 between LAQ versus HAQ;

P < 0.05 for LAQ versus HAQ and MQ versus HAQ.

We generated DNA methylation data for 850,000 CpG sites using the Illumina (San Diego, CA) Infinium MethylationEPIC BeadChip assay. After quality control filtering of the probes, we analyzed 759,154 CpG sites for methylation differences between study groups (fig. 2). To identify methylation differences associated with lifelong exposure to high altitude, we compared LAQ with HAQ. To identify methylation differences associated with early developmental exposure, we compared LAQ with MQ as well as HAQ to MQ. For all comparisons, we used the following linear model in limma: DNA methylation ∼ Group + Sex + Age + Lead levels + PC1 (blood cell types) + chip batch.

. 2

Manhattan plot of DMPs. Both LAQ versus MQ and LAQ versus HAQ analyses plotted here. The −log10(P value) is shown on the y axis for each CpG site according to its genomic location as indicated on the x axis. The two significant DMPs are marked by their CpG site probe ID and their gene IDs. The dotted black line indicates P value <5 × 10−5 and the solid black line indicates the P value <5 × 10−8.

Differentially Methylated Positions and Regions Associated with Lifelong High-Altitude Exposure

We identified one significantly hypermethylated differentially methylated position (DMP), cg02269207, at q value <0.05 and genome-wide inflation factor λ = 1.23 between LAQ and HAQ, after correcting for multiple comparisons. This DMP was located in intron 1 of the gene CCDC64/BICDL1 (table 2 and  fig. 2). Differentially methylated regions (DMRs) can be considered more representative of biologically relevant epigenetic changes than single methylation sites (i.e., DMPs), given they reflect several CpG sites in close proximity changing in the same direction. Moreover, CpG islands potentially have more biological relevance than individual sites, since they are commonly associated with promoter regions, and thus can have a direct influence on gene expression (Saxonov et al. 2006). Therefore, we determined DMRs using comb-p (Pedersen et al. 2012). We identified 26 significant DMRs after Sidak’s correction for multiple testing between LAQ and HAQ (supplementary table 1, Supplementary Material online). These included DMRs in genes involved in anaerobic glycolysis (LDHC) (Goldberg et al. 2010; Odet et al. 2013), hematopoiesis (MEIS) (Simsek et al. 2010; Unnisa et al. 2012), and the remethylation pathway (BHMT) (Blaise et al. 2005).

Table 2

Statistically Significant DMPs Identified by EWAS

Study Group Comparison	CpG Site	Gene	Position (HG19)	Relation to Islanda	Gene Region Featurea	450K Enhancer	P Valueb	q Value	% Δ (High–Low)
LAQ vs. MQ	cg11428724	PAX7	chr1:18957632	Island	5′-UTR/1st exon	TRUE	8.77E-09	6.66E-03	4.02
LAQ vs. HAQ	cg02269207	CCDC64/BICDL1	chr12:120433323	OpenSea	Body of gene	—	5.89E-08	4.47E-02	5.94

Note.—TSS, transcription start site.

Based on annotation from Illumina.

Genome-wide inflation factor λ = 1.1 for the altitude of birth and 1.23 for the altitude of recruitment.

Differentially Methylated Positions and Regions Associated with Early Developmental High-Altitude Exposure

We identified a significantly hypermethylated DMP, cg11428724, at q value <0.05 and genome-wide inflation factor λ = 1.10 between LAQ and MQ (table 2 and  fig. 2) located in a CpG island found in exon 1 of the gene PAX7. We also identified 36 significant DMRs between LAQ and MQ (supplementary table 2, Supplementary Material online). These included DMRs in genes encoding antioxidant enzymes (SOD3), genes associated with blood pressure (EPHA6, GPR35) (Min et al. 2010; Mackenzie et al. 2011; Kim et al. 2017; Divorty et al. 2018), and high-altitude pulmonary edema (MACROD2) (Li et al. 2017), as well as HIF pathway target genes (SLC7A6, PKLR) (Benita et al. 2009). Seven significant DMRs overlapped between the LAQ versus MQ and LAQ versus HAQ analyses. These DMRs showed increased methylation in MQ or HAQ compared with LAQ. No statistically significant DMPs were identified between MQ and HAQ (data not shown). Therefore, we did not look for significant DMRs between these two groups. To test for significant associations between years lived at high altitude and DNA methylation, we tested for an association between genome-wide DNA methylation and the number of years an individual lived at high altitude, accounting for Group, Sex, Age, Lead levels, PC1 (blood cell types), and chip batch. No significant associations were identified (supplementary table 3, Supplementary Material online).

Given SOD3’s role in the hypoxia response, we validated this DMR in a larger cohort of Peruvian Quechua (n = 498, where HAQ = 240, MQ = 124, and LAQ = 134). SOD3 encodes the superoxide dismutase-3 enzyme that converts reactive oxygen species (ROS) superoxide anion to H2O2 and O2, and is a major extracellular scavenger of superoxide anion (Zelko et al. 2002; Suliman et al. 2004). The DMR is composed of four CpG sites located just upstream of the SOD3 promoter region (fig. 3). We targeted four CpG sites (cg03063309, cg16037589, cg17923553, cg11372436) that were identified in our differential methylation analysis using quantitative pyrosequencing. Based on the pyrosequencing results, SOD3 methylation was significantly increased in MQ compared with LAQ (β = 1.17, P = 0.04, fig. 3), replicating the results of our epigenome-wide association study (EWAS) analysis. Furthermore, SOD3 methylation was significantly increased in HAQ compared with LAQ (β = 1.76, P value <0.01), which we did not observe in our EWAS of lifelong high-altitude exposure. This is most likely the result of having increased statistical power to detect small differences in a sample of n = 374 individuals compared with n = 74. MQ and HAQ were not significantly different.

. 3

SOD3 DMR. (A) The SOD3 DMR between LAQ and MQ. Percent DNA methylation is represented on the x axis for each CpG site. Genetic coordinates for chromosome 4 shown on the y axis correspond to hg19 build. The SOD3 DMR (green box) between MQ+HAQ (blue line) and LAQ (orange line) is composed of four significant CpG sites marked with an *. The DMR lies just upstream of the SOD3 promotor indicated by the pink arrow. (B) Boxplot of SOD3 methylation by altitude of birth for each CpG position included in the DMR. CpG sites 1 through 4 are shown across the top. DNA methylation level is shown on the y axis. Significance levels between the study group comparisons are shown for each CpG site separately. LAQ, low-altitude Quechua; MQ, migrant Quechua; HAQ, high-altitude Quechua; NS, not significant; **P value <0.01; ***P value <0.001.

DMPs Associated with Altitude-Adaptive Phenotypes

To explore the impact of DNA methylation on altitude-adaptive phenotypes, we performed associations between genome-wide DNA methylation levels and FeNO, Hb, BMI, FVC, and SaO2 at rest. We chose these phenotypes because of known associations with high-altitude exposure and high-altitude adaptation (Brutsaert et al. 1999; Beall et al. 2001). We found 40 significant DMPs associated with FeNO (supplementary table 4, Supplementary Material online). Significant DMPs include ROBO2 involved in cardiomyocyte reoxygenation (Li et al. 2014) and P4HA1 upregulated by HIF-1α in breast tumors (Gilkes et al. 2013; Gilkes and Semenza 2013). We identified six DMRs associated with FeNO (supplementary table 5, Supplementary Material online). Genes in two of these DMRs, LOXL3 and P4HA1, have been shown to be hypoxia-inducible (Gilkes et al. 2013; Xie et al. 2017), and a third linked to RUNX1 has been shown to affect hematopoiesis (North et al. 2002; Growney et al. 2005; Friedman 2009). Importantly, FeNO is known to regulate hematopoiesis (Michurina et al. 2004). None of the DMP/DMRs associated with FeNO overlapped with the DMPs/DMRs associated with study group. Additionally, no significant associations were found between other phenotypes and DNA methylation (data not shown).

Associations between DNA Methylation and Genome-Wide Polymorphisms

Genetic variation is known to influence DNA methylation patterns (Zhang et al. 2010; Bell et al. 2011). Furthermore, local SNP variation has been shown to be associated with nearby DNA methylation variation (Bell et al. 2011). To better characterize the influence of SNP variation on DNA methylation, we determined associations between genome-wide SNPs and DNA methylation levels of two genes, SOD3 and EPAS1, for which we had targeted pyrosequencing data. We chose SOD3 given our present EWAS finding that MQ had significantly higher methylation at this locus compared with LAQ. Pyrosequencing data for this locus were generated as part of the current study. We selected EPAS1 based on three factors: 1. EPAS1 plays an essential role in the response to hypoxia. 2. It is a well-known high-altitude adaptive gene, and it has been shown to be under selection in Tibetans (Beall et al. 2010; Yi et al. 2010) and Andeans (Foll et al. 2014; Eichstaedt et al. 2017). Among Andeans, EPAS1 has been identified as a candidate locus for altitude adaptation among Aymara and Quechua speakers from Bolivia and Peru, respectively (Foll et al. 2014), and a novel EPAS1 SNP has been shown to be under early-stage selection among Argentine Collas (Eichstaedt et al. 2017). 3. In a previous study using the same study participants analyzed here, we demonstrated that EPAS1 methylation levels significantly differed between HAQ and LAQ (Childebayeva et al. 2019). EPAS1 pyrosequencing data analyzed here were generated for Childebayeva et al. (2019).

Autosomal SNP genotype data were generated using the Affymetrix (Santa Clara, CA) Axiom Biobanking Array for 460 Peruvian Quechua individuals from our cohort. SNPs that passed QC filters, including Affymetrix best practices thresholds, call rates >95%, and a minor allele frequency >0.01 were included in our analysis (n = 260,289 SNPs). We performed separate genome-wide association studies (GWAS) for EPAS1 and SOD3. Our models were adjusted for age and sex and corrected for multiple tests using the Benjamini–Hochberg correction. We identified two significant SNPs (q value <0.05, genome-wide inflation factor λ = 1), rs7579899 and rs2044456, that were associated with EPAS1 methylation (fig. 4A). Both SNPs reside in EPAS1 intron 1 located downstream (10–20 kilobases) from the EPAS1 methylation locus assayed by pyrosequencing. For both EPAS1 SNPs, the minor allele (rs7579899(A < G) and rs2044456(G < A)) was associated with lower EPAS1 methylation (table 3; fig. 4B and C). The effect size of EPAS1 SNPs rs7579899 (β = −1.27, P value <0.005; altitude of recruitment: β = −0.63, P value <0.05) and rs2044456 (β = −1.19, P value <0.005; altitude of recruitment: β = −0.65, P value <0.05) on EPAS1 methylation was greater than the effect size of altitude group. The allele frequencies of these two SNPs were not significantly different between the three altitude exposure groups (rs7579899 MAF LAQ = 0.30, MQ = 0.25, and HAQ = 0.28, rs2044456 MAF LAQ = 0.25, MQ = 0.22, and HAQ = 0.23). According to Ensembl (Yates et al. 2020), rs7579899 and rs2044456 EPAS1 SNPs are not enriched variants among Peruvians from Lima (table 3) nor are they in LD with SNP rs570553380 reported to be under early-stage selection in Andean Collas (Eichstaedt et al. 2017). No statistically significant SNPs were identified for SOD3 (genome-wide inflation factor λ = 1) (table 3). However, of interest, the second most significant SNP in this analysis resided 145 kb upstream of the SOD3 gene.

. 4

GWAS with EPAS1 methylation. (A) Manhattan plot of GWAS results testing for the effect of SNP genotype on EPAS1 methylation. Models were adjusted for age and sex. The two significant SNPs after correcting for multiple tests are indicated by their dbSNP IDs. Both significant SNPs are located EPAS1 intron 1. Boxplots of EPAS1 methylation by SNP genotype for the two, significant SNPs, (A) EPAS1 rs7579899 and (B) EPAS1 rs2044456. ***P value <0.001; **P value <0.01; *P value <0.05.

Table 3

GWAS Results for EPAS1 and SOD3

SNP	Location (hg19)	Nearest Genes	Type	MA	MAF	1KGa	AMRa	EURa	EASa	AFRa	n	β effect on DNA Methylation	P Value	q Value
GWAS with EPAS1 (λ = 1)
rs7579899	2:46537604	EPAS1	Intronic	A	0.28	0.57	0.37	0.38	0.81	0.66	454	−1.269	2.71E-09	0.0007
rs2044456	2:46546316	EPAS1	Intronic	G	0.23	0.28	0.27	0.3	0.56	0.06	454	−1.186	2.42E-07	0.0315
GWAS with SOD3 (λ = 1)
rs10480851	7:134285060	—	Intergenic	G	0.01	0.15	0.05	0.03	0.02	0.43	460	−7.875	5.70E-07	0.1484
rs10939004	4:24652142	SOD3	Intergenic	C	0.34	0.21	0.25	0.07	0.11	0.42	459	−1.486	2.06E-06	0.1911

Note.—MA, minor allele; MAF, minor allele frequency; 1KG, 1,000 genomes; AMR, American; EUR, European; EAS, East Asian; AFR, African.

Ensembl 1,000 genomes and gnomAD allele frequencies.

Accelerated Epigenetic Aging at High Altitude

The epigenome ages across an individual’s lifespan in a process known as epigenetic aging. An individual’s epigenetic age, estimated by an epigenetic clock algorithm based on DNA methylation data, can serve as a marker of the epigenetic maintenance system as well as of premature biological ageing. To establish how high altitude is associated with the epigenetic clock, we estimated epigenetic age using the bioinformatic method developed by Horvath (Horvath 2013). We determined the degree of epigenetic aging (ΔDNAmAge) as the difference between Horvath epigenetic age and chronological age. ΔDNAmAge was significantly higher in HAQ individuals compared with LAQ (β = 2.43, P = 0.005, fig. 5) and approached significance in HAQ compared with MQ (β = 1.39, P = 0.10). LAQ and MQ ΔDNAmAge were not significant.

. 5

Boxplot of ΔDNAmAge by group. ΔDNAmAge (difference between epigenetic age and calendar age) is shown across the three study groups where LAQ, low-altitude Quechua, MQ, migrant Quechua, and HAQ, high-altitude Quechua. HAQ have significantly higher Δage (β = 2.43, P value = 0.005) compared with LAQ when adjusted for chronological age, sex, batch, and the first principal component of blood cell types. **P value <0.01; +P value <0.10; NS P value >0.10.

Discussion

The role of epigenetic change, including DNA methylation, in high-altitude adaptation is not well characterized. We performed an EWAS in a cohort of Andean Quechua with differential exposure to high-altitude hypoxia. We identified: 1) DMPs and DMRs associated with lifelong high-altitude exposure, 2) DMPs and DMRs associated with early developmental high-altitude exposure, 3) DMPs associated with altitude-adaptive phenotypes, and 4) associations between genetic variation and local DNA methylation. Our results are the first to demonstrate genome-wide DNA methylation differences associated with differential high-altitude exposure among Andeans. These findings further support the notion that epigenetics may play a role in adaptive plasticity as a potential mediator of proximal environmental effects on the body (Jablonka and Lamb 2006). Furthermore, they support the notion that early developmental exposures can have a persistent impact on DNA methylation patterns, and add to the growing field of developmental origins of health and disease (DoHAD) (Waterland and Michels 2007; Kubota et al. 2015).

We identified one significant CpG site and 26 significant DMRs associated with lifelong exposure to high altitude (LAQ vs. HAQ). The DMP was located in the coiled-coil domain containing 64 (CCDC64) gene, also known as BICD family-like cargo adaptor 1, involved in neuronal development by regulating the transport of secretory vesicles. CCDC64 was identified as a candidate for genetic adaptation to high altitude among Tibetans (Simonson et al. 2010) and has been proposed to be under convergent adaptation to high altitude in Asia and South America (Foll et al. 2014). Among the 26 significant DMRs, MEIS, BHMT, and LHDC are potentially associated with high-altitude adaptation. Myeloid ecotropic viral integration site 1 (MEIS1) plays an important role in hematopoiesis, and is involved in the oxidant defense of hematopoietic stem cells (Kocabas et al. 2012). In mice, Meis1 regulates Hif-1a transcription, and its loss is associated with downregulation of Hif-1a and Hif-2a in hematopoietic stem cells (Simsek et al. 2010). Betaine-homocysteine S-methyltransferase (BHMT) is involved in the remethylation pathway of homocysteine into methionine, and has been implicated in cardiovascular disease (Weisberg et al. 2003). Lactate dehydrogenase C (LDHC) is a testis-specific glycolytic enzyme involved in anaerobic glycolysis (Goldberg et al. 2010; Odet et al. 2013). Changes to glucose metabolism genes are of interest in studies of high-altitude adaptation, and a shift from fatty acid oxidation to glucose oxidation and glycolysis has been proposed to be an adaptive mechanism at high altitude (Ge et al. 2012). This finding further supports the notion that glucose metabolism alterations may play a role in high-altitude adaptation, and supports the interpretation that LDHC DNA methylation is reprogrammed in the testis later in life.

We identified a single significant DMP and 36 significant DMRs associated with developmental exposure to high altitude (LAQ vs. MQ). The DMP was located in the gene PAX7 involved in skeletal development (Schafer et al. 1994; Relaix et al. 2006). Among high-altitude native populations, changes to skeletal muscle have been shown including decreased levels of muscular damage from ROS, lower intracellular lipid composition, and higher counts of type I fibers, as well as higher myoglobin levels (Hoppeler and Vogt 2001; Gelfi et al. 2004). Each of these changes suggests an adaptive response to low-oxygen conditions. Of the 36 DMRs, PKLR and SOD3 are potentially involved in developmental high-altitude adaptation. Pyruvate kinase L/R (PKLR) encodes R/L-type pyruvate kinase involved in red blood cell production and maintenance. It is under selection in Tibetans and high-altitude Mongolians (Yi et al. 2010; Xing et al. 2013; Wuren et al. 2014), has been associated with hereditary hemolytic anemia (Kanno et al. 1992; Huang et al. 2004), and is upregulated in the rat fetus in hypoxic conditions (Kanno et al. 1992; Huang et al. 2004). SOD3 promotes the expression of erythropoietin (EPO) in the kidney, which in turn increases red blood cell production. Moreover, SOD3 knockout mice show a lagging hematopoietic response to hypoxia (Suliman et al. 2004). We found increased methylation of SOD3, which we expect to be associated with decreased expression of its protein, since the region we assessed is upstream of the SOD3 promoter. This is in line with previous research demonstrating decreased levels of SOD3 among high-altitude natives of India, the Ladakhi, compared with acclimatized lowlanders (Padhy et al. 2016).

Together, our results support the hypothesis that lifelong high-altitude exposure as well as pre/and early postnatal (developmental) exposures to high altitude are associated with differences in DNA methylation. On an average, high-altitude exposure resulted in hypermethylation. However, both hyper- and hypomethylation were observed among HAQ and MQ when compared with LAQ. Seven of the DMRs were significant for both the HAQ/LAQ and the LAQ/MQ analyses. Two of the overlapping DMRs were associated with the genes CRTAC1 and TRIOBP linked to idiopathic pulmonary fibrosis. This condition is characterized by irreversible fibrosis of the lung, and known to be influenced by hypoxia in its etiology (Tzouvelekis et al. 2007; Huang et al. 2014; Gangwar et al. 2017; Aquino-Galvez et al. 2019; Wang et al. 2019).

Previous studies have proposed that NO may be involved in high-altitude adaptation given its role in regulating blood pressure and oxygen delivery (Hoit et al. 2005; Beall et al. 2012). However, studies disagree on whether Andeans have elevated (Beall et al. 2001) or decreased (Ghosh et al. 2019) levels of NO at high altitude. Overall, studies agree that NO levels are relevant for high-altitude adaptation, since higher exhaled NO is associated with increased blood vessel diameter and elevated blood flow (Hoit et al. 2005). We found 40 DMPs and six DMRs significantly associated with FeNO among Peruvian Quechua. Several significant DMRs including LOXL3, P4HA1, and RUNX1, have been associated with angiogenesis (Gilkes et al. 2013; Xie et al. 2017) and hematopoiesis (North et al. 2002; Growney et al. 2005; Friedman 2009), both of which are known to be regulated by NO (Michurina et al. 2004).

It is well known that local genetic variation (Zhang et al. 2010; Bell et al. 2011), as well as variation in one-carbon metabolism SNPs (Childebayeva et al. 2019) can effect DNA methylation. This is the first study to demonstrate that local variation in the EPAS1 gene at two SNPs, rs2044456 and rs7579899, is associated with EPAS1 methylation. Neither SNP’s MAF was significantly different between the three study groups. Therefore, differences in allele frequencies between the groups were not driving the association. Furthermore, the SNPs were not enriched in Andeans compared with 1000 Genomes Populations, and they were not in LD with the coding region SNP previously identified as under early-stage selection in Andean Collas (Eichstaedt et al. 2017). Together, this indicates that they are not under selection in Andeans and not linked to variation that has been shown to be under selection. Our finding suggests that genetic variation can influence nearby DNA methylation regions, highlighting the importance of underlying genetic variation in determining DNA methylation levels and emphasizing the need to consider this variation when comparing DNA methylation signatures within and between populations. However, the mechanism by which local variation can influence DNA methylation remains unclear.

We found that accelerated ΔDNAmAge was associated with being born at high altitude and potentially with lifetime high-altitude exposure. Life history theory predicts that organisms will prioritize faster growth and earlier reproduction when environmental conditions during the critical early period of development signal adversity. When this same period of development signals few indicators of adversity, life history theory predicts a “slow” life history strategy wherein slower growth and later reproduction are prioritized. Species with faster life histories exhibit a faster epigenetic clock (Lowe et al. 2018). Therefore, accelerated epigenetic aging may be a mechanism through which a “fast” life history strategy proceeds. In fact, accelerated epigenetic aging has been associated with increased stress in early life as well as with phenotypic outcomes including increased height for age, elevated diurnal cortisol, among others (Zannas et al. 2015; Davis et al. 2017; Suarez et al. 2018). We have found that epigenetic age is accelerated at high altitude, thus supporting the hypothesis that hypoxia places a strain on the cellular machinery, which is a known phenomenon (Semenza 2011).

Our study benefited from a migrant study design. An important feature of this design was that all of the study participants’ parents and grandparents were born and raised at high altitude. This held environmental exposure to high altitude in the parental and grandparental generations constant across all study participants. In so doing, we were able to control for the effects of intergenerational inheritance, wherein maternal exposures (F0) will affect the fetus (F1) and the gonads of the fetus (F2) (Lacal and Ventura 2018). Therefore, differences in methylation between study groups were the result of early developmental exposures of the fetus, and not the exposures of previous generations.

There are several limitations in our study. First, a potential explanation for the differences in DNA methylation that we observed here is a dose effect of high-altitude exposure given HAQ have full high-altitude exposure, MQ have various degrees of high-altitude exposure, and LAQ have no high-altitude exposure. Dose effect in this case cannot be separated from developmental programming. However, early developmental exposure is hypothesized to have a more significant effect on DNA methylation patterns than dose effect as a result of DNA methylation reprogramming (Maccari et al. 2017). Our results support this hypothesis. We did not observe significant DMPs between MQ and HAQ, and did not identify significant DMPs when DNA methylation was tested for an association with overall years lived at high altitude. However, our inability to detect significant differences between MQ and HAQ as well as significant DMPs with years at altitude could be the result of small sample size in the MQ group.

Second, we used DNA methylation data from whole blood, which is composed of various cell types, each with its own methylation landscape. Differences in methylation seen here could be a result of differences in cell-type composition between high and low altitudes. To mitigate this effect on our results, we controlled for blood cell type bioinformatically following the procedure proposed by Houseman et al. (2012).

Third, lead exposure may be associated with DNA methylation (Goodrich et al. 2016). Individuals recruited in Cerro de Pasco had higher levels of lead exposure than individuals recruited in Lima. This difference in toxicant exposure could have impacted DNA methylation levels. To control for this difference, our statistical models of differential methylation controlled for lead exposure. Moreover, differences in SES, prenatal health, psychosocial, or other physical exposures not collected as part of the study could be associated with DNA methylation patterns identified between the groups. It is also important to note that in some cases, such as with SOD3, we observed a very small difference in DNA methylation between the groups <2%. Small effect size differences are common in epigenetic studies, and can still be important given % methylation difference refers to the number of cells in the sample that have a different methylation status (methylated vs. unmethylated) at a specific CpG (Breton et al. 2017).

Fourth, we did not recruit a fourth developmental group to our migrant design study, that is, Peruvian Quechua individuals born at low altitude who moved to high altitude, since it is rare for Peruvians from generally more urban low-altitude locations to move to high-altitude locations that tend to be more rural. Data from the fourth developmental group could help characterize epigenetic differences between altitude exposure groups. Similarly, we did not characterize DNA methylation in a low-altitude Indigenous American population. Therefore, we were unable to compare LAQ with a low-altitude population with similar genetic ancestry to understand if LAQ exhibit unique DNA methylation patterns.

Finally, it is important to stress that the DNA methylation differences that we identified may be adaptive, maladaptive, or neutral. A more detailed exploration of DNA methylation differences together with RNA, protein, as well as phenotypic data will allow us to understand whether DNA methylation differences do in fact lead to changes in gene expression that affect altitude-adaptive phenotypes.

Conclusions

Andean high-altitude adaptation has been hypothesized to involve both genetic and epigenetic factors (Frisancho 2009; Julian 2017). We performed an EWAS on 113 individuals of Peruvian Quechua ancestry who were recruited using a migrant study design allowing us to separate the effects of developmental altitude exposure from lifetime altitude exposure. We identified DMPs and DMRs associated with birth and early development at high altitude, including genes involved in glycolysis (LDHC), blood pressure (EPHA6, GRP35), red blood cell production (PKLR, MEIS1), HIF pathway (SOD3), and skeletal muscle growth (PAX7). Moreover, we found DMPs and DMRs associated with FeNO. We performed separate genome-wide association studies between genome-wide polymorphism data and DNA methylation levels of EPAS1 and SOD3. We found that local genetic variation was significantly associated with DNA methylation levels for EPAS1 and approaching significance for SOD3, in some cases greater than the exposure of interest. Lastly, we found that individuals born at high altitude age faster compared with individuals born at low altitude suggesting that hypoxia influences cellular machinery at large.

Materials and Methods

Study Population

Our analysis included individuals from three high-altitude exposure groups: 1) high-altitude Quechua (HAQ), migrant Quechua (MQ), and low-altitude Quechua (LAQ) (fig. 1). Participants were unrelated, healthy, nonpregnant/lactating, nonsmokers, between 18 and 35 years old (table 1). All participants were of self-reported Quechua ancestry. Participants provided last names and birth locations for their parents and grandparents as another way to ascertain ancestry. At the time of enrollment, Hb, height, weight, BMI, and phenotypes associated with high-altitude adaptation, including forced vital capacity (FVC), FeNO, SaO2, and others were collected. Hb level was measured using a HemoCue Hb201+ analyzer (Angelholm Sweden). SaO2 was measured using a pulse oximeter. In LAQ and MQ participants, SaO2 was taken under simulated hypoxic conditions. FeNO was measured using a portable, hand-held nitric oxide analyzer, the NIOX MINO (AerocrineAB, Solna, Sweden). FVC was measured with a Vacumed Dry Seal Rolling Spirometer using disposable mouth pieces. The best of three expirations was selected for analysis. All study participants were screened for anemia using altitude specific cutoffs. Given that lead exposure has been shown to associate with DNA methylation (Wright et al. 2010; Senut et al. 2012, 2014), we measured WBLLs in a subset of participants using graphite furnace atomic absorption spectrometry.

At the time of enrollment, all study participants provided written informed consent approved by Syracuse University, Universidad Peruana Cayetano Heredia, and The University of Michigan Institutional Review Boards.

DNA Methylation Analyses

We generated DNA methylation data for ∼850,000 CpG sites using the Illumina Infinium MethylationEPIC BeadChip assay for 116 Quechua individuals between the ages of 18 and 35 years. These 116 DNA samples were randomly selected from a larger cohort of Peruvian Quechua participants (Brutsaert et al. 2019; Childebayeva et al. 2019) who were recruited using a migrant study design. We analyzed a total of 113 samples for genome-wide DNA methylation signatures using the Illumina Infinium MethylationEPIC (EPIC) BeadChip. We used the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, CA) to bisulfite convert each DNA sample following the standard protocol with alternative incubation conditions optimized for the Illumina Infinium MethylationEPIC BeadChip assay. Raw data were loaded into R in the form of idat files using the minfi package (version 1.20.0) (Aryee et al. 2014). QC was performed using minfi and ENmix to assess the quality of the bisulfite conversion and probe hybridization for each sample. In minfi, we plotted log2 unmethylated median intensities over log2 methylated median intensities, compared predicted sex based on X and Y methylation to recorded sample sex, and generated density bean plots of beta values. We also visualized methylated and unmethylated intensities of control probes using ENmix. Based on QC metrics, three HAQ samples were excluded from the analysis (supplementary fig. 2, Supplementary Material online). Two of the samples, self-reported as males, clustered with females, and the third sample was an outlier based on the PCA analysis of DNA methylation (supplementary fig. 2, Supplementary Material online). To check for potential sample mix-up, we extracted CpG sites that are SNPs using the MethylToSNP R package (LaBarre et al. 2019). The R package extracts CpG sites that are present in three discrete levels of methylation: fully methylated, fully unmethylated, and ∼50% methylation, which correspond to genotypes CC, TT, and CT. We then found an overlap between these CpG sites and SNPs present on the Affymetrix array. We compared five random sites between SNPs called using DNA methylation, and the Affymetrix gene chip array, and found perfect consistency. We consider this as proof that there was no sample mix-up during the analysis.

All statistical analyses were performed on 113 samples with reliable DNA methylation data.

Sample Preprocessing and Normalization

Data preprocessing and normalization are crucial steps in the array analysis. Background correction, dye-bias adjustment, and within-array normalization are commonly performed to adjust for technical variation in the sample (Marabita et al. 2013). We performed probe correction and normalization using preprocessQuantile followed by ChampBMIQ. Probes that were above the 0.00001 detection P value threshold in more than 5% of the samples, that is, failed to hybridize, were excluded from the analysis. Cross-reactive probes, probes associated with sex chromosomes, probes containing SNPs with MAF >5% at target CpG sites, single base extension sites of type I probes (n = 414), and in the body of the probe were removed (Chen et al. 2013; Pidsley et al. 2016). Probes containing SNPs were removed given that methylation levels detected by these probes can reflect the underlying genetic polymorphisms that could be misinterpreted as true signals (Chen et al. 2013). Samples were tested for batch effects using singular value decomposition (SVD) analysis in champ. SVD performed on n = 759,154 beta values. Pool_ID was determined to play a significant role on the variation of the genome-wide DNA methylation data based on the output from SVD, and corrected for in the linear modeling (supplementary fig. 3, Supplementary Material online).

Blood Cell-Type Heterogeneity Correction

DNA methylation is tissue and cell-type specific. Therefore, it is crucial to correct for cell-type composition, since cell-type variation can be correlated with phenotype. Blood is composed of seven to eight cell types: neutrophils, eosinophils, basophils, CD14+ monocytes, CD4+ T cells, CD8+T cells, CD19+ B cells, and CD56+ natural killer cells (Houseman et al. 2012; Teschendorff and Zheng 2017). Thus, we estimated cell counts in blood computationally using minfi following the method by Houseman et al. (2012). Houseman’s is a reference-based method that is based on an a priori defined DNA methylation reference for the cell types present in blood (Houseman et al. 2012). We performed a PCA on whole-blood cell counts, and chose only to include the epigenome-wide PC1 that was the most significant contributor to the blood-type variation to correct for differential methylation associated with cell-type heterogeneity in blood (supplementary fig. 4, Supplementary Material online). A linear regression analysis between PC1 from blood-cell types and hemoglobin levels was not significant (P value = 0.16) indicating that a change in hemoglobin at high altitude is not significantly associated with the composition of the DNA containing cells.

Differential Methylation Analysis

Fully processed M values were tested for differential methylation using the package limma functions lmFit and eBayes (Ritchie et al. 2015). The following model was used to test for the DMPs: DNA methylation ∼ Group + Sex + Age + Lead levels + PC1 (blood cell types) + chip batch. In this analysis, LAQ was coded as the reference group, and LAQ methylation was compared with MQ and with HAQ. We then reran this model with HAQ as the reference group in order to mimic the ANOVA with post hoc test in the regression framework. Using linear modeling allowed us to control for covariates in the differential methylation testing. The covariates to include in the model were determined based on the analysis of overall variation in DNA methylation based on singular value decomposition (SVD) analysis (supplementary fig. 3, Supplementary Material online). The resulting P values were adjusted for multiple testing using the false-discovery rate (FDR) following the Benjamini–Hochberg procedure (Hochberg and Benjamini 1990). DMPs were annotated with the R package IlluminaHumanMethylationEPICanno.ilm10b3.hg19 (Hansen 2017). We performed the association between DMPs and phenotypes for all samples combined. These models were corrected for group. The following model was used: DNA methylation ∼ Phenotype + Group + Sex + Age + Lead levels + PC1 (blood cell types) + chip batch. We tested height and weight as covariates, however, they were not significantly associated with genome-wide DNA methylation and were excluded them from the final model analysis.

Comb-p was used on the results of the limma analysis to test for DMRs (Pedersen et al. 2012). DMRs were defined as genomic regions with three or more CpG sites significant at unadjusted P value <0.05 within 3,000 bp of each other. All CpG sites with unadjusted P value <0.05 were used in this analysis. DMRs were annotated with the R package IlluminaHumanMethylationEPICanno.ilm10b3.hg19 (Hansen 2017). This annotation was used to determine whether a DMR was located in a CpG island. Sidak correction for multiple tests was performed as the default correction in comb-p (Pedersen et al. 2012).

Pyrosequencing

Quantitative pyrosequencing was performed to assess DNA methylation levels of SOD3 in n = 603 samples. EPAS1 methylation was available from Childebayeva et al. (2019). The following primer sequences were used to amplify the region of interest: SOD3 forward primer sequence: 5′-[biotin]-TTGTGTGTTGAAGGTTATTGGTTATAAG-3; SOD3 reverse primer sequence: 5-CCTCCTCTACCCCTCCCATTTTT-3; SOD3 sequencing primer: 5-ACCCCTCCCATTTTTA-3. We assessed the DNA methylation levels at four CpG sites in SOD3 using the Pyromark Q96 pyrosequencer (Qiagen, Valencia, CA), for more details, see Virani et al. (2012). Of the 603 samples analyzed using pyrosequencing, 104 either failed pyrosequencing or had high coefficient of variance (CV > 10) between duplicates and were excluded from the statistical analyses for SOD3.

Epigenetic Age

We estimated epigenetic age using the method developed by Horvath (2013). Our analysis included 336 out of the 353 CpG sites originally included by Horvath. Seventeen loci were not included on the Illumina MethylationEPIC array and therefore were unable to be included in our analysis. We used linear models to determine if epigenetic age (ΔDNAmAge, where ΔDNAmAge = Horvath Epigenetic Age − Chronological Age) was associated with high altitude. All models were adjusted for chronological age, sex, first blood cell-type principal component, and batch.

Microarray Genotyping

We generated microarray genotype data for all Andean Quechua participants using the Affymetrix (Santa Clara, CA) Axoim Biobanking array that features ∼610,000 markers including haplotype tagging SNPs, exonic SNPs, and eQTL variants. We performed statistical analyses on autosomal variants that passed Affy best practices, had call rates >95%, and minor allele frequency >0.01 (n = 260,289). Genome-wide association studies were performed on % methylation of EPAS1 and SOD3 in PLINK with sex and age as covariates, quantile–quantile plots, and Manhattan plots were made using qqman package in R (Turner 2018). We performed the Benjamini–Hochberg correction for multiple tests, and the genome-wide inflation factor lambda was calculated as median(chisq)/qchisq (0.5,1). A PCA was performed on genome-wide SNP data for the participants in this study together with three populations from 1,000 genomes (90 East Asians including 45 Han Chinese from Beijing and 45 Japanese from Tokyo, 60 Centre de Polymorphism Human [CEPH] Europeans, and 60 Yorubans from Ibadan, Nigeria) (supplementary figures 1A–C as in Childebayeva et al. [2019]).

Click here for additional data file.