The genetic architecture of adaptations to high altitude in Ethiopia

Although hypoxia is a major stress on physiological processes, several human populations have survived for millennia at high altitudes, suggesting that they have adapted to hypoxic conditions. This hypothesis was recently corroborated by studies of Tibetan highlanders, which showed that polymorphisms in candidate genes show signatures of natural selection as well as well-replicated association signals for variation in hemoglobin levels. We extended genomic analysis to two Ethiopian ethnic groups: Amhara and Oromo. For each ethnic group, we sampled low and high altitude residents, thus allowing genetic and phenotypic comparisons across altitudes and across ethnic groups. Genome-wide SNP genotype data were collected in these samples by using Illumina arrays. We find that variants associated with hemoglobin variation among Tibetans or other variants at the same loci do not influence the trait in Ethiopians. However, in the Amhara, SNP rs10803083 is associated with hemoglobin levels at genome-wide levels of significance. No significant genotype association was observed for oxygen saturation levels in either ethnic group. Approaches based on allele frequency divergence did not detect outliers in candidate hypoxia genes, but the most differentiated variants between high- and lowlanders have a clear role in pathogen defense. Interestingly, a significant excess of allele frequency divergence was consistently detected for genes involved in cell cycle control, DNA damage and repair, thus pointing to new pathways for high altitude adaptations. Finally, a comparison of CpG methylation levels between high- and lowlanders found several significant signals at individual genes in the Oromo.


INTRODUCTION
examined in different parts of the world. Specifically, the Ethiopian highlands offer a unique opportunity to study HA adaptation because individuals from distinct, but closely related ethnic groups have communities at HA and LA, thus allowing more informative genetic and phenotypic comparisons. In this study, we extended genomic analysis to two Ethiopian ethnic groups, Amhara and Oromo, with the goal of determining whether Tibetans and Ethiopian highlanders share the same adaptations and of elucidating the genetic bases of adaptive HA phenotypes in Ethiopia. We also measured genome-wide methylation levels to explore the contribution of epigenetic modifications to HA adaptations.

The Ethiopian Amhara and Oromo differ in adaptive phenotypes
We obtained phenotype data in two distinct, but closely related ethnic groups, the Amhara and the Oromo (Texts S1-S2; Figures S1-S5), that include communities of HA and LA residents. All individuals were born and raised at the same altitude where they were sampled. These samples allow comparing phenotypes across altitudes within ethnic groups as well as across ethnic groups. While historical records indicate that the Oromo have moved to HA only in the early 1500s [22,23], the Amhara have inhabited altitudes above 2500m for at least 5ky and altitudes around 2300-2400m for more than 70ky [24,25]. Therefore, sufficient time has elapsed for the Amhara to have evolved genetic adaptations to hypoxia.
As shown in Figure 1, the HA samples of both ethnic groups had higher Hb than the LA samples, however the Oromo had twice as much elevation in Hb as the Amhara. The elevation in Hb levels is particularly evident for the measurements in males, raising the possibility that other factors (e.g. menstrual cycle) in females affect the power to detect significant phenotypic differences between groups. With regard to O 2 sat, HA Amhara had a 5.6% lower O 2 sat compared to LA Amhara while HA Oromo had 10.5% lower O 2 sat than their LA counterparts. Therefore, we detected significant phenotypic differences not only between populations from the same ethnic group that live at different altitudes, but also across populations from closely related ethnic groups (Oromo and Amhara) that live at the same altitude. Given the low genetic divergence between these two ethnic groups at the genomewide level (mean F ST = 0.0098), the phenotypic differences between Amhara and Oromo highlanders are unlikely to be due to independent genetic adaptations in these ethnic groups; rather they are likely to reflect genetic adaptations that evolved in the Amhara, due to their longer residence at HA.
We also measured pulse and calculated arterial oxygen content, but these phenotypes did not show significant differences across ethnic groups or altitudes and were omitted from further analyses. For details on the phenotypic variation in Amhara and Oromo, see Text S3, Tables S1-S2 and Figure S6.

Genome-wide association signals in Ethiopia
To learn about the genetic bases of variation in Hb level and O 2 sat in Ethiopia, we tested for an association between SNP genotype in the 260 unrelated Ethiopian samples and Hb levels or O 2 sat. We considered the total Ethiopian sample (i.e. HA and LA Amhara and Oromo) as well as each ethnic group and each altitude separately (Figures S7-S15). No genome-wide significant signal was observed for either Hb levels or O 2 sat in the Oromo and in the total Ethiopian sample and for O 2 sat in the Amhara (Figures S7-S15 and Tables S3-S20). Likewise, no excess of low p-values was observed in these association analyses relative to null expectations obtained by permutations ( Figures S7-S15). In contrast, the Amhara showed an excess of low association p-values in the analysis of Hb levels compared to expectations obtained by permutations (Figure 2A), indicating a genetic contribution to variation in Hb levels. In addition, one SNP (rs10803083) on chromosome 1 was associated with variation in Hb levels (p= 4.96x10 -8 ) in Amhara; this association is genome-wide significant after correction for the 985,385 tests performed ( Figure 2B and Table S3). In addition, the next six most strongly associated SNPs are in strong LD with rs10803083 (r 2 ≥0.69). There are no known genes within 600kb of these SNPs, and the closest genes, i.e. phospholipase D family member 5 (PLD5) and centrosomal protein 170kDa (CEP170), are not obvious candidate genes for variation in Hb levels. These SNPs do not reside within an ultra conserved sequence element. Notably, the effect size of SNP rs10803083 (~0.83 g/dL, Figure  2C) is half that of the EGLN1 SNPs [16] but comparable to that of Hb associated EPAS1 SNPs in Tibetans [17]. Although SNP rs10803083 reaches genome-wide significance levels, replication studies will be needed to further assess the evidence for an association with Hb levels.
Our sample size is small relative to traditional GWAS, thus liming the power of our analysis. In addition, due to the correlation between linked SNPs, the Bonferroni correction we applied is overly conservative. Therefore, many true Hb concentration variants may not reach genome-wide levels of significance. In this regard, it is interesting to note that the second strongest association signal (rs2899662) is located within the established hypoxia-candidate gene retinoid-related orphan receptor alpha (RORA). RORA encodes a protein that induces the transcriptional activation of hypoxia-induciblefactor-1alpha (HIF-1α) [26], thus it plays a significant role in the same pathway where adaptations were detected in Tibetans. This SNP is, therefore, an excellent candidate for variation in Hb levels. Additional SNPs of interest for follow up analyses (Table S3) include: the solute carrier family 30 member 9 (SLC30A9), which is regulated by hypoxia [27], the collagen type VI alpha 1 (COL6A1), which is associated with performance during endurance cycling [28] and is a HIF response gene [29], and the hepatocyte growth factor (HGF), which is induced by hypoxia [30], activates HIF1 DNA binding [31], plays a role in angiogenesis and protects against hypoxia induced cell injury [32,33].

Comparing the genetic architecture of Hb levels between Ethiopians and Tibetans
The above association analyses allow comparing the genetic architecture of Hb levels between Ethiopians and Tibetans [16,17,18].
First, we focused on the EPAS1 and EGLN1 SNPs that were previously associated with variation in Hb levels in Tibetans, with effect sizes of 0.8 g/dL and 1.7 g/dL, respectively [16,17]. None of these SNPs were significantly associated with Hb levels in Ethiopians (Table 1). Because we have complete or nearly complete power to detect a genotype-phenotype association in our Ethiopian samples (see Text  S4 and Table S21), we infer that the SNPs associated with variation in Hb levels in the Tibetans do not make a contribution in Ethiopians.
Second, because the association signal in the Tibetans may be due to an untyped variant that is tagged by different SNPs in Tibetans and Ethiopians, we also considered all SNPs within 10kb of the EPAS1 and the EGLN1 genes and repeated this analysis applying a Bonferroni correction for the number of tests performed. None of the EPAS1 or EGLN1 SNPs was significantly associated with Hb levels in Ethiopians. Based on power analyses (Figure 3), we can exclude associated variants with the same effect size as in Tibetans if the MAF in Ethiopia is greater than 10% and 5%, respectively, for the EPAS1 and EGLN1 genes (see Figure S16 for Amhara and Oromo). Therefore, this more comprehensive analysis suggests that genes shown to contribute to variation in Hb levels in Tibetans either do not influence variation in the Ethiopian populations or if they do, their effect sizes are lower than those reported for the Tibetans.
Third, we considered all SNPs within 10kb of the candidate genes in the "Response to Hypoxia" Gene Ontology (GO) category (26 genes). None of these SNPs is significantly associated with Hb levels after multiple test correction (p < 0.05/1309 = 3.81x10 -5 ). Because of the larger number of SNPs tested, this analysis has a relatively high multiple testing burden. Nonetheless, we find that we have greater than 80% power to detect a SNP significantly associated with Hb levels and effect size 0.8g/dL if its MAF is at least 20% and 100% power if its effect size is 1.7g/dL Hb ( Figure 3C and Figure S16 for Oromo and Amhara). Therefore, we conclude that variation within the "Response to Hypoxia" GO category genes is unlikely to have the same effect on Hb levels in Ethiopians as that observed in Tibetans.

Assessing genetic differences between high and low altitude populations
A widely used family of approaches for the detection of beneficial alleles uses information about the haplotype structure around the selected site [34,35]. However, these approaches have adequate power only in the case of new advantageous alleles that were driven to high frequency by natural selection, i.e. ≥70% [34,36]. Because the largest allele frequency differences observed between HA and LA among Amhara or Oromo is less than 40%, these approaches are unlikely to be powerful in this setting. Therefore, to identify alleles that contribute to genetic adaptations to HA in Ethiopians, we used two complementary approaches that focus on the divergence of allele frequency between HA and LA populations. One of these approaches was previously used to successfully identify adaptive alleles in Tibetan highlanders [18].
The first approach is based on the population branch statistic (PBS) [18], which summarizes information about the allele frequency change (PBS A_BC ) at a given locus in the history of a population (population A) since its divergence from two populations (population B and C) so that a high PBS A_BC value represents a marked change in allele frequency on the branch leading to population A. This approach was previously used to detect advantageous alleles in Tibetans relative to Han Chinese and Europeans [18]. We tested for an excess of high allele frequency differentiation (i.e. large PBS values) on the branch leading to the Ethiopian populations in SNPs within candidate genes for response to hypoxia (i.e., genes within "Response to Hypoxia" GO category) relative to SNPs in all other genes (Table S22 lists all the population trios tested). Specifically, we calculated the ratio of the proportion of SNPs in hypoxia genes vs. the proportion of SNPs in all other genes in the top 0.5%, 1% and 5% of the distribution of PBS values and used bootstrap resampling to assess the significance of the excess of large PBS values. Although an excess was observed in most population trios (Table S22), this excess was rarely statistically significant; this finding suggests that levels of linkage disequilibrium in hypoxia genes tend to be higher than in other genes and that this feature may be a confounder in tests for selection [18]. A significant excess of large PBS values in hypoxia genes was observed only in the HA Amhara and the entire Amhara sample (Table S22 and Figures S17A-B), thus suggesting that HA Amhara indeed evolved genetic adaptations to hypoxic environments. When we extended this analysis of these same population trios to additional gene classifications (i.e. BioCarta, KEGG, Gene Ontology), we found significant enrichments for SNPs in gene sets related to cell cycle control, response to DNA damage and DNA repair (Table 2).
Interestingly, however, the SNPs with the highest PBS values are found in genes with a wellestablished role in pathogen response (Table S23-S24). More specifically, the SNPs with the highest PBS values are located within the major histocompatibility complex class II DR alpha (HLA-DRA). Moreover, the null allele (FY*0) at the Duffy blood group locus, which protects against Plasmodium vivax malaria [37] and predicts white blood cell and neutrophil counts [38], has the second highest value. Consistent with expectations based on the protective effects of the FY*0 allele against malaria, its frequency is lower at HA compared to LA, where malaria is endemic (51.5% vs. 74.1%). Therefore, these results suggest that, in Ethiopian populations, differences in pathogen loads between LA and HA environments result in stronger selective pressures compared to differences in oxygen levels.
SNPs with large, even though not extreme PBS scores and lying within genes known to play an important role in hypoxia are of potential interest for follow up studies. These genes include: Cullin3 (CUL3), which potentiates HIF- In a complementary analysis, we developed a multiple regression (MR) approach to identify SNPs that show high allele frequency differentiation in HA populations relative to predictions based on a large set of worldwide population samples. This method should also be able to predict allele frequencies appropriately in a situation where the target population is admixed. In this approach, we used allele frequency data from 61 LA populations (including the HGDP and several other populations) to predict the expected allele frequencies in the HA Amhara. We focused on the HA Amhara because they have lived at HA for a longer period of time and exhibit distinct patterns of Hb and O 2 sat levels compared to the Oromo (Figure 1). In addition, we omitted the LA Amhara in an attempt to reduce the effect of gene flow between altitudes, which could potentially reduce our power to detect adaptive divergence. We used all SNPs to estimate the best-fitting regression coefficients for each population: that is, these are the coefficients that generate the lowest mean square error in predicting the HA Amhara allele frequencies. The populations with the largest regression coefficients in the Amhara regression model are from geographically proximate populations in East Africa (Maasai, Luhya and LA Oromo) and from the Middle East and Southern Europe (see Figure S18). We reasoned that changes in allele frequencies due to high altitude adaptation would be detectable as departures (i.e. large residuals) from the predicted allele frequencies based on all other populations As for the PBS analysis, we tested for an excess of SNPs with high allele frequency differentiation using the MR statistic for genes within "Response to Hypoxia" GO category relative to SNPs in all other genes and we used a bootstrap procedure to assess the significance of the observed excess. An excess was observed for all tail cut-offs, but only one reached statistical significance ( Table  2). Other gene sets that showed a significant enrichment of SNPs with strong MR signals include chromosome organization and biogenesis, DNA repair, histone modification. These findings are consistent with the pattern observed in the PBS analysis, indicating that they are robust to the choice of populations used in the test (Table 2).
Among the SNPs with the largest MR scores, there are several SNPs in hypoxia genes, which may be of potential interest for follow up studies (Table S27). The SNP (rs12510722), which shows the 4 th highest MR scores, lies within the alcohol dehydrogenase 6 (ADH6) gene, whose expression is affected by pseudohypoxia [47]. The SNPs with the 8 th and the 15 th highest MR score (rs2660342 and rs2660343, respectively) are within 100kb from solute carrier family 30 member 9 (SLC30A9) and transmembrane protein 33 (TMEM33). SLC30A9 is up-regulated by hypoxia [27] while TMEM33 is down-regulated under hypoxia and up-regulated after knockdown of HIF1A [41].

Assessing epigenetic differences between high and low altitude populations
Methylation is an epigenetic modification that is known to play a crucial role in the cellular response to hypoxia [48]. Since HA adaptation could be in part maintained by methylation, we measured methylation levels at 27,578 CpG sites in 17 HA and 17 LA Amhara and 17 HA and 17 LA Oromo. CpG methylation levels were tested in DNA extracted from blood in the Amhara and from saliva in the Oromo. To avoid confounding due to differences in methylation across tissues, we performed the comparison across altitudes within each ethnic group.
In Oromo, four CpG sites reached significance after multiple test correction (p <1.85x10 -6 ), but the closest genes are not known hypoxia candidate genes: apolipoprotein B mRNA editing enzyme catalytic polypeptide-like 3G (APOBEC3G), metallothionein 1G (MT1G), paired-like homeodomain 2 (PITX2) and olfactory receptor family 2 subfamily K member 2 (OR2K2) (Table S28). Interestingly, APOBEC3G codes for a well-established cellular antiviral protein and a specific inhibitor of human immunodeficiency virus-1 (HIV-1) infectivity [49,50]. The MT1G gene also has a role in HIV-1 infection because it upregulates MT1G expression in immature dendritic cells, which in turn facilitates the expansion of HIV-1 infection [51]. Although the prevalence of HIV was not surveyed in our fieldwork, HIV/AIDS is known to be a major health problem in Ethiopia [52,53]. While these functions for APOBEC3G and MT1G point to a role for methylation in defense against pathogens, MT1G also plays a role in the response to hypoxia as its promoter is induced by vascular endothelial growth factor (VEGF), which in turn contributes to the prosurvival and angiogenic functions of VEGF [54]. Likewise, expression of PITX2 is required for normal hematopoiesis [55,56], raising the interesting scenario that methylation of this gene may influence beneficial phenotypes in the response to hypoxia. Some of the CpG sites with nominally significant (p<6.7x10 -5 ) differences in methylation between HA and LA are close to genes that are differentially expressed in response to hypoxia; these genes include: toll-like In Amhara, no CpG site showed a methylation difference that reached significance after multiple test correction (Table S29). However, we note that the 3 rd most significant differentially methylated CpG site (p = 8.03x10 -05 ) was closest to Glutathione-S-Transferase (GSTP1) whose expression is increased by prolonged hypoxia [59] and whose loss of expression correlates with methylation in prostate cancer [60]. Hypoxia also regulates the expression of genes close to other differentially methylated CpG sites in Amhara (p<1. Finally, no significant excess of methylation differences between LA and HA populations was observed at the genome-wide level in Oromo or Amhara ( Figure S19) nor did we find a significant enrichment of methylation differences between LA and HA populations for gene sets defined by BioCarta or KEGG pathways and by Gene Ontology categories (data not shown).

DISCUSSION
HA human populations across the world allow studying independent realizations of the adaptive process in response to the same selective pressure, i.e. hypoxia, thus providing an excellent opportunity to investigate how natural selection shapes the genetic architecture of adaptive traits. To make progress on these enduring questions, we have sampled two closely related ethnic groups in the Ethiopian highlands that include both HA and LA residents, thus allowing comparisons across altitudes within and between ethnic groups. Of these two groups, the Oromo have moved to HA only 500 years ago [61,62], thus making it unlikely that genetic adaptations evolved in this group. In contrast, the Amhara have a history of HA residence of at least 5ky and possibly as far as 70ky [24,25]. Because previously identified selection signals [63,64,65,66] occurred within a similar period of time, including HA adaptations [18], we conclude that enough time has elapsed since the Amhara moved to HA for genetic adaptations to have taken place. Consistent with this idea, we observe significant phenotypic differences between Amhara highlanders and the more recent HA residents, i.e. the Oromo. While HA Amhara are characterized by mildly elevated Hb levels (similar to Tibetans) and no or mildly reduced O 2 sat [13], the HA Oromo sample resembles acclimatized lowlanders with a response characterized by elevated Hb concentration and marked reduction in O 2 sat. Our data indicates that Amhara and Oromo are very similar at the genome-wide level, therefore, the observed phenotypic differences are likely due to the different histories of HA occupation. In addition to these phenotypic comparisons, our genomic analyses of these two ethnic groups resulted in several important observations that shed new light on the biology of HA adaptations and that are discussed in detail below.
First, in a GWAS of Hb levels in Amhara, we find a genome-wide significant signal of association as well as an excess of low p-values. In addition, the second most strongly associated SNP is found within the RORA gene, which belongs to HIF1 pathway and is, therefore, an excellent candidate gene for hypoxia response phenotypes. Additional strong associations were observed in other candidate hypoxia genes, such as COL6A1, SLC30A9, and HGF. We looked at the Amhara data of Scheinfeldt et al [21] to test for replication of the association signal at these genes. Two of them showed suggestive associations with Hb levels (p=0.06 and p=0.15 for SLC30A9 and RORA, respectively), with β values (-0.60 and 1.31, respectively) consistent with ours (-0.67 and 0.92, respectively). It should be noted that the replication test was performed in only 21 Amhara samples in Scheinfeldt et al [21] for which age and BMI data were available and who had Hb levels within the normal range; thus, the lack of replication may well be due to the very low power of the replication sample. We note that, though variation in EPAS1 and EGLN1 has been consistently associated with Hb levels in Tibetans, no genomewide significant association signal and no excess of low p-values were observed (see Figure S5 in Simonson et al [16,17]). While the signals we detected await replication, it is interesting to note that their effect sizes are as high as those found in Tibetans for SNPs in EPAS1, thus raising the interesting scenario that selection may have favored alleles with similar effect sizes on Hb levels even though the specific loci contributing to the trait are different.
Some interesting patterns are beginning to emerge with regard to the genetic contribution to variation in Hb levels and O 2 sat, the two phenotypes that have been most widely studied in highlander populations. No evidence of a genetic contribution to O 2 sat in Amhara, Oromo, and the combined Ethiopian sample could be detected ( Figures S7-S15). This is true also in the Tibetans, even though segregation analysis detected a major O 2 sat locus, which is also associated with reproductive success [15]. Therefore, the data so far suggest that while genetic factors contribute to variation in Hb levels, their importance in O 2 sat is lower. This is consistent with studies in Tibetans and Andeans showing a markedly lower heritability for O 2 sat compared to Hb levels; indeed, the O 2 sat heritability in Andeans was not significantly different from zero [67,68,69]. More data, especially at the genome-wide level, are needed to elucidate the contribution of genetic factors to these two phenotypes.
Second, to compare the genetic bases of Hb variation with the Tibetans, we tested for an association between SNP genotypes and Hb levels within the Ethiopians. Although we had appropriate power, none of the SNPs within 10 kb of the EPAS1 and EGLN1 genes or the genes in the hypoxia pathway, including SNPs previously associated with Hb variation (and signatures of natural selection) in Tibetans, associated with Hb. Therefore, we can rule out that the SNPs and loci contributing to Hb variation and showing selection signals in Tibetans affect the same trait in Ethiopians even though Tibetans and Amhara have lower Hb levels compared to all other highlanders. Alternatively, if the same variant affects Hb levels in both populations, their effect sizes in Amhara must be markedly lower than those reported for the Tibetans.
Third, by using approaches based on allele frequency divergence, we find that outlier SNPs in the HA Amhara sample are not in hypoxia response genes, but in loci known to play a role in immune defense. These include variation in the HLA-DRA locus and the null allele at the Duffy blood group locus; none of these variants is associated with Hb levels or O 2 sat. Interestingly, malaria and schistosomiasis were prevalent in the LA, but not in the HA Amhara communities sampled in this study (Text S1), reflecting important differences in pathogens between HA vs. LA environments. Indeed, epidemiological studies in the areas near the LA sampling sites for the Amhara and Oromo reported malaria prevalence of 39.6% and as much as 25% of malaria morbidity is due to P. vivax [70,71,72]. Therefore, the immune defense variants with extreme frequency divergence represent excellent candidates as selection targets. These findings are important in several respects. First, they indicate that, because HA and LA habitats differ by multiple environmental stresses, signals of allele frequency divergence cannot be unambiguously attributed to hypoxia without additional information about the gene function or the specific phenotypic effects of the alleles. Second, they further corroborate that the Amhara populations are indeed adapted to spatially-varying selective pressures, despite likely high levels of gene flow between HA and LA communities. Third, the fact that variants in hypoxia response genes are not outliers in these analyses suggests that adaptations to pathogens and to hypoxia have a different genetic architecture or that the intensity of pathogen-related selective pressures is stronger than those due to hypoxia. Overall, these findings highlight the opportunities and challenges of ecological genomic studies and point to the power of approaches that use environmental information combined with phenotypic data collected in the field.
Although selection has not created extreme HA vs. LA frequency shifts in hypoxia genes, we find the SNPs within candidate genes for response to hypoxia as a group show an excess of allele frequency differentiation based on the PBS analysis performed in the Amhara (Table 2). This approach, however, requires specification of a set of three populations to be tested. When we used our MR approach, which uses information from all populations simultaneously, the response to hypoxia genes did not show a significant excess ( Table 2). In contrast, both approaches detected a significant excess of allele frequency divergence in genes involved in cell cycle control, DNA repair, DNA damage, chromatin structure and modification, consistent with the known role of oxygen sensing in the regulation of cell proliferation [73]. Therefore, our findings raise the possibility that genetic variation in these pathways can contribute to adaptations to HA.
Fourth, though we did not observe a genome-wide excess of methylation differences between HA and LA samples, we found genome-wide significant signals in the Oromo in genes with a known function in pathogen defense or in the biology of hypoxia, i.e. VEGF signaling and hematopoiesis. Interestingly, no genome-wide significant differential methylation was observed in the Amhara; this may be due to the fact that DNA from different tissues was analyzed in the two ethnic groups. However, given the difference in the history of HA residence between Amhara and Oromo, it is also possible to speculate that epigenetic modifications play a role in the early phases of adaptations to new environments and that this role is replaced over time by adaptations at the genetic level.
It has been proposed that epigenetic modifications are important in ecological adaptations [74] and methylation is known to play a crucial role in the cellular response to hypoxia [48]. A previous study showed that gene expression differences observed between Moroccan populations were not explained by methylation differences on the tested 1,505 CpG sites [75]. Our study is more comprehensive having interrogated 27,578 CpG sites in a larger sample size. Though our results do not unambiguously point to a major role for methylation in the adaptations to HA, they suggest that further studies using DNA extracted from different tissues and including additional epigenetic modifications, in addition to methylation, are warranted.
In conclusion, studies of genetic variation in indigenous populations with long-time residence at HA are giving rise to a composite picture regarding the genes contributing to and the genetic architecture of HA adaptations. Clearly, different loci contribute to Hb levels in Ethiopians and Tibetans. Moreover, while some aspects (i.e. phenotypic effect sizes) of the genetic architecture of Hb levels may be similar, others (i.e. the allele frequency shifts due to selection) are different. Additional examples of environmental pressures that acted on different human populations include malarial endemia, low UVB radiation levels, and an adult diet rich in dairy products. Detailed genome-wide studies of parallel adaptations to these selective pressures are needed to elucidate the impact of natural selection on the genetic architecture of complex adaptive traits.

MATERIALS AND METHODS Ethics statement
All participants in the study gave informed consent. The studies were approved by the Institutional Review Boards of Case Western Reserve University and of the University of Chicago and by the Ethiopian Science and Technology Council Ethics Review Board.

Population samples
Samples were drawn from native residents of the Semien Mountains area of northern Ethiopia inhabited mainly by Christian Amhara and the Bale Mountains area of southern Ethiopia inhabited mainly by Muslim Oromo (also referred to as Galla, Boran and Gabbra). For each ethnic group, we sampled individuals at HA and LA. The LA samples were chosen to achieve the maximum altitude contrast, yet avoid confounding due to the presence of endemic malaria; all LA samples were from agropastoral people reporting the same ethnicity as the HA samples and no visits to altitudes above 2500 m in the past 6 months. For details on the sampled populations and their ecology, see Text S1.
DNA was extracted from blood samples provided by 192 Amhara individuals living at 3700 m in the Simien Mountains National Park or at 1200 m in the town of Zarima. Forty-seven of these individuals were sampled in 1995 and previously described [13,76], the remaining 145 individuals were sampled in a separate expedition in 2005. DNA was extracted from, saliva samples provided by 118 Oromo individuals and collected using Oragene DNA sample collection kits; 79 individuals lived at 4000 m in the Bale Mountains National Park while 39 individuals lived at 1560 m in the town of Melkibuta. The study participants were healthy (refer to Text S3 for details). Hb and O 2 sat of Hb were measured in all individuals. Hb was determined in duplicate using the cyanmethemoglobin technique (Hemocue Hemoglobinometer, Hemocue AB, Angelholm, Sweden), immediately after drawing a venous blood sample. O 2 sat was determined by pulse oximetry (Criticare Models 503 and SpO2) as the average of six readings taken 10 seconds apart.

SNP Genotyping and imputation
The samples were genotyped using Hap650Yv3 (n=46), Human1M-duoV3 (n=112), and Human Omni-Quad1 (n=160) Illumina arrays at Southern California Genotyping Consortium. Nineteen samples with less than 93% genotype call rate were omitted from the analysis. We used the program RELPAIR 2.0.
[77] to test for hidden relatedness in the samples by using 3 independent sets of 1800 autosomal and 200 X-linked SNPs; individuals with relationships closer than first cousins to any other individual in the sample were omitted from the analysis. After applying these filters, 260 unrelated samples remained: 102 HA Amhara, 60 LA Amhara, 63 HA Oromo and 35 LA Oromo. Principal component analysis (PCA) of the genotype data did not detect any major outlier in either population (see Text S2 and Figure  S3. Due to the incomplete overlap between SNPs on the three genotyping arrays, genotypes for 1,819,369 HapMap3 SNPs were imputed by using the program IMPUTE2 [78] and the HapMap populations as a reference panel: Utah residents with Northern and Western European ancestry from the CEPH collection (CEU), Han Chinese in Beijing, China (CHB), Japanese in Tokyo, Japan (JPT), Luhya in Webuye, Kenya (LWK), Maasai in Kinyawa, Kenya (MKK), Toscani in Italy (TSI) and Yoruba in Ibadan, Nigeria (YRI) samples. A total of 1,297,134 autosomal SNPs with minor allele frequency (MAF) higher than 0 and imputation accuracy higher than 90% were used in the downstream analyses ( Figure S20).

Genotype-phenotype association
Phenotype-genotype associations were tested at each SNP with MAF>0.1 by linear regression using the whole-genome association analysis toolset in PLINK [79]. To determine whether there was an excess of low p-values taking into account the large number of tests performed, the distribution of observed p-values from the linear regression tests were compared to a null distribution obtained by permuting 100 times the phenotype (and corresponding covariate variablessee below) value across individuals and running the same linear model. The results were visualized by means of quantilequantile (QQ) plots and the 95% confidence interval (CI) was estimated by permutations. Gender, body mass index (BMI), altitude, ethnic group and year of sampling were used as covariates in the linear regression. As an alternative approach, we grouped the samples based on their gender, altitude, ethnic group and year of sampling and we quantile normalized the observed phenotypes within each group; these phenotypes were then pooled across groups before testing for an association with genotype by linear regression. This latter approach preserves information about the individual ranks without introducing a bias due to the effects of covariates. However, it does not preserve information about the phenotype values, thus leading to a reduction in power. The results of these two approaches correlated highly (r 2 > 0.92); therefore, only the results for the linear regression with the covariates are shown.

Population genetics analyses
The population branch statistic (PBS A_BC ) was used to summarize the amount of allele frequency change in the history of population A since its divergence from two related populations, B and C [18]. In a complementary approach based on multiple regression (MR), we used the allele frequencies for the LA Oromo and world-wide population samples (i.e., Human Genome Diversity Project (HGDP) panel populations [80] and 4 HapMap Phase III populations -LWK, MKK, TSI and Gujarati Indians in Houston, Texas (GIH) -(www.hapmap.org)) to detect SNPs whose HA Amhara frequencies deviate most from estimated frequencies. Briefly, denoting the observed allele frequencies in the p populations by X 1 , X 2 , …, X p and the expected allele frequency within the HA Amhara by y, we can predict the allele frequency for the i th SNP using the linear model y i = β 0 + β 1 X i1 + β 2 X i2 + … +β p X ip + ε i . Using the data from all genotyped SNPs that overlap among the p populations and applying multivariate linear regression, we found b 0 , b 1 , b 2 , …, b p , which are estimates of parameters β 0 , β 1 , β 2 , …, β p and estimated the magnitude of the residual for the i th SNP as |ε i |= |y i -(b 0 + b 1 X 1 + b 2 X 2 + … + b p X p )|. We refer to this residual as the MR score.
To test for an enrichment of allele frequency differentiation in candidate hypoxia genes, we focused on the 28 genes belonging to the "Response to hypoxia" category in the Gene Ontology database (GO:0001666). We then compared the proportion of SNPs within 10 kb of each gene in this category to that of SNPs within 10 kb of all other genes in the tail of the PBS and MR score distributions. Given the arbitrary nature of choosing a single cutoff for the tail of the distribution, we set three cutoffs (5%, 1% and 0.5%). In other words, we looked at the top 5%, 1% and 0.5% of all PBS and MR values and asked whether there is an enrichment of hypoxia gene SNPs relative to all other genic SNPs for each tail cut-off. A value of 1 represents no excess and a value greater than 1 represents enrichment in the tail of the distribution. SNPs are likely to cluster along the genome due to linkage disequilibrium, thus reducing the number of independent signals contributing to an observed enrichment. To account for this possibility, we found the confidence interval for the enrichment using a bootstrap approach described in Hancock et al [81]. An enrichment of hypoxia SNPs was considered significant (with a one-tailed test) if at least 95% of the 1000 bootstrap replicates were enriched (i.e., had a ratio above 1).

Methylation
Methylation levels were measured at 27,578 CpG sites in 17 HA and 17 LA Amhara and Oromo DNA samples, for a total of 68 individuals, using six Infinium HumanMethylation27 arrays at Southern California Genotyping Consortium. Two LA Amhara sample data were discarded due to low data quality. In each ethnic group, the HA methylation level of each CpG site was compared to corresponding LA levels using the following linear model correcting for age, gender, and the methylation array: lm (% methylation ~ altitude + gender + age + [array]). Two comparisons were performed: (1) HA vs. LA Amhara and (2) HA vs. LA Oromo. As for the genotype-phenotype association, excess of differential methylation between HA and LA was tested by comparing the observed p-value distribution to the null distribution obtained by permuting 100 times the altitude label across individuals and running the same linear model. Permutation based 95% CI was estimated.

Text S1. Sampled populations and their ecology
The Simien Plateau in Northwest Ethiopia is cut by streams and rivers and has cliffs and land areas sloping down to canyons. The high altitude (HA) Amhara are agropastoralists living in a temperate Afro-alpine ecosystem in the Simien Mountains National Park at altitudes ranging from 3500-4100 meters (m). Altitudes above 2500m on the East African Plateau have been inhabited for at least 5 thousand years (ky) and altitudes around 2300-2400m for more than 70ky [24,25]. Linking historical and modern ethnic groups in Ethiopia with prehistoric sites is not feasible with current knowledge. This paper considers that 5ky is a reasonable lower estimate for human habitation of this area. The Bale Plateau in Southeast Ethiopia is a flat mesa with sharp relief down to the lowlands. The HA Oromo are pastoralists herding cattle, sheep and goats and living in a temperate Afro-alpine ecosystem in the Bale Mountains National park and reside on the Sanetti Plateau at 4000-4100m. The HA areas of the Bale Plateau have been inhabited by Oromo since the early 1500s according to historical records [22,23]. The two plateaus are about 500 miles apart and the intervening terrain is not a continuous plateau as the Tibetan or Andean, instead it is a mosaic of terrain generally above 1500m and punctuated by lowlands The Amhara reported occupations as farmers and, at LA, farmers and traders. The HA sample reported eating a diet based on barley injera (a flat crepe-like bread made of fermented dough) with sauce (wot) of potato or beans. The LA Amhara community was exposed to malaria and schistosomiasis and had a high prevalence of visible goiter and iron deficiency. Such individuals were excluded from the sample. The HA Oromo sample reported eating a diet of teff injera with potato or bean wot and some meat and dairy products. The LA Oromo were farmers, traders, and government officials such as teachers; they were also exposed to malaria.

Text S2. The genetic structure of the Amhara and Oromo populations
In order to investigate the genetic relationship between Amhara and Oromo and other populations, we used genotype data in the two Ethiopian populations for SNPs that had previously been typed also in a set To look more closely at the relationship between the Amhara and the Oromo, we performed principal component analysis (PCA) using the genotype data from 13,000 random autosomal SNPs for all Ethiopian individuals. In this analysis, PC1, which explains 0.78% of the variance, separates the majority of the Amhara from the majority of the Oromo individuals; in addition, within each ethnic group, the HA and LA individuals tend to cluster separately ( Figure S3. Accordingly, the mean F ST between LA Amhara and Oromo is very low (0.0098), though it is higher than the mean F ST between LA and HA populations within the same ethnic group (HA vs. LA Amhara: F ST =0.0047; HA vs. LA Oromo high F ST =0.0074). The STRUCTURE analysis did not reveal detectable differences between the HA and LA populations (10 6 iterations; following a burn-in period of 30,000 iterations; k=3) ( Figure S4. However, when we compared the F ST distribution between HA and LA populations (Amhara or Oromo) to the distribution obtained by permuting the altitude labels (i.e. under random mating for HA and LA populations within each ethnic group) we detected a significant excess of high F ST values, indicating that some population structure between altitudes does exist ( Figure S5). In order to test for an excess of allele frequency divergence between populations, the observed distribution of F ST was compared to a null distribution obtained by permuting the population labels 100 times, and for each comparison, a different permuted F ST distribution was created. Permutations were used to estimate the 95% confidence interval.

Text S3. Phenotypic variation in Amhara and Oromo
260 unrelated people provided the phenotype and genotype data presented here. HA Amhara were sampled in 1995 (n=47) and 2005 (n=65) and the LA Amhara in 2005-6 (n=59). All Oromo were sampled in 2007 (n=88). These unrelated people were identified from a larger sample of partially related individuals.
In the Semien Mountain area, two models of Criticare oximeter were used, the Criticare Model 503 was used for 152 people and the SpO2 was used for 150 people. Twenty-five people were measured using both instruments and the average difference (SpO2 minus 503) was =0.031 + 2 (SD) %. One observation was just outside the limits of agreement of the two models (+ 2 SDs). The Model 503 values were used as the O 2 sat measured by Criticare oximeter for the 127 who did not have SpO2 measurements. Table S1 summarizes height, weight, body mass index (BMI), and pulse. Considering altitude differences, HA Amhara native men and women were shorter, lighter and had lower BMI than their lowland Amhara counterparts. HA Oromo native men had higher BMI than their lowland Oromo counterparts although there were no significant differences in height or weight. HA Oromo native women were shorter but did not differ in weight or BMI from their lowland counterparts. Considering ethnic differences within altitude zones, HA Amhara men and women were shorter, lighter and had lower BMI than Oromo. LA Amhara men were heavier and had higher BMIs. Amhara lowland women did not differ from Oromo lowland women in these characteristics although there was a trend toward shorter height. Table S2 summarizes the three HA adaptation phenotypes, hemoglobin (Hb) concentration, percent of oxygen saturation (O 2 sat) of Hb and calculated arterial oxygen content (AOC). Arterial oxygen content in ml O 2 /dL was calculated as 1.39*(Hb*O 2 sat)/100 [86]. Phenotypes were assessed using the same equipment and protocols in all areas [13]. To remove many possible confounding factors that could have added spurious variation to the phenotypes, the analyses were confined to healthy people (based on self-report and a review of systems by an Ethiopian physician), who were free of respiratory symptoms, not hypertensive, not pregnant, had not delivered an infant in the past year and had not visited at another altitude (below 2500m for highlanders or above 2500m for lowlanders) in the past six months. Samples collected in 2005, 2006 and 2007 were also screened for normal lung function and for infection with all four species of human-associated malaria [87]. Furthermore, because poor iron status could limit the Hb response to HA hypoxia, the analyses were limited to those with normal iron status. Iron deficiency in the 1995 sample was identified on the bases of zinc erythrocyte protoporphyrin > 70 gm/dL, plasma ferritin concentration <12 ng/ml or transferrin receptor concentration <8.3 g/L [13]. Iron deficiency in the samples collected later was identified on the bases of body iron stores calculated using the log of the ratio of transferrin receptor to ferritin concentration [88]. If an individual were missing data for any of these variables, then his or her Hb concentration data were not analyzed for association with genetic variants. Stringent criteria for phenotypes improve the likelihood of finding phenotypic variation associated with health, normal genetic variation and add confidence to findings of a lack of association.
With respect to altitude differences, Hb was significantly higher among all four HA samples compared with the age-sex groups at LA but size of the effect was smaller among the Amhara (Table S2 and Figure S6. Similarly, the percent of O 2 sat was lower among all four HA samples compared with the age-sex groups at LA, but the size of the effect was smaller among the Amhara (Table S2 and Figure S6. The altitude sub-samples did not differ from one another in calculated AOC. With respect to ethnic group differences within altitude zones, Amhara had lower Hb concentrations at both HA and LA (Table  S2 and Figure S6 show that the differences were more than one gm/dL of Hb and more than one standard deviation). In contrast, Amhara had higher percent of O 2 sat at HA than their Oromo counterparts. The result was similar for AOC among highlanders but a trend toward lower AOC among LA Amhara males that was statistically significant among Amhara females was observed.
In summary, there were altitude and ethnic group differences in phenotypes with the Oromo showing a larger altitude response than the Amhara.

Text S4. Comparing the genetic architecture of Hb levels between Ethiopian and Tibetans
To compare the genetic architecture of Hb levels between Ethiopians and Tibetans [16,17,18], we performed association tests and power analyses at three levels: 1) the SNPs in EPAS1 and EGLN1 associated with the trait in Tibetans, 2) all the SNPs within the same genes, and 3) all the SNPs within all the genes in the same pathway, i.e. response to hypoxia.
The EPAS1 and EGLN1 SNPs that were previously associated with variation in Hb levels in Tibetans have effect sizes of 0.8 g/dL and 1.7 g/dL, respectively [16,17]. Of these SNPs, we considered those that were genotyped or imputed with greater than 90% accuracy in our Ethiopian data. Because the lack of association in the Ethiopians could be due to incomplete power, we calculated the probability of observing a significant association between SNP genotype and Hb levels in our data, assuming that the  coefficient for each SNP is as high as or higher than that observed in Tibetans and using the corresponding MAF in the Ethiopians. As shown in Table 1, we have complete or nearly complete power to detect a genotype-phenotype association in our Ethiopian samples. This suggests that either the SNPs associated with variation in Hb levels in the Tibetans do not make a contribution in Ethiopians or if they do their phenotypic effect is smaller in Ethiopia. However, given their corresponding MAFs in Amhara, Oromo and Ethiopian, we have nearly complete power to detect a significant association even if the phenotypic effect was half as large as that observed in the Tibetan studies (Table S21).
Next, we considered all SNPs within 10kb of the EPAS1 gene and repeated this analysis applying a Bonferroni correction for the number of tests performed (i.e., 72 SNPs). We find that we have greater than 95% power to detect a SNP with an effect size of at least 0.8 g/dL Hb concentration if the MAF in Ethiopia ( Figure 3A) is greater than 10% (see Figure S16 for Amhara and Oromo). When we performed the same analysis for SNPs within 10kb of the EGLN1 gene (38 SNPs), we find that we had 100% power to detect SNPs with effect size of at least 1.7 g/dL Hb concentration if the SNP MAF in Ethiopia is greater than 5% ( Figure 3B and Figure S16 for Oromo and Amhara). Because none of the EPAS1 or EGLN1 SNPs was significantly associated with Hb levels in Ethiopians after multiple test correction, this analysis suggests that genes shown to contribute to variation in Hb levels in Tibetans either do not influence variation in the Ethiopian populations or if they do, their effect sizes are lower than those reported for the Tibetans.
Finally, we considered the SNPs within 10kb of the candidate genes in the "Response to Hypoxia" GO category (26 genes). Because of the larger number of SNPs tested, this analysis has a relatively high multiple testing burden. Nonetheless, we find that we have greater than 80% power to detect a SNP significantly associated with Hb levels and effect size 0.8g/dL if its MAF is at least 20% and 100% power if its effect size is 1.7g/dL Hb ( Figure 3C and Figure S16 for Oromo and Amhara). Figure S1.                     P < 0.01 t-test comparing same sex and ethnic group at high and low altitude + 0.05 < p < 0.01 t-test comparing same sex and ethnic group at high and low altitude a p < 0.05 t-test comparing Amhara and Oromo of the same sex at one altitude b P < 0.01 t-test comparing Amhara and Oromo of the same sex and ethnic group at one altitude c 0.05 < p < 0.01 t-test comparing Amhara and Oromo of the same sex and ethnic group at one altitude 2 + 0.50 *p < 0.05 t-test comparing same sex and ethnic group at high and low altitude **P < 0.01 t-test comparing same sex and ethnic group at high and low altitude + 0.05 < p < 0.01 t-test comparing same sex and ethnic group at high and low altitude a p < 0.05 t-test comparing Amhara and Oromo of the same sex at one altitude b P < 0.01 t-test comparing Amhara and Oromo of the same sex and ethnic group at one altitude c 0.05 < p < 0.01 t-test comparing Amhara and Oromo of the same sex and ethnic group at one altitude COL6A2,FTCD,C21orf 56,COL6A1 Only SNPs with MAF <10% and imputation accuracy > 0.9 were tested. In addition to age, sex, BMI (body mass index) and altitude, collection year was also used as covariate since samples were collected ten years apart.       Only SNPs with MAF <10% and imputation accuracy > 0.9 were tested. Age, sex, BMI (body mass index), collection year and ethnicity were used as covariates.     [16] while those for EPAS1 were obtained from Beall et al [17]. 2  indicates the observed linear coefficient for the relationship between SNP genotype and Hb levels. 3 half  indicates half of the observed linear coefficient for the relationship between SNP genotype and Hb levels. 3 Power refers to the probability of detecting a significant association (p<0.05) between SNP genotype and Hb level given the MAF and the sample size in the Ethiopian populations assuming that the  coefficient is higher than half of the observed in Tibetan  Table S22. Proportions of Response to Hypoxia genic SNPs relative to the proportion of all other genic SNPs in the tails of the PBS distributions for each population trio tested. * and ** denote support from ≥ 95% and 99% of bootstrap replicates, respectively.