The Genetic Architecture of Adaptations to High Altitude in Ethiopia

Summary

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 and 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

Hypoxia is a major stress on human physiological processes and a powerful homeostasis system has evolved in animals to cope with fluctuations in oxygen concentration [1]. High altitude (HA) hypoxia, such as that experienced at 2500 m of altitude or greater, engages this system and elicits physiological acclimatization when lowlanders become exposed to hypoxia. In addition to lower oxygen levels, lower biodiversity and extreme day-to-night temperature oscillations challenge HA living. The classic response is an increase in hemoglobin (Hb) concentration that during acclimatization compensates for the unavoidable lowered percent of oxygen saturation (O₂ sat) of Hb due to ambient hypoxia. Acclimatization is not completely effective, however. For example, birth weights are lower than at low altitude (LA) as is physical exercise performance [2], [3]. In addition, lowlanders residing at altitudes higher than 2500 meters (m) are at risk for chronic health problems arising in part from acclimatization processes. For example, long-term high Hb levels increase blood viscosity as well as the risk of thrombosis and stroke [4], [5], [6] and poorer pregnancy outcomes [7]. These results taken together suggest that the acclimatization response does not assure that fitness is unaltered at HA. Several distantly related human populations have survived for 5–50 thousand years (ky) [8], [9], [10] at altitudes above 2500 m. Indeed, in most cases, sufficient time has elapsed since HA settlement for natural selection to have changed the frequency of adaptive alleles. Interestingly, even though all HA residents are exposed to the same, constant, ambient hypoxia, indigenous highlander populations show distinctive physiological characteristics thought to offset HA stress: Andeans show some reduction in O₂ sat, but a marked increase in Hb levels [11], Tibetans present markedly low O₂ sat, but relatively little increase in Hb levels [12], and Amhara in Ethiopia present little reduction in O₂ sat or increase in Hb levels [13], [14]. Whether these phenotypic contrasts reflect different genetic adaptations across populations remains an open question.

Substantial evidence in Tibetan highlanders suggests that variation in Hb levels and O₂ sat is adaptive. In the Tibetan population, a major gene effect on O₂ sat was detected, with the inferred genotypes associated with higher O₂ sat also associated with higher reproductive success [15]; though the locus underlying this effect has not yet been identified, its effect on phenotypes directly related to fitness points to the presence of adaptive variation. With regard to Hb levels, which are surprisingly similar between Tibetan highlanders and lowlanders at sea level [11], population genetics and genotype-phenotype association analyses have identified alleles at two loci (endothelial PAS domain protein 1 (EPAS1) and egl nine homolog 1 (EGLN1)) that are consistently associated with signatures of positive natural selection and with lower Hb levels, suggesting that natural selection in Tibet favored variants that counteract the deleterious effects of long-term acclimatization [16], [17], [18], [19].

An analysis of genome-wide genotype data in Tibetan and Andean highlanders suggested that natural selection acted on largely distinct loci in the two populations [20]. In addition, a recent study comparing Ethiopian Amhara highlanders with other ethnic groups at LA identified yet another set of candidate targets of selection [21]. However, the Tibetans remain unique with regard to the strength of the evidence for natural selection and the marked genetic effects on the Hb level phenotype.

HA populations offer a rare opportunity to investigate the impact of natural selection on the genetic architecture of adaptation because independent realizations of the adaptive process can be 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.

Results

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, S2, S3, S4, 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 1500 s [22], [23], the Amhara have inhabited altitudes above 2500 m for at least 5 ky and altitudes around 2300–2400 m for more than 70 ky [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₂ sat, HA Amhara had a 5.6% lower O₂ sat compared to LA Amhara while HA Oromo had 10.5% lower O₂ 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 genome-wide 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 and S2 and Figure S6.

Genome-wide association signals in Ethiopia

To learn about the genetic bases of variation in Hb level and O₂ sat in Ethiopia, we tested for an association between SNP genotype in the 260 unrelated Ethiopian samples and Hb levels or O₂ 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, S8, S9, S10, S11, S12, S13, S14, S15). No genome-wide significant signal was observed for either Hb levels or O₂ sat in the Oromo and in the total Ethiopian sample and for O₂ sat in the Amhara (Figures S7, S8, S9, S10, S11, S12, S13, S14, S15 and Tables S3, S4, S5, S6, S7, S8, S9, S10, S11, S12, S13, S14, S15, S16, S17, S18, S19, S20). Likewise, no excess of low p-values was observed in these association analyses relative to null expectations obtained by permutations (Figures S7, S8, S9, S10, S11, S12, S13, S14, 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.96×10⁻⁸) 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²≥0.69). There are no known genes within 600 kb of these SNPs, and the closest genes, i.e. phospholipase D family member 5 (PLD5) and centrosomal protein 170 kDa (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-inducible-factor-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 10 kb 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 10 kb 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.81×10⁻⁵). 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.8 g/dL if its MAF is at least 20% and 100% power if its effect size is 1.7 g/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 versus 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 Figure S17A and S17B), 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 well-established role in pathogen response (Tables S23 and 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% versus 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-1 signaling [39], as well as adrenergic beta receptor kinase 1 (ADRBK1) [40], coronin actin binding protein 1B (CORO1B) [41], anti-silencing function 1 homolog B (ASF1B) [42] and MAPK-activated protein kinase MK2 (MAPKAPK2) [43], which are all down-regulated under hypoxia (Tables S23 and S24). None of those large PBS SNPs were significantly associated with Hb or O₂ sat (Tables S25 and S26), but a SNP within utrophin A (UTRN) - rs7753021 - reached nominal levels of significance with O₂ sat (p = 0.005; Table S26). UTRN expression correlates with oxidative capacity [44] and increases with chronic physical training [45]. Slow-twitch muscles, which are associated with endurance performance, have high levels of UTRN [46].

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₂ 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 100 kb 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.85×10⁻⁶), 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.7×10⁻⁵) differences in methylation between HA and LA are close to genes that are differentially expressed in response to hypoxia; these genes include: toll-like receptor 6 (TLR6) [57], mif two 3 homolog 1 (SUMO1) [27]; phosphodiesterase 4A (PDE4A) [58] and human immunodeficiency virus type I enhancer binding protein 2 (HIVEP2) [42].

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.03×10⁻⁰⁵) 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.5×10⁻³): protein regulator of cytokinesis 1 (PRC1) [42], protein tyrosine phosphatase receptor type O (PTPRO) [41], ring finger protein 146 (RNF146) [27] and Ras-related GTP binding D (RRAGD) [41].

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 5 ky and possibly as far as 70 ky [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₂ sat [13], the HA Oromo sample resembles acclimatized lowlanders with a response characterized by elevated Hb concentration and marked reduction in O₂ 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 genome-wide 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₂ sat, the two phenotypes that have been most widely studied in highlander populations. No evidence of a genetic contribution to O₂ sat in Amhara, Oromo, and the combined Ethiopian sample could be detected (Figures S7, S8, S9, S10, S11, S12, S13, S14, S15). This is true also in the Tibetans, even though segregation analysis detected a major O₂ 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₂ sat is lower. This is consistent with studies in Tibetans and Andeans showing a markedly lower heritability for O₂ sat compared to Hb levels; indeed, the O₂ 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₂ 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 versus 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 versus 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₂ 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₂ 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 variables – see below) value across individuals and running the same linear model. The results were visualized by means of quantile-quantile (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²>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₁, X₂, …, 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 = β₀+β₁X_i1+ β₂X_i2+…+β_pX_ip+ε_i. Using the data from all genotyped SNPs that overlap among the p populations and applying multivariate linear regression, we found b₀, b₁, b₂, …, b_p, which are estimates of parameters β₀, β₁, β₂, …, β_p and estimated the magnitude of the residual for the i^th SNP as |ε_i| = |y_i−(b₀+b₁X₁+b₂X₂+…+b_pX_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 versus LA Amhara and (2) HA versus 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.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21

Attachment 22

Attachment 23

Attachment 24

Attachment 25

Attachment 26

Attachment 27

Attachment 28

Attachment 29

Attachment 30

Attachment 31

Attachment 32

Attachment 33

Attachment 34

Attachment 35

Attachment 36

Attachment 37

Attachment 38

Attachment 39

Attachment 40

Attachment 41

Attachment 42

Attachment 43

Attachment 44

Attachment 45

Attachment 46

Attachment 47

Attachment 48

Attachment 49

Attachment 50

Attachment 51

Attachment 52

Attachment 53