Genetic Variations in the TP53 Pathway in Native Americans Strongly Suggest Adaptation to the High Altitudes of the Andes

Abstract

The diversity of the five single nucleotide polymorphisms located in genes of the TP53 pathway (TP53, rs1042522; MDM2, rs2279744; MDM4, rs1563828; USP7, rs1529916; and LIF, rs929271) were studied in a total of 282 individuals belonging to Quechua, Aymara, Chivay, Cabanaconde, Yanke, Taquile, Amantani, Anapia, Uros, Guarani Ñandeva, and Guarani Kaiowá populations, characterized as Native American or as having a high level (> 90%) of Native American ancestry. In addition, published data pertaining to 100 persons from five other Native American populations (Surui, Karitiana, Maya, Pima, and Piapoco) were analyzed. The populations were classified as living in high altitude (≥ 2,500 m) or in lowlands (< 2,500 m). Our analyses revealed that alleles USP7-G, LIF-T, and MDM2-T showed significant evidence that they were selected for in relation to harsh environmental variables related to high altitudes. Our results show for the first time that alleles of classical TP53 network genes have been evolutionary co-opted for the successful human colonization of the Andes.

About the Authors:

Vanessa Cristina Jacovas

Contributed equally to this work with: Vanessa Cristina Jacovas, Diego Luiz Rovaris

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

Diego Luiz Rovaris

Contributed equally to this work with: Vanessa Cristina Jacovas, Diego Luiz Rovaris

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

Orlando Peréz

Affiliation: Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET), Puerto Madryn, Argentina

Soledad de Azevedo

Affiliation: Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET), Puerto Madryn, Argentina

Gabriel Souza Macedo

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

José Raul Sandoval

Affiliation: Facultad de Medicina Humana, Universidad de San Martin de Porres (USMP), Lima, Peru

Alberto Salazar-Granara

Affiliation: Facultad de Medicina Humana, Universidad de San Martin de Porres (USMP), Lima, Peru

Mercedes Villena

Affiliation: Instituto Boliviano de Biología de Altura (IBBA), Universidad Mayor de San Andres, La Paz, Bolivia

Jean-Michel Dugoujon

Affiliation: Anthropologie Moléculaire et Imagerie de Synthèse, CNRS UMR 5288, Université Paul Sabatier Toulouse III, Toulouse, 31000, France

Rafael Bisso-Machado

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

Maria Luiza Petzl-Erler

Affiliation: Laboratório de Genética Molecular, Departamento de Genética, Universidade Federal do Paraná, Curitiba, PR, Brasil

Francisco Mauro Salzano

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

Patricia Ashton-Prolla

Affiliations Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil, Serviço de Genética Medica, Hospital de Clínicas de Porto Alegre, Porto Alegre, RS, Brazil

Virginia Ramallo

Affiliation: Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET), Puerto Madryn, Argentina

Maria Cátira Bortolini

* E-mail: [email protected]

Affiliation: Departamento de Genética, Instituto de Biociências, Universidade Federal do Rio Grande do Sul, Porto Alegre, Rio Grande do Sul, Brasil

Introduction

The product of the TP53 gene is a transcription factor (p53) that activates or represses a large number of target genes that regulate a broad array of extremely important cellular functions, such as cell cycle, metabolism, DNA repair, senescence, and apoptosis. This factor is therefore essential for maintaining genome integrity [1]. In humans, p53 has 393 amino acids and the TP53 gene is located in the short arm of chromosome 17 [2]. Alterations of the TP53 gene or perturbations in the TP53 pathway are frequently correlated with carcinogenesis; more than 50% of human tumors carry mutations in this gene [3].

The steady-state levels of p53 are primarily determined by the rate at which it is degraded, rather than the rate at which it is produced. The TP53 gene is constitutively expressed in all cell types, but p53 does not accumulate in non-stressed cells, since it is rapidly degraded by the proteasome via ubiquitination [4, 5]. On the other hand, the p53 levels increase in response to various stress signals, such as UV irradiation, low oxygen concentrations (hypoxia), and exposure to high temperatures [6, 7, 8, 9].

There are many polymorphisms described for TP53, but a C→G non-synonymous substitution (rs1042522: c.215C>G, p. Pro72Arg; [10]) that promotes the amino acid change Pro→Arg at codon 72 of p53 is one of the most widely studied. This polymorphism has been described to be associated with an increased risk for developing cancer, since the p53-Pro allele is less active than p53-72Arg in inducing apoptosis, among other characteristics [11, 12].

Proper p53 transcriptional function is strongly linked to the activity of several other proteins encoded by the genes MDM2 (Mouse double minute 2 homolog; OMIM 164785), MDM4 (Mouse double minute 4 homolog, OMIM 602704), and USP7 (Ubiquitin-specific protease 7; OMIM 602519), also known as HAUSP (Herpesvirus-associated ubiquitin-specific protease). Another important gene in the so-called classical TP53 network [13] is LIF (Leukemia-inhibitory factor; OMIM 159540), which plays an essential role in the early phases of embryonic development in humans, and is regulated by p53 (S1 Fig).

The E3 ubiquitin-protein ligase, MDM2, mediates the activity of p53 by directing it to degradation by the proteasome [5, 14, 15]. MDM2 expression is also tightly regulated by p53 [16]. This auto-regulatory loop allows for the precise regulation of protein levels and activities of both p53 and MDM2 proteins [4, 17, 18].

The most well studied polymorphism in the MDM2 gene (rs2279744: c.14+309T>G) is located in its internal promoter. It consists of a single-nucleotide change from T→G, which increases the affinity of a sequence in MDM2 for the Sp1 transcription factor (Specificity protein 1; OMIM 189906). As a result, homozygotes for the G allele express more MDM2 than homozygotes for the T allele [19, 20]. In the presence of high levels of MDM2, there is a corresponding decrease of p53, causing a reduced response to cellular stress, impaired DNA repair, decreased apoptosis, and senescence [19]. Some studies have demonstrated that MDM2-309T and MDM2-309G alleles have different distributions in human populations [18, 21]. For instance, derived allele MDM2-309-G has higher frequency in European and Asian than African populations (average values: ∼0.35, ∼0.70, and ∼0.03, respectively; [22, 23, 24]). This allele may compensate for the higher apoptotic frequencies caused by the prevalence of allele p53-72Arg in Eurasians (∼0.56; [22, 23, 24]), suggesting adaptation [18].

The MDM4 protein, encoded by the MDM4 gene acts as a negative regulator of p53, inhibiting its transcriptional activity [25, 26, 27]. MDM2 and MDM4 form heterodimers with a high capacity for ubiquitination of target proteins, thus leading to degradation of targets, like p53 [28]. Deletion of either MDM2 or MDM4 induces p53-dependent early embryonic lethality in an animal model [16, 29]. The AA genotype for the single-nucleotide MDM4 polymorphism (rs1563828:g.204547449A>G) was associated with an increased risk for breast cancer [30].

Another important regulator of p53 is USP7, encoded by the USP7 gene, which deubiquitylates p53 and protects it from proteasome degradation [31]. The USP7 gene has a G→A substitution in intron 25 (rs1529916: g.8897333G>A), and derived allele A has been associated with endometriosis, female infertility, and prostate cancer [13, 32].

LIF is a cytokine expressed in various cell types, and its main function is to strengthen the blastocyst training of human embryos. In the very first days post-fertilization, LIF expression increases in the endometrium, creating a favorable environment for blastocyst implantation. Allele G of LIF (T→G transversion at the 3′ UTR region of the gene; rs929271: g.30242237T>G) is associated with female infertility [13]. LIF expression level is also known to be 2 times lower in cells bearing the p53-72Pro allele, compared to p53-72Arg, which can lead to the decrease of the implantation and fertility rate. In summary, several studies have strongly suggested that polymorphisms in the p53 signaling pathway play an important role in blastocyst implantation and are associated with recurrent pregnancy loss [13, 33, 34].

The genetic variability observed in contemporary human populations and the functionalities associated with the polymorphisms described above allow us to infer that a simple neutral model of mutation and drift is insufficient to explain the allelic distributions observed. Thus, it has been suggested that positive selection contributed to adaptation of Homo sapiens in different ecosystems. For example, the p53-72Arg allele (rs1042522) is more common in Europeans than in Africans, leading to the hypothesis that its distribution is dependent on latitude and maintained by selective pressures [35, 36]. On the other hand, Shi et al. [23] found that winter temperatures and UV radiation correlated significantly with the TP53 (rs1042522) and MDM2 (rs2279744) allele distributions in East Asian populations, indicating the possibility of adaptation to distinct environments.

America was the last continent occupied by humans in pre-colonial times. González-José et al. [37] and Bortolini et al. [38] suggested that an initial major dispersal began after 21,000 years before present, and that the biological and cultural characteristics of the first Americans that emerged, in part, were reshaped by recurrent trans-Beringian/circum-Arctic gene flow and important local population dynamics during a standstill period in Beringia. For example, Native Americans have experienced dramatic episodes of genetic drift and successive bottleneck events during migration across the continent. Furthermore, signals of positive natural selection associated to autochthonous environmental and cultural conditions have also been described [39, 40, 41].

Based on these findings, we hypothesize that the allele distributions of the classical TP53 pathway genes in Native American populations reflect adaptation, not only demographic and/or random events. To test our hypothesis, we determined the genotypes of the five above-mentioned SNPs in 282 unrelated individuals and compared the results to a large number of climate-related environmental variables, such as altitude, temperature, and seasonal mean UV radiation. Additional data regarding two of these SNPs (TP53-rs1042522 and MDM2-rs2279744) were compiled from the literature for a more extensive population analysis.

Materials and Methods

Samples and ethical procedures

Five SNPs (rs929271, rs1042522, rs1563828, rs2279744, and rs1529916) were genotyped in 282 volunteers characterized as Native American or as having large (> 90%; [42]) Native American ancestry. Volunteers were from 12 populations located in different ecoregions, namely highland (populations located at altitudes ≥ 2,500 m; [43]) and lowland (populations located at altitudes below 2,500 m). Highland populations were Aymara (n = 18) and Quechua (n = 17) from Bolivia, and Chivay (n = 18), Cabanaconde (n = 17), Yanke (n = 10), Taquile (n = 43), Amantani (n = 29), Anapia (n = 15), and Uros (n = 22) from Peru. All highland populations were located in the Andean region, including on Lake Titicaca islands or in their vicinity. Lowland populations were Andoas (n = 61), a Native Amazonian population living in North Peru, and Guaraní Indians from Brazil (Tupian speakers from two sub-groups: Ñandeva, n = 16; and Kaiowa, n = 16). Details about these populations have been summarized elsewhere [42, 44, 45, 46]. To facilitate the presentation of the results and discussion, we will collectively refer to all communities as “Native Americans”. The geographical coordinates (latitude and longitude) of all populations are presented in S1 File (Table A in S1 File).

Ethical approval for the use of these samples was obtained from the National Ethics Committee of Brazil (Resolution No. 123/98 CONEP) for individuals from Brazilian tribes; and by the Ethics Committee of Universidad San Martín de Porres, Lima, Peru (Peruvian samples) and Université Paul Sabatier Toulouse, Toulouse, France (Bolivian samples). Written informed consent or verbal informed consent (illiterate persons) was obtained individually from tribal participants. Verbal informed consent was registered in the field, and the institutional review ethics committees approved this procedure. This study was carried out in accordance with the Declaration of Helsinki.

Data from literature

Data from 100 additional individuals from five other Amerindian populations (Surui and Karitiana (Brazil), Piapoco (Colombia), Maya and Pima (México) were included in this study. For more details on these samples, please refer to http://www.cephb.fr/HGDP-CEPH-Panel/[47]. The environmental conditions evaluated for all populations (present study and literature sample) are compiled in S1 File (Table A in S1 File). The environmental data were collected for each population using the SoDa Service and WorldClim (http://www.soda-is.com/[48] and http://www.worldclim.org/[49], respectively; last access: December 19, 2014).

All analyses were performed with two sets of data: (A) 12 South American populations, for which original data regarding five SNPs (rs1042522, rs2279744, rs1529916, rs1563828, and rs929271) were obtained in the present study, and (B) all populations genotyped in this study plus five additional populations, for which TP53 rs1042522 and MDM2 rs2279744 data are available in the HGDP-CEPH panel [22].

Laboratory methods

Genomic DNA was obtained from saliva, whole blood, or plasma, using the QIAamp DNA extraction Mini kit (Qiagen; https://www.qiagen.com/br/[50]) according to manufacturer’s instructions. Genotyping of the TP53-rs1042522, MDM4-rs1563828, USP7-rs1529916, LIF-rs929271, and MDM2-rs2279744 SNPs was performed by allelic discrimination using the TaqMan Genotyping Assays (Applied Biosystems; http://www.lifetechnologies.com/br/en/home/brands/applied-biosystems.html [51]). Genotyping of MDM2 rs2279744 was performed using a customized (assay-by-design) assay using probes FAM-TCCCGCGCCGCAG and VIC-CTCCCGCGCCGAAG, with primers 5′-CGGGAGTTCAGGGTAAAGGT-3′ (forward) and 5′-ACAGGCACCTGCGATCATC-3′ (reverse).

PCR reactions were carried out in 48-well plates, with each reaction containing: 10 ng of genomic DNA, 2× TaqMan® genotyping Master Mix (Applied Biosystems), specific probes for each SNP (40×), and ultra-pure water for a final reaction volume of 10 μL. The PCR conditions were as follows: 95°C for 10 min, followed by 45 cycles of 95°C for 15 s and 63°C for 60 s. MDM2 rs2279744 genotyping was also done in 48-well plates, with each reaction containing: 10 ng of genomic DNA, 2× TaqMan® genotyping Master Mix, 5 μM of each primer and probe, and water to reach a final volume of 10 μL. MDM2 PCR conditions were as follows: 50°C for 2 min, 95°C for 10 min, and 45 cycles of 95°C for 15 s and 60°C for 60 s. All reactions were performed in an Illumina Eco Real-Time PCR System, (http://www.uniscience.com/ [52]) and results were analyzed using an Eco Real-Time PCR System and the Software v5.0 associated with that system. All wet-lab analyses were performed in the Laboratory of Human and Molecular Evolution of the Department of Genetics at Federal University of Rio Grande do Sul in Brazil.

Statistical analyses

Hardy-Weinberg equilibriums were calculated using a web-based program (http://www.oege.org/software/hwe-mr-calc.shtml [53]), and the statistical significance was assessed by Chi-square tests (p < 0.01). Analysis of molecular variance (AMOVA) using Arlequin 3.5.1.227 was applied to assess the variance among and within the investigated Native American populations [54, 55, 56].

Allele distributions were tested for possible associations with three groups of environmental conditions: 1) geographic: altitude, latitude, and longitude; 2) annual and seasonal mean UV radiation, and 3) Nineteen climate-related variables (Table A in S1 File). Principal component analysis (PCA) was performed to convert the nineteen possibly correlated bioclimatic variables into a smaller number of artificial variables (PCs) accounting for most of the variance in the previously observed variables. The correlation analysis between allele frequencies in each population and the environmental conditions was performed using Spearman´s rho correlation coefficient. The association between SNPs and altitude was assessed through binary logistic regression using two geographic categories (highlands: ≥ 2,500 m, and lowlands: < 2,500 m [43]) as the outcome and SNPs as predictors. Since this analysis was not intended to infer causality relationships, the odds ratio was reported as an estimate of size effect. For these analyses a Bonferroni correction was performed and the alpha was set at 0.01 (α_Bonferroni = 0.05/5 SNPs tested). Additionally, we performed the nonparametric Multifactor Dimensionality Reduction (MDR, v3.0.2; [57]) approach to detect potential gene–gene interactions. Thus, we used MDR to incorporate information from our 5 and 2 selected loci (data sets A and B, respectively) and an environmental condition as the outcome (altitude: highland and lowland geographic categories). The percentage of information gain (IG) by each SNP is visualized for each node, while the IG for each pairwise connection between SNPs is visualized for each branch. Thus, the independent main effects of each SNP can be compared to the interaction effect. The p-value was calculated based on 10,000 permutations.

Results

Table 1 shows the derived allele frequency for each SNP investigated (individual genotypes can be seen in S2 File). Wide variations were observed in some allele frequencies in both population groups (highland and lowland). For instance, the frequency of MDM2-309-G is about five times higher in Guaraní Ñandeva than Guaraní Kaiowa, which may reflect genetic drift since the split of these two Guaraní partialities occurred less than 2,000 years ago [45]. On the other hand, several highland populations from Peru and Bolivia present similar distributions of MDM2-309-G. Most of these highland populations show deviations from the Hardy-Weinberg equilibrium (HWE), especially in Peruvian samples for the MDM2-309 locus (Table B in S1 File).

[Figure omitted. See PDF.]

Table 1. Derived¹ allele frequencies and AMOVA results.

https://doi.org/10.1371/journal.pone.0137823.t001

AMOVA analysis, using both data sets (Table 1), indicated that homogeneity and population structures could be seen in both highland and lowland populations. For instance, population structure measured by F_ST statistics (i.e. the among-populations component of genetic variance) is observed in the two groups considering TP53 rs1042522 (F_ST = 0.068 and 0.054, for highland and lowland, respectively), while for MDM2 rs2279744 homogeneity is observed in highland populations (F_ST = −0.020; p = 0.801) while high heterogeneity is observed in lowland populations (F_ST = 0.274; p < 0.001). Only the F_ST value observed for LIF rs929271 in the highland group (11.8%) is similar to the average estimated across the human genome (12%; [58]). For F_CT (between-groups component of variance), the variance is high (11%) for MDM2 rs2279744 data, indicating a remarkable and significant difference between the allelic distributions of the highland and lowland populations.

Principal component analysis

In data set A, the first principal component (PC1) accounted for 73% of total variance, comprising the following bioclimatic variables: annual mean temperature, mean diurnal range, maximum temperature of warmest month, minimum temperature of coldest month, temperature annual range, mean temperature of wettest quarter, mean temperature of driest quarter, mean temperature of warmest quarter, mean temperature of coldest quarter, annual precipitation, precipitation of wettest month, precipitation in the driest month, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter, and precipitation of coldest quarter. The second principal component (PC2) represented 13% of variance, and comprised temperature seasonality, which is a measure of standard deviation × 100 of average annual daily temperatures.

When we expanded our analysis to data set B, PC1 represented 59% of total variance and comprised sixteen bioclimatic variables: annual mean temperature, mean diurnal range, maximum temperature of warmest month, minimum temperature of coldest month, mean temperature of wettest quarter, mean temperature of driest quarter, mean temperature of warmest quarter, mean temperature of coldest quarter, annual precipitation, precipitation of wettest month, precipitation in the driest month, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter, and precipitation of coldest quarter. The second principal component (PC2) represented 23% of variance, and comprised isothermality (the ratio of mean diurnal range to temperature annual range), temperature seasonality, and temperature annual range, all of which are connected with climatic changes by seasonality.

Correlation analyses

Correlation coefficients and their statistical significances are given in Table 2. In data set A, there were significant associations between the G allele of USP7 (rs1529916) and the annual mean of ultraviolet irradiance (rho = 0.760 p = 0.004) and PC1 (rho = −0.741, p = 0.006). This allele was also nominally associated with the mean of ultraviolet irradiance in the coldest semester (rho = 0.681 p = 0.015) and in the warmest semester (rho = 0.618 p = 0.032). MDM2 (rs2279744) T allele was nominally associated to longitude (rho = −0.587 p = 0.045) and the mean of ultraviolet irradiance in the coldest semester (rho = 0.605 p = 0.037), while LIF (rs929271) T allele was nominally associated to annual mean of ultraviolet irradiance (rho = 0.693 p = 0.013) and PC1 (rho = −0.664, p = 0.018).

[Figure omitted. See PDF.]

Table 2. Correlation analysis results.

https://doi.org/10.1371/journal.pone.0137823.t002

In data set B, there were significant associations between the T allele of MDM2 (rs2279744) and altitude (rho = 0.673, p = 0.003), the mean of ultraviolet irradiance in the coldest semester (rho = 0.827 p < 0.001), and PC1 (rho = −0.610 p = 0.009). This allele was also nominally associated to PC2 (rho = −0.567, p = 0.018).

Binary logistic regression analyses

We performed a binary logistic regression analysis to search for possible associations between SNPs and two geographic categories (Highlands: ≥ 2,500 m; Lowlands: < 2,500 m) using altitude as dependent variable (Table C in S1 File). In data set A, we observed statistically significant associations for USP7 rs1529916 and LIF rs929271 SNPs. Individuals who inhabit the highlands were less likely to carry USP7-GA (OR = 0.417, p = 0.002) and USP7-AA (OR = 0.135, p < 0.001) genotypes. A similar association was observed for LIF-TG (OR = 0.324, p = 0.001 and LIF-GG (OR = 0.270, p < 0.001) genotypes. Regarding data set B, an association between the MDM2 rs2279744 SNP and altitude was detected. Individuals who inhabit the highlands were less likely to carry MDM2-TG (OR = 0.218, p < 0.001) and MDM2-GG (OR = 0.175, p < 0.001) genotypes.

Gene-gene interaction analyses

We used the MDR approach to search for gene-gene interactions (Table D in S1 File). These analyses were intended to explore differences between highland and lowland populations in genotype combinations among the SNPs investigated since gene networks, such as those investigated here, can be sources of epistasis. Significant two- (p = 0.004) and three-locus interactions (p = 0.004) were identified in data set A. However, an analysis of IG based on entropy measures revealed that these effects were not explicated by epistasis (negative values in the branches among nodes; Fig 1A). On the other hand, IG values of both USP7 (7.26%) and LIF (4.23%) indicated that both genes have a large main effect in a scenario where altitude is considered, corroborating our previous analysis. Regarding data set B, the largest main effect was observed for MDM2 (IG = 9.51%), which contrasts with the low value for TP53 (IG = 0.41%). A potential synergism (epistasis) between the two loci was also found, but it is apparently weak (IG value of only 1.54%; Fig 1A), at least when it is compared with the potential mechanism of action on MDM2. On the other hand, it is 3.7 times greater than the main effect of TP53. It is noteworthy that independent of the TP53 genotype, the genotype MDM2-TT is always favorable and most commonly found in highlands (Fig 1B). In other words, MDM2 showed the greatest contribution to adaptation to hostile environments, such as those found in the highlands.

[Figure omitted. See PDF.]

Fig 1. Summary of the multifactor dimensionality reduction (MDR) interaction models.

(A) Interaction graphs comprised of nodes with pairwise connections between them. Values in nodes represent information gain (IG) of individual genes (main effect), while values between nodes are the IG of each pairwise combination (interaction effects). Positive entropy (plotted in red) indicates interaction (epistasis) and negative entropy (plotted in green or blue) indicates redundancy. Independence is represented by the gold color. (B) The MDM2-TP53 interaction associated with altitude in data set B. High-frequency genotype combinations in individuals who inhabit highlands (≥ 2,500 meters) are depicted as darkly shaded cells and low-frequency combinations in those individuals as lightly shaded. For each cell, the left bar indicates the absolute number of individuals who inhabit highlands and the right bar the absolute number of individuals who inhabit lowlands (< 2,500 meters).

https://doi.org/10.1371/journal.pone.0137823.g001

Discussion

More than 60,000 scientific studies have been published in the last 30 years concerning the roles of TP53 network genes, as well as of their variants, in human susceptibility to cancer and other pathological conditions. Special issues in scientific journals, dedicated to these topics, have also been published (see, as example [59]). This overwhelming number of studies contrasts with the rarity of studies of an evolutionary context, which are indispensable for explaining differences in the TP53 network allele distributions along human populations, which often cannot be understood as simply a result of stochastic processes. Our goal here was to help fill this gap, providing information about five polymorphisms of the classical TP53 network in Native American populations and how their variability patterns could be explained.

Our analysis of data set A, which included original information of 5 SNPs in 12 Native American populations, suggests a well-known role of genetic drift in those groups, illustrated by wide difference in MDM2-G allele frequencies between the two Guaraní sub-groups. However, other instigating results can be associated to adaptation to environmental conditions in Native American populations. Alleles USP7-G (rs1529916) and LIF-T (rs929271) were correlated with ultraviolet irradiance and index of temperature and precipitation, variables comprising PC1. Additionally, examining variables with the highest representation in the PC1 components (> 0.90), it is possible to see that in regions where the annual mean temperatures, minimum temperatures of the coldest month, mean temperatures of the driest quarter, mean temperatures of the coldest quarter, and precipitation are low, the presence of ancestral alleles G and T are significantly higher. In other words, our analysis as whole reveals that alleles USP7-G and LIF-T are more highly represented in stressful environments (low temperature, arid climate, wide temperature range during the day, and high levels of UV radiation), which is typical of high altitudes. It is noteworthy that derived alleles of these SNPs have been associated with cancer susceptibility, infertility, and endometriosis [13, 32], so that the alleles USP7-G and LIF-T could be considered as protective factors against the consequences of harsh environmental stress.

Human populations living at high altitudes are likely to have developed specific adaptations to support both the harsh conditions described above and low oxygen concentrations (hypoxia; [41]). Monge in 1948 [60] proposed that the hypoxia could reduce fertility in humans. However, recent studies have shown that the reproductive functioning of populations indigenous to high altitudes is adapted to hypoxia and other extreme conditions [61]. Our results with USP7 (rs1529916) and LIF (rs929271) polymorphisms could be connected with adaptation of the reproductively successful ancestors of modern Andes populations.

In examining data set B, we found the ancestral MDM2-T allele is strongly correlated with winter mean UV radiation, altitude, and PC1. The highest representations in the PC1 components (> 0.90) are annual mean temperature, minimum temperature of coldest month, minimum temperature of coldest quarter, and annual precipitation. Allele T is significantly more frequent in communities located at high altitudes experiencing extreme environmental conditions, such as high UV radiation and dry and cold climate. In addition, the binary logistic regression analysis showed that MDM2-TT individuals are more frequently found in highlands. MDM2-TT homozygotes express typical steady-state levels of MDM2, maintaining an adequate level of p53 [20], and consequently can appropriately respond environmental stresses. An important confounding factor could be admixture with Europeans, which is more important in Andean than in the lowland populations considered here [42, 62, 63]. However, any effect of admixture would be in the opposite direction, since MDM2-G frequency is relatively high in Spaniards (0.39; [24]).

The inverse correlation between MDM2-T frequencies and winter UV radiation is consistent with the findings of Shi et al. [23], which showed that low levels of UV are significantly correlated with genotype MDM2-GG in Han Chinese populations, similarly deviating from HWE. These authors suggested that MDM2-GG is selected for in areas of low UV activity (at high altitudes, the thinner atmosphere will filter less UV radiation; consequently for every 1000 m increase in altitude, the UV radiation level will increase ∼12%; http://www.weather.gov.hk/radiation/tidbit/201012/uv_e.htm [64]). Natural selection can be evoked to explain these results, although the HWE test is considered too weak to detect this phenomenon.

As mentioned above, native Andean populations have successfully adapted to environments with low oxygen concentrations. One gene that contributes to hypoxia adaptation is EPAS1 (Endothelial PAS domain-containing protein 1, also known as HIF-2α, Hypoxia-inducible factor—alpha 2 (OMIM 603349)), which acts by preventing toxicity promoted by hypoxia. This gene plays an important role in both the classical and the expanded TP53 network. For instance, the alpha subunit of EPAS1 regulates p53 activity, including through prevention of damage-induced degradation and nuclear export of MDM2, stabilizing nuclear p53 [65]. Foll et al. [41] confirmed the action of positive selection on EPAS1 in both Tibetans and Andeans. Furthermore, several studies have revealed a role for p53 and its regulation in physiological and metabolic processes resulting from environments with low oxygen concentrations [8, 66, 67]. Recently, Eichstaedt et al. [68] studied an indigenous population living in the Argentinean Andes (Colla) and identified signatures of positive selection in genes involved in cellular hypoxia, including TP53. Importantly, hypoxia induces p53 accumulation through down-regulation of MDM2 [66]. These results reinforce our suggestion that individuals with the MDM2-TT genotype represent an adaptation to the environmental stresses of high altitudes. In addition, the interaction analysis performed by the MDR method using both data sets (A and B) revealed the potential for the MDM2, LIF, and USP7 genes to play an additional central role in a high altitude setting. Thus, taken together, our results demonstrate that variation of the p53-activating stressors could not be directly correlated with p53-Pro72Arg alleles, but with frequencies of the other functional polymorphisms examined, such as USP7-G (rs1529916), LIF-T (rs929271), and MDM2-309, as well as synergic interactions between them.

Under neutral model conditions, South Amerindians living in lowlands present higher levels of population structure when compared to those seen in indigenous Andean communities [62, 69]. However, not all F_ST values obtained in our study were consistent with this expectation (Table 1). Positive selection disturbs the patterns of genetic variation expected under a standard neutral model [70]. Additionally, it is possible to see that some derived alleles, such as MDM2-G, have high frequencies in Asian populations with putative common ancestry (0.57–0.82; [22, 23, 24]), but a surprisingly low distribution in Andeans (average value: ∼0.13). An excess of unexpectedly low and/or high frequencies of derived alleles can also be considered a marker of positive selection [70]. Thus, the distributions of the classical TP53 pathway alleles in Native American populations could be under selective pressure. Sucheston et al. [22] investigated 52 worldwide populations from the HGDP-CEPH-panel for the prevalence of p53-Pro72Arg and MDM2-309 polymorphisms, but found no significant association with climate variables. However, the Native American samples in the Sucheston et al.’ study [22] were much smaller than the present study (see Table 1), which may explain the divergent results.

Finally, government surveys in Peru indicate that the rate of gestational and postpartum complications in Aymara regions is lower than the national average (1.8% and 5% respectively; http://www.dge.gob.pe/publicaciones/pub_asis/asis26.pdf, p. 165; [71] http://www.dge.gob.pe/portal/docs/intsan/asis2012.pdf, p. 76 [72]). These same official sources also indicate differences in the cancer incidences between lowland localities and some regions situated at high altitude (for example in the Puno state, where the Anapia community is located; http://www.dge.gob.pe/portal/docs/asis_cancer.pdf, p. 64 [73]). These findings are in agreement with our genetic results. However, only additional and specific studies can accurately relate our evolutionary findings with those related to the health of contemporary Andean populations.

A well-regulated p53 network is crucial for maintaining genomic integrity. Several polymorphisms in this pathway have been described, and the different allele frequencies among human populations have been interpreted as the result of selective pressure. Humans occupied high-altitude locations in the Andes as early as 12,800 years ago, providing a sufficient period of time for the initiation of organismal selection and developmental functional adaptation ([74] and references therein). Here we are suggesting that natives from Andes, who are subjected to low temperatures, arid climates, wide temperature ranges during the day, high levels of UV radiation, and hypoxia, among other environmental insults, are protected by a selected combination of alleles/genotypes of the TP53 pathway. The present study identifies for the first time the potential role of the MDM2, LIF, and USP7 in the adaptation of the Andean populations.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. The p53 network.

Network view of p53 pathway analyzed by STRING 10.0 (http://string-db.org/). Interaction confidence score cutoff was 900 (highest confidence). Each color arrow represents a predicted functional partner: green (activation), red (inhibition), blue (binding), purple (catalysis), pink (post-translational modification), black (reaction), and yellow (expression). TP53 = tumor protein p53, USP7 = ubiquitin specific peptidase 7 (herpes virus-associated), MDM4 = Mouse double minute 4 homolog, MDM2 = Mouse double minute 2 homolog, and LIF = leukemia inhibitory factor.

https://doi.org/10.1371/journal.pone.0137823.s001

(TIFF)

S1 File. Additional Results.

Climatic variables evaluated in population of this study (Table A). Allelic frequencies and Hardy-Weinberg Equilibrium results (Table B). Binary logistic regression analyses results (Table C). Locus interaction by the multifactor dimensionality reduction (MDR) approach (Table D).

https://doi.org/10.1371/journal.pone.0137823.s002

(DOC)

S2 File. Individual Genotype Database.

In this database, we found individual genotypes of all polymorphisms discussed in the study.

https://doi.org/10.1371/journal.pone.0137823.s003

(XLSX)

Acknowledgments

We are very grateful to the individuals who donated the samples analyzed here and to the Fundação Nacional do Indio (Brazil) for logistic support. We thank René Vasquez for his assistance in the sample collection in Bolivia and Sidia Maria Callegari Jacques for statistical support. We are grateful to professor David Comas for his careful review of this manuscript.

Author Contributions

Conceived and designed the experiments: VCJ VR MCB. Performed the experiments: VCJ VR. Analyzed the data: VCJ DLR OP VR. Contributed reagents/materials/analysis tools: SA GSM JRS AS-G MV J-MD RB-M MLP-E FMS PA-P. Wrote the paper: VCJ DLR VR FMS MCB.

Citation: Jacovas VC, Rovaris DL, Peréz O, de Azevedo S, Macedo GS, Sandoval JR, et al. (2015) Genetic Variations in the TP53 Pathway in Native Americans Strongly Suggest Adaptation to the High Altitudes of the Andes. PLoS ONE 10(9): e0137823. https://doi.org/10.1371/journal.pone.0137823

References

1. Botcheva K. p53 binding to human genome: crowd control navigation in chromatin context. Front Genet. 2014 Dec 22;5:447. pmid:25566329

2. Linzer DI, Levine AJ. Characterization of a 54K dalton cellular SV40 tumor antigen present in SV40-transformed cells and uninfected embryonal carcinoma cells. Cell. 1979 May;17(1):43–52. pmid:222475

3. Leroy B, Anderson M, Soussi T. TP53 mutations in human cancer: database reassessment and prospects for the next decade. Hum Mutat. 2014 Jun;35(6):672–88. pmid:24665023

4. Nag S, Qin J, Srivenugopal KS, Wang M, Zhang R. The MDM2-p53 pathway revisited. J Biomed Res. 2013 Jul;27(4):254–71. pmid:23885265

5. Chao CC. Mechanisms of p53 degradation. Clin Chim Acta. 2015 Jan 1;438:139–47. pmid:25172038

6. Ljungman M. Dial 9-1-1 for p53: mechanisms of p53 activation by cellular stress. Neoplasia. 2000 May-Jun;2(3):208–25. pmid:10935507

7. Latonen L, Taya Y, Laiho M. UV-radiation induces dose-dependent regulation of p53 response and modulates p53-HDM2 interaction in human fibroblasts. Oncogene. 2001 Oct 11;20(46):6784–93. pmid:11709713

8. Sermeus A, Michiels C. Reciprocal influence of the p53 and the hypoxic pathways. Cell Death Dis. 2011 May 26;2:e164. pmid:21614094

9. Chen J. The Roles of MDM2 and MDMX Phosphorylation in Stress Signaling to p53. Genes Cancer. 2012 Mar;3(3–4):274–82. pmid:23150760

10. Matlashewski GJ, Tuck S, Pim D, Lamb P, Schneider J, Crawford LV. Primary structure polymorphism at amino acid residue 72 of human p53. Mol Cell Biol. 1987 Feb;7(2):961–3. pmid:3547088

11. Dumont P, Leu JI, Della Pietra AC 3rd, George DL, Murphy M. The codon 72 polymorphic variants of p53 have markedly different apoptotic potential. Nat Genet. 2003 Mar;33(3):357–65. pmid:12567188

12. Whibley C, Pharoah PD, Hollstein M. p53 polymorphisms: cancer implications. Nat Rev Cancer. 2009 Feb;9(2):95–107. pmid:19165225

13. Kang HJ, Feng Z, Sun Y, Atwal G, Murphy ME, Rebbeck TR, et al. Single-nucleotide polymorphisms in the p53 pathway regulate fertility in humans. Proc Natl Acad Sci U S A. 2009 Jun 16;106(24):9761–6. pmid:19470478

14. Momand J, Zambetti GP, Olson DC, George D, Levine AJ. The mdm-2 oncogene product forms a complex with the p53 protein and inhibits p53-mediated transactivation. Cell. 1992 Jun 26;69(7):1237–45. pmid:1535557

15. Marine JC, Francoz S, Maetens M, Wahl G, Toledo F, Lozano G. Keeping p53 in check: essential and synergistic functions of Mdm2 and Mdm4. Cell Death Differ.2006 Jun;13(6):927–34. pmid:16543935

16. Eischen CM, Lozano G. The Mdm network and its regulation of p53 activities: a rheostat of cancer risk. Hum Mutat. 2014 Jun;35(6):728–37. pmid:24488925

17. Kohn KW, Pommier Y. Molecular interaction map of the p53 and Mdm2 logic elements, which control the Off-On switch of p53 in response to DNA damage. Biochem Biophys Res Commun. 2005 Jun 10;331(3):816–27. pmid:15865937

18. Atwal GS, Bond GL, Metsuyanim S, Papa M, Friedman E, Distelman-Menachem T, et al. Haplotype structure and selection of the MDM2 oncogene in humans. Proc Natl Acad Sci U S A. 2007 Mar 13;104(11):4524–9. pmid:17360557

19. Bond GL, Hu W, Bond EE, Robins H, Lutzker SG, Arva NC, et al. A single nucleotide polymorphism in the MDM2 promoter attenuates the p53 tumor suppressor pathway and accelerates tumor formation in humans. Cell. 2004 Nov 24;119(5):591–602. pmid:15550242

20. Alazzouzi H, Suriano G, Guerra A, Plaja A, Espín E, Armengol M, et al. Tumour selection advantage of non-dominant negative P53 mutations in homozygotic MDM2-SNP309 colorectal cancer cells. J Med Genet. 2007 Jan;44(1):75–80. Epub 2006 Jul 6. pmid:16825434

21. Millikan RC, Heard K, Winkel S, Hill EJ, Heard K, Massa B, et al. No association between the MDM2–309 T/G promoter polymorphism and breast cancer in African-Americans or Whites. Cancer Epidemiol Biomarkers Prev. 2006 Jan;15(1):175–7. pmid:16434608

22. Sucheston L, Witonsky DB, Hastings D, Yildiz O, Clark VJ, Di Rienzo A, et al. Natural selection and functional genetic variation in the p53 pathway. Hum Mol Genet. 2011 Apr 15;20(8):1502–8. pmid:21266458

23. Shi H, Tan SJ, Zhong H, Hu W, Levine A, Xiao CJ, et al. Winter temperature and UV are tightly linked to genetic changes in the p53 tumor suppressor pathway in Eastern Asia. Am J Hum Genet. 2009 Apr;84(4):534–41. pmid:19344876

24. 1000 Genomes Project Consortium, Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, et al. An integrated map of genetic variation from 1,092 human genomes. Nature. 2012 Nov 1;491(7422):56–65. pmid:23128226

25. Shvarts A, Steegenga WT, Riteco N, van Laar T, Dekker P, Bazuine M, et al. MDMX: a novel p53-binding protein with some functional properties of MDM2. EMBO J. 1996 Oct 1;15(19):5349–57. pmid:8895579

26. Shvarts A, Bazuine M, Dekker P, Ramos YF, Steegenga WT, Merckx G, et al. Isolation and identification of the human homolog of a new p53-binding protein, Mdmx. Genomics. 1997 Jul 1;43(1):34–42. pmid:9226370

27. Liu J, Tang X, Li M, Lu C, Shi J, Zhou L, et al. Functional MDM4 rs4245739 genetic variant, alone and in combination with P53 Arg72Pro polymorphism, contributes to breast cancer susceptibility. Breast Cancer Res Treat. 2013 Jul;140(1):151–7. pmid:23793604

28. Pei D, Zhang Y, Zheng J. Regulation of p53: a collaboration between Mdm2 and Mdmx. Oncotarget. 2012 Mar;3(3):228–35. pmid:22410433

29. Xiong S. Mouse models of Mdm2 and Mdm4 and their clinical implications. Chin J Cancer. 2013 Jul;32(7):371–5. pmid:23327795

30. Song CG, Fu FM, Wu XY, Wang C, Shao ZM. Correlation of polymorphism rs1563828 in MDM4 gene with breast cancer risk and onset age. Zhonghua Wai Ke Za Zhi. 2012 Jan 1;50(1):53–6. pmid:22490292

31. Shan J, Brooks C, Kon N, Li M, Gu W. Dissecting roles of ubiquitination in the p53 pathway. Ernst Schering Found Symp Proc. 2008;(1):127–36. pmid:19202598

32. Sun T, Lee GS, Oh WK, Pomerantz M, Yang M, Xie W, et al. Single-nucleotide polymorphisms in p53 pathway and aggressiveness of prostate cancer in a Caucasian population. Clin Cancer Res. 2010 Nov 1;16(21):5244–51. pmid:20855462

33. Dd Paskulin, Paixão-Côrtes VR, Hainaut P, Bortolini MC, Ashton-Prolla P. The TP53 fertility network. Genet Mol Biol. 2012 Dec;35(4 (suppl)):939–46.

34. Fraga LR, Dutra CG, Boquett JA, Vianna FS, Gonçalves RO, Paskulin DD, et al. p53 signaling pathway polymorphisms associated to recurrent pregnancy loss. Mol Biol Rep. 2014 Mar;41(3):1871–7. pmid:24435975

35. Beckman G, Birgander R, Själander A, Saha N, Holmberg PA, Kivelä A, et al. Is p53 polymorphism maintained by natural selection? Hum Hered. 1994; 44:266–270. pmid:7927355

36. Själander A, Birgander R, Saha N, Beckman L, Beckman G. p53 polymorphisms and haplotypes show distinct differences between major ethnic groups. Hum Hered. 1996 Jan-Feb;46(1):41–8. pmid:8825462

37. González-José R, Bortolini MC, Santos FR, Bonatto SL. The peopling of America: craniofacial shape variation on a continental scale and its interpretation from an interdisciplinary view. Am J Phys Anthropol. 2008 Oct;137(2):175–87. pmid:18481303

38. Bortolini MC, González-José R, Bonatto SL, Santos FR. Reconciling pre-Columbian settlement hypotheses requires integrative, multidisciplinary, and model-bound approaches. Proc Natl Acad Sci U S A. 2014 Jan 14;111(2):E213–4. pmid:24398530

39. Hünemeier T, Amorim CE, Azevedo S, Contini V, Acuña-Alonzo V, Rothhammer F, et al. Evolutionary responses to a constructed niche: ancient Mesoamericans as a model of gene-culture coevolution. PLoS One. 2012;7(6):e38862. pmid:22768049

40. Hünemeier T, Gómez-Valdés J, Ballesteros-Romero M, de Azevedo S, Martínez-Abadías N, Esparza M, et al. Cultural diversification promotes rapid phenotypic evolution in Xavánte Indians. Proc Natl Acad Sci U S A. 2012 Jan 3;109(1):73–7. pmid:22184238

41. Foll M, Gaggiotti OE, Daub JT, Vatsiou A, Excoffier L. Widespread signals of convergent adaptation to high altitude in Asia and america. Am J Hum Genet. 2014 Oct 2;95(4):394–407. pmid:25262650

42. Sandoval JR, Salazar-Granara A, Acosta O, Castillo-Herrera W, Fujita R, Pena SD, et al. Tracing the genomic ancestry of Peruvians reveals a major legacy of pre-Columbian ancestors. J Hum Genet. 2013 Sep;58(9):627–34. pmid:23863748

43. Moore LG. Human genetic adaptation to high altitude. High Alt Med Biol. 2001 Summer;2(2):257–79. pmid:11443005

44. Tsuneto LT, Probst CM, Hutz MH, Salzano FM, Rodriguez-Delfin LA, Zago MA, et al. HLA class II diversity in seven Amerindian populations. Clues about the origins of the Aché. Tissue Antigens. 2003 Dec;62(6):512–26. pmid:14617035

45. Marrero AR, Silva-Junior WA, Bravi CM, Hutz MH, Petzl-Erler ML, Ruiz-Linares A, et al. Demographic and evolutionary trajectories of the Guarani and Kaingang natives of Brazil. Am J Phys Anthropol. 2007 Feb;132(2):301–10. pmid:17133437

46. Gayà-Vidal M, Dugoujon JM, Esteban E, Athanasiadis G, Rodríguez A, Villena M, et al. Autosomal and X chromosome Alu insertions in Bolivian Aymaras and Quechuas: two languages and one genetic pool. Am J Hum Biol. 2010 Mar-Apr;22(2):154–62. pmid:19593738

47. Available: http://www.cephb.fr/HGDP-CEPH-Panel/. Accessed October 2014.

48. http://www.soda-is.com/. Accessed 19 December 19 2014.

49. http://www.worldclim.org/. Accessed 19 December 2014.

50. https://www.qiagen.com/br/

51. www.lifetechnologies.com/br/en/home/brands/applied-biosystems.html

52. http://www.uniscience.com/

53. http://www.oege.org/software/hwe-mr-calc.shtml,, last accessed January 2015.

54. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution 1984; 38: 1358–1370.

55. Weir BS. The second National Research Council report on forensic DNA evidence.bAm J Hum Genet. 1996 Sep;59(3):497–500.

56. Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992 Jun;131(2):479–91. pmid:1644282

57. Moore JH, Gilbert JC, Tsai CT, Chiang FT, Holden T, Barney N, et al. A flexible computational framework for detecting, characterizing, and interpreting statistical patterns of epistasis in genetic studies of human disease susceptibility. J Theor Biol. 2006 Jul 21;241(2):252–61. pmid:16457852

58. International HapMap Consortium. A haplotype map of the human genome. Nature. 2005 Oct 27;437(7063):1299–320. pmid:16255080

59. Soussi T. The TP53 gene network in a postgenomic era. Hum Mutat. 2014 Jun;35(6):641–2. pmid:24753184

60. Monge C. Acclimatization in the Andes. Baltimore: Johns Hopkins. University Press; 1948.

61. Vitzthum VJ. Fifty fertile years: anthropologists' studies of reproduction in high altitude natives. Am J Hum Biol. 2013 Mar-Apr;25(2):179–89. pmid:23382088

62. Wang S, Lewis CM, Jakobsson M, Ramachandran S, Ray N, Bedoya G, et al. Genetic variation and population structure in native Americans. PLoS Genet. 2007 Nov;3(11):e185. pmid:18039031

63. Verdu P, Pemberton TJ, Laurent R, Kemp BM, Gonzalez-Oliver A, Gorodezky C, et al. Patterns of admixture and population structure in native populations of Northwest North America. PLoS Genet. 2014 Aug 14;10(8):e1004530. pmid:25122539

64. http://www.weather.gov.hk/radiation/tidbit/201012/uv_e.htm. Accessed January 2015.

65. Chen D, Li M, Luo J, Gu W. Direct interactions between HIF-1 alpha and Mdm2 modulate p53 function. J Biol Chem. 2003 Apr 18;278(16):13595–8. pmid:12606552

66. Alarcón R, Koumenis C, Geyer RK, Maki CG, Giaccia AJ. Hypoxia induces p53 accumulation through MDM2 down-regulation and inhibition of E6-mediated degradation. Cancer Res. 1999 Dec 15;59(24):6046–51. pmid:10626788

67. Vousden KH, Ryan KM. p53 and metabolism. Nat Rev Cancer. 2009; 9:691–700. pmid:19759539

68. Eichstaedt CA, Antão T, Pagani L, Cardona A, Kivisild T, Mormina M. The Andean adaptive toolkit to counteract high altitude maladaptation: genome-wide and phenotypic analysis of the Collas. PLoS One. 2014 Mar 31;9(3):e93314. pmid:24686296

69. Tarazona-Santos E, Carvalho-Silva DR, Pettener D, Luiselli D, De Stefano GF, Labarga CM, et al. Genetic differentiation in South Amerindians is related to environmental and cultural diversity: evidence from the Y chromosome. Am J Hum Genet. 2001 Jun;68(6):1485–96. pmid:11353402

70. Biswas S, Akey JM. Genomic insights into positive selection. Trends Genet. 2006; 22:437–446. pmid:16808986

71. http://www.dge.gob.pe/publicaciones/pub_asis/asis26.pdf, p. 165;

72. http://www.dge.gob.pe/portal/docs/intsan/asis2012.pdf, p. 76;

73. http://www.dge.gob.pe/portal/docs/asis_cancer.pdf, p. 64;

74. Rademaker K, Hodgins G, Moore K, Zarrillo S, Miller C, Bromley GR, et al. Paleoindian settlement of the high-altitude Peruvian Andes. Science. 2014 Oct 24;346(6208):466–9. pmid:25342802