2.1 Study populations and community engagement

The Andean study groups have been described in detail [10,12]. The individuals from the Argentinean Andes were recruited for a study on genetics for arsenic metabolism and sampled from San Antonio de los Cobres (SAC; 3,800 m above sea level in the Puna region) in 2008. Before the recruitment, an information meeting with the health personnel at the local hospital was held. During interviews, the participants were asked, as an inclusion criterion, whether they were originally from SAC and whether their families had lived in the same region for at least 2–3 generations. Directly asking about ethnicity was avoided since it was a potentially sensitive topic in this region at the time of recruitment. Based on later research in the same area, however, people from SAC identified themselves as of Atacameño-Kolla origin [13]. The Atacameños have lived in this region for 11,000 years [14], and traces of human settlements in this area date back to 1,500 years BCE [15]. Since the participants did not identify themselves as belonging to a specific ethnicity at the time of the recruitment, we refer to this population based on their geographic location.

In the Bolivian Andes, we recruited individuals for a study on genetics for arsenic metabolism from ten villages around Lake Poopó, located in the southern part of the Bolivian Altiplano (3,686 m above sea level). Recruitment took place between 2015 and 2018 and involved a series of pre-information meetings with the health personnel and the local authorities, followed by information meetings with the health personnel at the time of recruitment. The study participants were identified as belonging to the Aymara-Quechua or Uru ethnicity based on the location of residency. The Aymara and Quechua are the predominant indigenous groups in this region. They are fairly similar from a genetic standpoint [16,17], and especially around Lake Poopó, the two ethnicities cohabit. Therefore, both ethnicities were assessed as a common study group. In contrast, the Uru communities around Lake Poopó, specifically of Uru-Murato descent, have been historically isolated from other communities in that region, and the number of Uru individuals living in this region is low [18]. We recruited participants from all Uru communities living around Lake Poopó, but not from all Aymara-Quechua communities.

The study areas in Bolivia and Argentina were revisited in 2020 to obtain oral and written informed consent from the participants for this study on whole-genome selection. In total, 81 individuals from the Bolivian Altiplano (65 Aymara-Quechua women, 13 Uru women, and 3 Uru men) and 32 women from the Argentinean Puna agreed to participate in the current study. Most study participants were women, as the men were often away from home to work. All samples and genetic data had been previously obtained and described [9,10]. This study was approved by the Comisión Provincial de Ética del Ministerio de Salud de Salta (Argentina), the Comité Nacional de Bioética (Bolivia), and the Regional Ethic Committee of Karolinska Institutet (Sweden, No. 2020-00495 and 2020-00493).

To return results to the indigenous communities who have contributed to the present study, a series of workshops will be held in several communities. In addition to describing the purpose and methodologies used, the workshops will also be devoted to contextualizing how this genetic research can be understood within their Andean culture.

2.2 Genome-wide SNP genotyping

Study participants were genotyped as described previously [9,10]. Briefly, DNA was extracted from whole blood samples or buccal cells with EZNA Blood DNA Mini kit (Omega Bio-teck, USA) or Qiagen Blood Mini kit (Qiagen, Germany). Genome-wide genotyping was performed at the SNP&Seq Technology Platform in Uppsala (Sweden) on the Illumina Infinium Omni5Exome and on the Illumina Infinium Omni5M bead chips for the Bolivian and Argentinean study groups, respectively. The data was aligned to the human reference genome build, version 37 (hg19).

Using PLINK v.1.90 [19], we filtered for a 5% genotype missingness threshold for SNPs, a 15% genotype missingness threshold for individuals, indels, and duplicates. Further, A/T and C/G SNPs were removed to avoid strand issues when merging to comparative datasets. This quality filtering was performed independently on the Bolivian and Argentinean datasets. Only autosomal SNPs were included.

2.3 Population structure analyses

Principal component analysis and ADMIXTURE were previously described [9] and included in this work in Figure S1 as a reference. Briefly, this dataset included the three Andean populations (Aymara-Quechua, Uru, and SAC) merged with other South American populations reported by Barbieri et al. (2019) [20] (n = 175) and complemented with selected populations of Native American ancestry from North, Central, and South America (n = 312) from a previously published dataset assembled by Lazaridis et al. (2014) [21]. A Hardy-Weinberg equilibrium filter (p-value < 0.001) was used to account for genotyping errors in the Bolivian and Argentinean datasets. First-degree relatives, as determined with the software KING v.2.1.4 [22], were excluded. In total, this dataset had 236,127 autosomal SNPs.

For principal component and admixture analyses [9], we further pruned SNPs in high linkage disequilibrium (plink --indep-pairwise 200 25 0.4), resulting in 120,384 autosomal SNPs. Principal component analysis was done using the smartpca program included in EIGENSOFT [23,24] with default settings. Admixture fractions were inferred using ADMIXTURE [25], with 2–15 clusters (K = 2–15) and 50 replicates, and visual inspection was performed with PONG [26].

Runs of homozygosity were detected in the study participants from the current work. We measured the lengths of homozygosity with PLINK using the following parameters: --homozyg --homozyg-window-kb 5000 --homozyg-window-het 1 --homozyg-window-threshold 0.05 --homozyg-kb 500 --homozyg-snp 25 --homozyg-density 50 --homozyg- gap 100. We summed the total length of homozygous blocks (binned by fragment sizes) for each individual, and then averaged the results by population.

2.4 Selection scans

2.5 Gene annotation and overlapping signals

Selection scan signals were identified by averaging the –log₁₀(p-value) of iHS, as well as the XP-EHH and LSBL values, on a sliding window of 10 SNPs. We first evaluated selection scan signals common to all Andean study populations by annotating the top 0.5% of selected SNPs obtained with each positive selection method and population to genes using Variant Effect Predictor [32]. This 0.5% cut-off is a conservative approach compared to that used in other Andean studies in the literature, which used 1% or 5% cut-offs because their genotype data was more sparse than the data used in the current study [33,34]. We searched for genes that were within 10 kb on either side of the input SNP to investigate potential distant regulatory effects of SNPs. A gene list was obtained for each selection scan method (n = 3) and population of interest (n= 3), resulting in 9 gene lists in total. Genes overlapping across populations and selection scan methods were identified with GeneVenn [34]. We analyzed the overlap across populations based on genes instead of SNPs since we used the averaged selection scans to minimize single-SNP outliers.

We then evaluated traits that had undergone positive selection at a population level by selecting the top five selection regions for each population and method. From each of these regions, we selected the top SNP and identified the genes within 200 kb of either side of the SNP using Genome Browser.

2.6 Pathway enrichment analyses

We performed pathway enrichment analyses to determine the biological implications of the putatively selected genes shared across the three Andean study groups (Aymara-Quechua, Uru, and SAC). First, we used the Ingenuity Pathway Analysis software (content version 68752261; Ingenuity Systems, USA), using Fisher’s exact test. We did not correct for multiple comparisons in the IPA pathway analysis as the genes selected were based on 0.5% of SNPs obtained with each positive selection method and population, which is a conservative cut-off. Then, we carried out over-representation analyses (ORA) of the list of putatively selected genes across Andean study populations with the WebGestalt platform [35], using the KEGG and PANTHER Pathway databases. Here, we corrected for multiple testing to be stricter than the more explorative analysis by using IPA. WebGestalt uses the FDR method of Benjamini and Hochberg [36]. Pathway enrichment analysis is relevant when the study genes have cell signaling functions and therefore likely participate in several molecular and biological pathways. However, it does not necessarily pinpoint the selective pressure driving such positive selection sweeps in the genome.

3.1 Population structure and ancestry

We previously evaluated the population structure of the Bolivian and Argentinean study groups [9]. For this, we assessed the genetic affinities of the three Andean study groups in comparison to indigenous populations from North and South America [20,21], Europe, and Africa (to account for admixture). The Andean populations clustered together in principal component analysis (Figure S1A). Furthermore, in the population structure analysis, the ancestry fraction at K = 3 for European origin was 5.6% for SAC, 1.9% for Aymara-Quechua, and 0.2% for Uru, indicating limited European admixture in the three Andean study groups (Figure S1B).

Groups with a history of isolation often show high levels of homozygosity across their genomes. By assessing homozygosity for the three study groups with comparative data from the Americas (as assembled previously [21]; Figure S2), we found that the Andean study groups had levels of extended homozygosity comparable to other Native American populations, and in particular to groups that have limited non-Native American admixture. We also observed a low number of long homozygous regions (>8 Mb) in the Andean groups (including Uru), implying limited recent endogamy.

3.2 Shared signals of positive selection

3.2.1 Selective signals across Andean populations

To study the extent and nature of the shared signals of positive selection between three populations from the Andes we used the lists of variants and genes by three selection tests reported in De Loma et al. 2022 study: two based on haplotype length (integrated haplotype score [iHS] and cross-population extended haplotype homozygosity [XP-EHH]) and one based on pairwise population differentiation at each given locus (locus-specific branch length [LSBL]). To minimize the detection of spurious single nucleotide polymorphism (SNP) outliers, we calculated the average p-values for iHS (Figure 1), and the average XP-EHH and LSBL values across the whole genome by sliding-window averaging (Figure S3 and Figure S4). For each selection scan method and population, we then focused on the top 0.5% candidate SNPs. Subsequently, we identified genes within 10 kb of either side from each top SNP. The number of putatively selected genes detected in the SAC, Aymara-Quechua, and Uru populations were 436, 507, and 454 for iHS; 340, 343, and 312 for XP-EHH; and 718, 720, and 882 for LSBL, respectively.

Figure 1. Genome-wide iHS selection scans for the three Andean study populations. To avoid spurious results due to single SNPs, we averaged the –log₁₀(p-value) of iHS with a sliding window of 10 SNPs. The top five peaks in each population are marked with red lines, as well as the locations of OPRM1, CSMD1, and GSG1L. The top 0.5% (solid line) and top 0.01% (dashed line) are shown.

We evaluated adaptation-related traits common to all three Andean populations by identifying putatively selected genes (from any of the three selection scan methods) that overlapped across the three populations. This approach accounts for the fact that different selection scan methods are based on different principles, and therefore not all signals of positive selection are expected to be identified by all three methods. A total of 116 putatively selected genes were shared across all three Andean study groups (Figure 2A, Table S1). Among genes shared between two of the three groups, more genes were shared between the Aymara-Quechua and SAC populations (n = 307) than between either of those and the Uru population (n = 151 shared with SAC, n = 160 shared with Aymara-Quechua; Figure 2A, Table S1). For example, the arsenic methylating gene AS3MT was shared between the Aymara-Quechua and Uru populations, and the opioid receptor gene OPRM1 was shared between the SAC and Aymara-Quechua groups (Figure 3A).

Figure 2. Genes potentially under positive selection identified by the three independent selection scan methods in the three Andean study populations (SAC, San Antonio de los Cobres; AQ, Aymara-Quechua; U, Uru). Genes included in the Venn diagrams are those annotated as being within ±10 kb from the top 0.5% SNPs identified by each selection scan method and in each population. (A) Venn diagram of putatively selected genes (identified by any of the selection scan methods) that are common among the three Andean study populations. (B) Venn diagrams per population of putatively selected genes shared across selection scan methods.

Figure 3. Zoom-ins of selection scans for OPR1M (chromosome 6), CSMD1 (chromosome 8), and GSG1L (chromosome 16) in the Aymara-Quechua study population (n = 65). The width of the red line indicates the size of each gene. The most representative selection scan method for each gene was included. The genome-wide average (solid line) and top 0.5% (dashed line) are shown.

To identify whether the putatively selected genes shared across the Andean study populations were over-represented in specific biological pathways and/or functions, we performed gene-enrichment pathway analyses (Table 1). The Ingenuity Pathway Analysis (IPA) software mapped 90 of the 116 shared gene IDs to specific loci in the IPA database; 26 gene IDs that encoded non-coding RNAs were not in the IPA database (Table S2). IPA canonical pathway analysis showed significant enrichment for the sperm motility pathway, as well as multiple pathways related to immune defense against pathogens, cardiovascular functions (including the renin-angiotensin pathway, although not statistically significant), and neurological signaling (Table 1, Table S3). In addition, there was an enrichment for the opioid signaling pathway (Table 1, Table S3). In the disease or function analysis, 71 out of 90 shared genes were associated with cancer (Table S4).

Table 1. Gene-enrichment analyses for putatively selected genes shared among Andean populations.

Putatively selected genes were those at ± 10 kb from the top 0.5% SNPs in each selection scan method. Enrichment pathway analyses were performed on the 116 genes that were shared among the Andean study populations according to any of the three selection scan methods. Two different programs were used: Ingenuity Pathway Analysis (IPA) and WebGestalt. The number of mapped genes in the datasets was 90 for IPA, 21 for KEGG, and 7 for PANTHER. Rank is within each dataset.

We performed complementary over-representation analyses using the web-based application WebGestalt and the KEGG and PANTHER Pathway databases. The morphine addiction and renin secretion pathways in KEGG, although not statistically significant after adjustment for false discovery rate (FDR, p-value < 0.05), were enriched among the shared putatively selected genes (Table 1, Table S5). In the PANTHER Pathway analysis, the histamine H1 receptor-mediated signaling pathway, related to inflammation, was enriched, along with the muscarinic acetylcholine receptor 1 and 2 signaling pathway (Table 1, Table S6). Furthermore, the endothelin signaling pathway was included within the top 10 results of the PANTHER Pathway analysis (Table S6).

3.3 Signals of positive selection per method and population

We next investigated the top five regions identified by each selection scan method in each population (Figure 1, Figure S3, and Figure S4 and Tables S10–S12). Since this is the approach most commonly used in adaptation studies in the Andes, where fewer methods and populations have been analyzed simultaneously, we included this approach to facilitate comparison with the literature deriving from this geographic area. Several genes in the top five selection regions were previously identified in selection studies on adaptation to high altitude, for example CPNE4, CSMD1, FOXB1, IL17F, and SP100, and in studies related to arsenic adaptation, such as AS3MT and LHFP (Tables S10–S12). In addition, several of the genes (n = 14, for example CPNE4, CSMD1, LHFP, SLC26A3, and TSPAN5) within the top selection peaks were present in the list of 116 putatively selected candidate genes shared across all Andean study populations (marked in bold in Tables S10–S12). Furthermore, two protamine genes (PRM1 and PRM2), potentially related to sperm motility, that scored highest in the IPA enrichment analysis, were the third highest XP-EHH signal in the SAC study group.

4.1 Opioid signaling

We identified signals of positive selection around OPRM1 and other genes related to opioid signaling (ITPR2, KCNJ5, and PRKCH) and morphine addiction (KCNJ5, PDE3A, and PDE4D). In addition, we found enrichment for the muscarinic acetylcholine receptor signaling pathway, which is linked to the opioid signaling pathway [41]. The opioid signaling pathway and its receptors interact with drugs such as morphine, methadone, and cocaine, and influence dependence on these substances [42]. Rats given cocaine showed a dose-response induction of OPRM1 in the brain (Soderman and Unterwald 2009) [43]. Genetic variants of OPRM1 have been linked to cocaine, heroin, and alcohol consumption in populations of European ancestry [44,45]. Interestingly, genetic variants of GSG1L (the only putatively selected gene in all three selection scan methods in both Bolivian populations) are associated with the physiological response to methadone in heroin-dependent patients [46].

Chewing coca leaves, which contain cocaine, is a well-known and longstanding practice in the Andean region as part of rituals and to prevent or mitigate altitude sickness, hunger, fatigue, and pain [47]. Traces of cocaine in mummies from Chile suggest that this practice dates back at least 3,000 years [47]. Opioids and cocaine have a central role in regulating pain, and sequence variants of OPRM1 seem to modify how humans handle pain [48]. A study in Peru highlighted that children with altitude-induced chronic hypoxia, living in Cusco (3399 m above sea level) and mostly from Andean indigenous communities, required 40% less opioids for analgesia than those living in Lima at sea level [49]. In this context, our identification of putatively selected genes related to opioid signaling potentially implicates an ethnicity-specific mechanism of pain perception, driven by adaptation to the use of coca leaves to overcome altitude sickness, hunger, fatigue, and pain.

Selection of genes in opioid-related pathways have been reported before, but not related to adaption to a selection pressure, such as cope with living in the Andes. Rockman et al. (2005) found signals of selection at the human prodynorphin gene, which is an opioid polypeptide hormone [50].

4.2 Immune defense against pathogens

At least four signaling pathways related to the immune system were significantly enriched within the 116 putatively selected genes common across the three Andean populations. Furthermore, our complementary enrichment analysis identified genes involved in the histamine H1 receptor signaling pathway. The histamine H1 receptor participates in the immune response, and its antagonists are regularly used against allergies. The histamine H1 receptor expression is modulated by Mycobacterium tuberculosis, which causes tuberculosis [51]. Recently, a study on indigenous people living in the Ecuadorian highlands showed positive selection for immune function related to tuberculosis [52]. Other antagonists of the histamine H1 receptor are effective against parasitic infections, such as leishmaniasis [53] and Chagas diseases [54].

Another gene that was putatively selected across all Andean populations and is involved in the defense against pathogens was Mucin 17 (MUC17). Mucins are involved in the immune response and host–parasite interactions [55]. Previous selection studies in humans also identified mucin genes, such as MUC19, in a population from central Mexico [56]. Positive selection of mucin genes could lead to resistance to infection by parasites and other pathogens.

4.3 Reproductive function

The top pathway enriched in the IPA analysis was the sperm motility pathway. In extreme environmental conditions, such as in the Andes with high altitude, harsh temperatures, and water scarcity, reproductive fitness is crucial. This is reflected in the selection scans of the SAC study group, where we identified within the top selection peaks signals near PRM1 and PRM2, which encode protamines involved in the condensation of chromatin during sperm formation [57]. Protamine genes are rapidly evolving in humans because of positive selection [58,59]. Although the evolution of protamine genes has been widely studied across species, few have reported positive selection in human indigenous populations in relation to specific environmental factors. A study including the native American populations from the Human Origins project [21], identified signs of positive selection in the region including the protamine genes [60].

In addition, CSMD1 was identified through all selection scan methods in the Aymara-Quechua and Uru populations, as well as in the LSBL and iHS scans in the SAC group. Although its most characterized function is related to inflammation, rare mutations in CSMD1 have been linked to infertility in males and females [61]. Recently, higher reproductive fitness of native Tibetans at high altitude was linked to polygenic adaptation [62].

4.4 Arsenic adaptation

We previously identified signals of positive selection near AS3MT, which encodes the main methyltransferase enzyme that metabolizes arsenic in humans, linked to efficient arsenic metabolism in the study groups from the Argentinean Andes [10] and the Bolivian Andes [9]. In the current study, despite having fewer individuals in each study group, AS3MT was still within the top selection signal for the LSBL method in the Aymara-Quechua population (Table S12), and was shared between the Aymara-Quechua and Uru in all of the methods (Table S1). Since the first description of arsenic adaptation in humans [10], selection for AS3MT has been identified in study groups from Argentina, Chile, Peru, and Bolivia [9,63-65]. Together, these results suggest that there has been strong selective pressure for metabolizing arsenic in several groups across the Andean Mountain range. However, the causal variants driving this adaptation to arsenic have not been identified but a recent study suggests multiple potential causal variants, some ethnically distinct, around AS3MT that are influencing arsenic metabolism [66]. It should be stressed that even if these populations have a genetically determined efficient arsenic metabolism, they are not immune to arsenic toxicity, and further efforts should be made to lower arsenic exposure in these areas.

4.5 Altitude adaptation

Current evidence implicates multiple genes in the adaptation to high altitude and indicates that these genes differ depending on the study population and the region [67]. In our study, we found two genes potentially associated with high-altitude adaptation in the Andes: CSMD1 and ALDH5A1 (Tables S10–S12 show all genes previously linked to altitude adaptation within the top selection peaks). CSMD1 participates in the complement pathway involved in inflammation, and it is considered a tumor suppressor whose allelic loss is associated with poor prognosis in epithelial cancers [37]. This gene has also been detected in selection scans for altitude adaptation in Tibetans, although it was not among the top candidates in that population [68]. Still, how it could confer adaptation to high altitude is unknown.

ALDH5A1, a candidate selected gene shared among all populations in our study, encodes an aldehyde dehydrogenase that catalyzes a step in the degradation of the neurotransmitter gamma-aminobutyric acid (GABA). Deficiency of this enzyme is known as 4-hydroxibutyricaciduria, a rare autosomal metabolic disease that leads to an abnormal accumulation of 4-hydroxibutyric acid (GHB) [69]. Under hypoxia, ischemia, or excessive metabolic demands, GHB levels rise in order to protect central and peripheral tissues [70,71]. Since GHB is considered to protect against hypoxic conditions, increased ALDH5A1 expression could be beneficial during altitude-induced hypoxia.

Several putatively selected genes were enriched in pathways related to renin-angiotensin signaling (ITPR2 and PRKCH) and renin secretion (ITPR2 and PDE3A), which could potentially be connected to high altitude adaptation. Changes in fluid retention and regulation of blood pressure are physiological processes involved in high altitude acclimatization. Indeed, genetic variants of the angiotensin-converting enzyme (ACE) have been associated with adaptation to high altitude [72]. Interestingly, adaptation to high altitude in the Andes is linked to increased hemoglobin levels in blood, compared to high-altitude Tibetan populations where this phenotype is not seen to the same extent [73]. High levels of hemoglobin in blood are a risk factor for cardiovascular disease [74,75]. Therefore, we speculate that potentially selected pathways related to cardiovascular functions may reflect a need to compensate for the higher hemoglobin levels induced to cope with hypoxia in the Andes at high altitude.

We did not find signal of selection around the HIF2A locus [7] or other genes previously linked to high altitude adaptation. It should be noted that some of the PEL population living at low altitude in Lima may be recent migrants from high altitude and therefore altitude adaptation signals may not be that strong. Nonetheless, in the present study we still identify several genes linked to altitude adaptation, suggesting that using PEL as the reference group provides adequate power to identify altitude adaptation.

4.6 Cancer

The enrichment analysis identified 71 out of 90 putatively selected genes as associated with different types of cancer. A recent study of Native American ancestry and cancer mortality risk in Chile showed that a greater degree of Aymara descent is significantly associated with increased mortality due to skin, bladder, larynx, bronchus, and lung cancer [76]. Populations in the highlands of Ecuador also have an elevated incidence of several types of cancer, including colon, hematopoietic, and liver cancers, supporting the possibility of a higher cancer risk among people living at high altitude [77]. A recent metanalysis found that genes selected for adaptation to extreme environments, such as high altitude or cold, are enriched for cancer-associated genes [78]. These previously reported cancer patterns differ from the pattern of top cancer types we found in our enrichment analysis. However, when interpreting the findings for cancer and cancer profiles, it is worthwhile considering the current over-representation of cancer data in the literature and the existing bias for certain cancer forms. Alternatively, the results reflect that genetic variants selected to enhance survival under extreme conditions could be associated with pleiotropic effects including promoting cancer risk.

In conclusion, we identified novel selection signatures related to opioid signaling, sperm motility, and altitude adaptation common among Andean populations in the Bolivian Altiplano and the Argentinean Puna. There are a longstanding custom of chewing coca leaves in the Andes to mitigate altitude sickness, hunger, fatigue, and pain. Considering this, selection for opioid signaling could conceivably represent a mechanism for coping with stressors during human settlement in the harsh environment of the Andes. This potential connection between modified opioid signaling and coca-leaf chewing provides evidence of possible gene–cultural co-evolution in adaptation to high altitude. Our positive selection scans also found that the reproductive fitness of these populations, immune defense against pathogens, and high-altitude adaptation may have been driven by genetic factors.

Ethics Statement

Ethical approval for this study was obtained by the Comisión Provincial de Ética del Ministerio de Salud de Salta (Argentina), the Comité Nacional de Bioética (Bolivia), and the Regional Ethic Committee of Karolinska Institutet (Sweden, no. 2020-00495 and 2020-00493).

Consent for Publication

Informed consent for the publication of this study from the concerned individual has been obtained.

Availability of Data and Material

The genome-wide genotypes presented in this paper will be made available under controlled access via Federated European Genome-phenome Archive (FEGA) Sweden once this repository becomes operational. The datasets will then be findable through the European Genome-phenome Archive web portal (https://ega-archive.org). Until the dataset included in this study is available at SciLifeLab Data Repository (https://figshare.scilifelab.se/articles/dataset/_b_Human_adaptation_in_the_Andes_Mountains_b_/25323256)

Funding

This work was supported by the Institut de Recherche pour le Développement, Bolivia; the Hydrosciences Montpellier Laboratory, France; the Eric Philip Sörensens Stiftelse, Sweden; the Swedish Research Council, Sweden (2017-01239 and 2020-02827 to K. Broberg, and 2014-5211 to C. Schlebusch); Karolinska Institutet, Sweden; the Swedish International Development Cooperation Agency, Sweden (2016-05799 and 2019-03832 to K. Broberg); and the Genetics Institute at the Universidad Mayor de San Andrés, Bolivia.

Competing Interests

Carina Schlebusch is a member of the Editorial Board of the journal Human Population Genetics and Genomics. The author was not involved in the journal’s review of or decisions related to this manuscript. The authors have declared that no other competing interests exist.

Author Contributions

Conceptualization: K.B. and C.S.; Methodology: J.D.L., M.V. C.S., K.B.; Software: M.V.; Formal analysis: J.D.L and M.V.; Investigation: J.D.L. J.G. N.T. K.B., F.A. L.P.; Writing – Original Draft: J.D.L.; Writing – Review & Editing: J.D.L. J.G. N.T. K.B., F.A. L.P. C.S., M.V.; Visualisation: J.D.L., M.V.; Supervision: K.B., C.S.; Funding acquisition: K.B., C.S., J.G., N.T.

Acknowledgement

We warmly thank all the participants who took part in this study. We would also like to thank the Servicio Departamental de Salud de Oruro (SEDES) and the Salud Familiar Comunitaria e Intercultural (SAFCI) program for their contributions during the field activities in Bolivia. We also thank the fieldwork teams and the workers at local health centers in San Antonio de los Cobres and in the villages around Lake Poopó for assisting during recruitment. We also thank Gabriela Concha for her invaluable contribution during the recruitment in San Antonio de los Cobres. We thank Marta Toscano and Paola Zago for helping with the informed consents in San Antonio de Los Cobres. We acknowledge Anda Gliga’s help regarding the pathway enrichment analyses. We thank Diana Ekman at the National Bioinformatics Infrastructure Sweden at SciLifeLab for productive discussions about SNP-to-gene location. Also, we are thankful to Chiara Barbieri for sharing her data prior to publication.

The genotyping was performed by the SNP&SEQ Technology Platform in Uppsala (www.genotyping.se). The facility is part of the National Genomics Infrastructure supported by the Swedish Research Council for Infrastructures and Science for Life Laboratory, Sweden. The SNP&SEQ Technology Platform is also supported by the Knut and Alice Wallenberg Foundation. The computations were performed at the Swedish National Infrastructure for Computing (SNIC-UPPMAX; Projects SNIC 2018/8-107 and SNIC 2021/22-262).