Genetic Signatures of Positive Selection in Human Populations Adapted to High Altitude in Papua New Guinea

Abstract

Papua New Guinea (PNG) hosts distinct environments mainly represented by the ecoregions of the Highlands and Lowlands that display increased altitude and a predominance of pathogens, respectively. Since its initial peopling approximately 50,000 years ago, inhabitants of these ecoregions might have differentially adapted to the environmental pressures exerted by each of them. However, the genetic basis of adaptation in populations from these areas remains understudied. Here, we investigated signals of positive selection in 62 highlanders and 43 lowlanders across 14 locations in the main island of PNG using whole-genome genotype data from the Oceanian Genome Variation Project (OGVP) and searched for signals of positive selection through population differentiation and haplotype-based selection scans. Additionally, we performed archaic ancestry estimation to detect selection signals in highlanders within introgressed regions of the genome. Among highland populations we identified candidate genes representing known biomarkers for mountain sickness (SAA4, SAA1, PRDX1, LDHA) as well as candidate genes of the Notch signaling pathway (PSEN1, NUMB, RBPJ, MAML3), a novel proposed pathway for high altitude adaptation in multiple organisms. We also identified candidate genes involved in oxidative stress, inflammation, and angiogenesis, processes inducible by hypoxia, as well as in components of the eye lens and the immune response. In contrast, candidate genes in the lowlands are mainly related to the immune response (HLA-DQB1, HLA-DQA2, TAAR6, TAAR9, TAAR8, RNASE4, RNASE6, ANG). Moreover, we find two candidate regions to be also enriched with archaic introgressed segments, suggesting that archaic admixture has played a role in the local adaptation of PNG populations.

Introduction

The peopling of the island of New Guinea, located in the Pacific directly north of Australia, represents one of the earliest human migrations out of Africa, with the first settlers arriving 50 thousand years ago (kya) (O’Connell et al. 2018). Upon settling, they rapidly spread across different regions, with evidence of ancient occupation in both the highlands and lowlands dating back to ∼40 kya (Summerhayes et al. 2010). This early dispersal eventually led to genetic divergence between the inhabitants of the Highlands and Lowlands occurring approximately 10 to 20 kya (Bergström et al. 2017; Brucato et al. 2021). Today, Papua New Guinea (PNG) is renowned for its exceptional cultural and linguistic diversity, making it one of the most varied areas globally. This diversity was further enriched by subsequent migrations that included the arrival of ancestral Austronesian-speaking people around 3.5 kya, who primarily settled along the coastal regions (Kayser 2010; Skoglund et al. 2016). The latter part of the 19th century saw an increase in European settlers arriving to the island, also mainly in coastal areas, while populations in the Highlands remained relatively isolated from European contact until the early-mid 20th century. The Highlands and Lowlands of PNG are very distinct geographical areas. For example, the Highlands feature rugged terrain with some of the highest mountains in the Pacific, ranging between ∼1,200 and 4,500 meters above sea level, and cooler temperatures, whereas the Lowlands are characterized by a more humid tropical environment with higher temperatures and denser rainforests. Such environments have influenced past settlement patterns, with malaria likely playing a significant role in driving populations to settle in the Highlands, away from the mosquito-dense Lowlands where the spread of malaria is facilitated due to ideal conditions for the Plasmodium-carrying Anopheles mosquito to thrive (Senn et al. 2010; Trájer 2022). It is possible, but uncertain, that other prevalent infectious diseases, such as tuberculosis (TB), were present to a degree before the first European contacts in the early 16th century, potentially affecting indigenous populations to some extent long before wider global interactions (Riley 1983; Comas et al. 2013; Trájer 2022). Importantly, historical reports indicate that the arrival of Europeans favored the spread of new and previously localized diseases across the island (Ley et al. 2014). Differences between the Highlands and Lowlands ecoregions, and the history of relative population isolation across PNG, might have led to different selective pressures acting in each one. Previous research has found that high altitude is a major factor contributing to the phenotypic variation observed between highlanders and lowlanders in PNG, suggesting this to be the result of an adaptive response to the selective pressure exerted by hypoxia and highlighting the need for further studies to investigate the genetic basis of such adaptation (André et al. 2021). A more recent study on two populations living at distinct altitude levels in PNG found candidate genetic signals related to cardiovascular phenotypes: in highlanders related to red blood cell count, and in lowlanders related to white blood cell count (André et al. 2024). Additionally, long-term exposure to high altitude can induce chronic mountain sickness (CMS), angiogenesis, increased oxidative stress, and cardiovascular complications (Pena et al. 2022). Moreover, populations from New Guinea have a history of genetic admixture with archaic hominins, exhibiting some of the highest proportions of Denisovan introgression (∼4%) (Reich et al. 2010; Meyer et al. 2012; Vernot et al. 2016; Jacobs et al. 2019), and levels of Neanderthal introgression similar to those of other populations outside of Africa (∼2%) (Prüfer et al. 2014). Archaic introgression has been identified to have adaptive effects on human populations (Racimo et al. 2015; Gittelman et al. 2016). For instance, Neanderthal introgression has contributed adaptive traits related to the immune response in modern humans (Enard and Petrov 2018), whereas Denisovan introgression has contributed to traits related to high altitude adaptation in Tibetan populations (Huerta-Sánchez et al. 2014). Recent studies have shown that Denisovan introgression might have contributed to the immune system of populations from PNG (Vespasiani et al. 2022; André et al. 2024; Yermakovich et al. 2024), and that some signals of selection found in Pacific Islanders are located in Denisovan segments harbored by Papuan ancestry (Choin et al. 2021; Brucato et al. 2022). André et al further elucidated the role of archaic introgression in local adaptive processes among PNG populations, reporting an archaic haplotype potentially related to angiogenesis in highlanders, and another in lowlanders related to immunity (André et al. 2024). Whereas Yermakovich et al explored Denisova introgression and its potential adaptive roles in brain biology and immune response in PNG populations (Yermakovich et al. 2024). Despite these advancements, the full scope of adaptive archaic introgression continues to be an area for further investigation. The unique setting and history of PNG provide an opportunity to discover regions in the genome underlying human adaptation to contrasting environment conditions in the Highlands and the Lowlands, as well as to assess the involvement of archaic introgression in molecular signatures of selection.

To address this, we investigated signals of positive selection in 62 highlanders and 43 lowlanders sampled from 14 different locations across the main island of PNG using whole-genome genotype data gathered from the Oceanian Genome Variation Project (OGVP) (Quinto-Cortés et al. 2024) (Fig. 1a). Our selection scans included statistics suitable to detect both within and between population signals, followed by functional annotation clustering of candidate genes to identify potential biological processes under selection in the Highlands and Lowlands, and archaic ancestry estimation to detect signals possibly involving introgressed regions of the genome.

Results

Population Structure and Admixture

To explore patterns of global ancestry and population structure across the Oceanian region we performed Principal Component Analysis (PCA) and ADMIXTURE runs combining our PNG samples and the broader OGVP data set Quinto-Cortés et al. (2024) as well as the Human Genome Diversity Project (HGDP) Han, French, and Yoruba from (Bergström et al. 2020) as continental reference samples (supplementary fig. S1, Supplementary Material online). We then focused on populations from the main island of PNG for downstream analyses (Fig. 1b and c). Principal component 1 (PC1) separates African from non-African ancestries, whereas principal component 2 (PC2) separates PNG populations from African, European, and East Asian ancestries. Individuals from the Highlands form a tighter cluster than those from the Lowlands, which are more dispersed in PCA space, ranging from near the highlanders to the axis toward the Han cluster, reflecting the presence of other ancestries besides Papuan-related that are closer to the East Asian ancestry of the Han (Fig. 1b). In ADMIXTURE (Fig. 1c), we found the lowest cross-validation error at K = 9 (supplementary fig. S2, Supplementary Material online), identifying the following ancestry components: African, European, and East Asian, along with Papuan-related ancestry, and five other distinct ancestries mainly found in offshore PNG islands and other islands across the Bismarck Archipelago, the Solomon Islands, Vanuatu, and Fiji in the Southwest Pacific (Figure S3 in Quinto-Cortés et al. 2024). Individuals from the Lowlands show varying proportions of these Pacific-related ancestry components, consistent with other findings from the region (Choin et al. 2021), whereas individuals from the Highlands show a predominant Papuan-related ancestry. No European or African ancestries were observed in populations from the Highlands or Lowlands.

Candidate Regions Under Positive Selection

Since different environmental pressures might exist between the Highlands and Lowlands (oxygen availability, temperature variability, ultraviolet radiation, pathogen/s exposure, staple food types, etc.), we categorized samples from each province into two ecoregions (supplementary table S1, Supplementary Material online), and conducted four different tests for positive selection: the integrated haplotype score (iHS) (Voight et al. 2006), the number of segregating sites by length (nSL) (Ferrer-Admetlla et al. 2014), the cross-population number of segregating sites by length (XP-nSL) (Szpiech et al. 2021) and the population branch statistic (PBS) (Yi et al. 2010).

We then integrated results from compatible statistics into two Fisher's Combined Scores (FCS). We combined the ranked results of iHS and nSL as FCS_within, comprising signals obtained from within-population tests, and the ranked results of XP-nSL and PBS as FCS_between, comprising signals obtained from between-population tests (see Materials and Methods). While the FCS_within represents the integration of signals found in a given ecoregion without a reference population to contrast against, the FCS_between represents the integration of signals identified in a given ecoregion when contrasting against the other. In the Highlands, FCS_between denotes signals obtained from comparing against the Lowlands and conversely, FCS_between in the Lowlands represents signals obtained from comparing against the Highlands.

We initially identified the top 50 candidate genomic regions putatively under positive selection in each ecoregion, per FCS computation, using an empirical outlier approach (see “Analyses of Genomic Signatures of Positive Selection” in Materials and Methods) (Fig. 2). This yielded two separate sets of candidate regions: 50 regions identified through FCS_within and another 50 through FCS_between, accounting for a total of 100 candidate regions in each ecoregion. Subsequently, we annotated all genes contained in these candidate regions using ANNOVAR (Wang et al. 2010) and noncoding variants in regulatory elements to target genes using GeneHancer (Fishilevich et al. 2017). Additionally, we used UniProt (The UniProt Consortium 2017) and GeneCards (Stelzer et al. 2016) to further characterize gene functions (supplementary table S2 and S3, Supplementary Material online). We considered a candidate region under positive selection to be unique to an ecoregion if it did not overlap in 50 kb or more with a candidate region from the other ecoregion, or if it was not entirely contained within another from the other ecoregion. The top 10 unique candidate regions in each ecoregion, per FCS computation, are described in Table 1 for the Highlands, and Table 2 for the Lowlands. We also identified shared candidate regions between ecoregions among the initial top 50 signals, per FCS computation, and separated them for individual analysis as a distinct group (supplementary table S4, Supplementary Material online).

In each ecoregion, we identified the top 50 strongest candidate regions based on FCS computations, amounting to 100 candidate regions per ecoregion (50 from FCS_within and 50 from FCS_between for both Highlands and Lowlands). Out of these 200 signals (100 in each ecoregion), we found 23 to be common between the two ecoregions: 11 originating from the Highlands and 12 from the Lowlands. In these shared candidate regions, two were ranked among the top 10 in highlanders, and four among the top 10 in lowlanders (Table 3). After excluding these 23 shared candidate regions from the initial set of 100 in each ecoregion, we analyzed the remaining candidate regions, resulting in 41 regions from FCS_within and 48 from FCS_between in the Highlands (supplementary table S2, Supplementary Material online), and 40 from FCS_within and 48 from FCS_between in the Lowlands (supplementary table S3, Supplementary Material online).

Candidate Genes Under Selection in the Highlands

When analyzing the top 10 candidate regions under selection private to the Highlands (Table 1), we observed that both FCS approaches, FCS_within and FCS_between, identified genes associated with similar biological processes and pathways, albeit pinpointing different specific genes within these categories. The genes involved in glucose metabolism, PCSK2 and DIRC3 were identified by FCS_within, while FCS_between identified LDHA, known to not only play a role in glucose metabolism but also in the Hypoxia-Inducible Factor (HIF) pathway (Taylor and Scholz 2022). Similarly, genes participating in ion homeostasis and transport were detected (WNK2, KCNK16, KCNK17, and SLC13A4 by FCS_within, and LRRC38 by FCS_between). Notably, FCS approaches identified genes related to the Wnt pathway, a major developmental pathway that plays a role in hypoxia tolerance (Gersten et al. 2014) (DAAM2, MCC, LINC01119 identified by FCS_within and PNMA1 identified by FCS_between).

Genes in the top 10 candidate regions from the Highlands identified only by the FCS_within approach are related to cell structure and transport (KIF6, NUP205), neurotransmission and hormone release (CADPS), and functions related to the regulation of the nervous system and cell cycle control (BRINP2, STMP1). Through this approach, we also found MCFD2, a gene with a role in blood coagulation, and TTC7A, involved in the development of the intestinal epithelial and immune system.

The FCS approach designed to identify candidate regions under selection in the Highlands when compared against the Lowlands, FCS_between, identified genes linked to various biological processes, some of which could have implications for high-altitude adaptation. This is the case of candidate genes involved in oxidative stress and redox regulation (MUTYH, PRDX1, ASPHD2, HEATR4), heart and vascular cell function (MIAT, MIATNB, DPF3), Notch signaling (PSEN1, NUMB), fatty acid metabolism (ACOT2, ACOT4, ACOT6) and inflammatory response (SAA2-SAA4, SAA4, SAA2, SAA1). Other signals in the top 10 candidate regions are related to the immune system (PIK3R3, ZFYVE1, PDPN), epigenetic modifications and chromatin structure (PHC3, RIOX1, ELMSAN1, PRDM2, RPSAP52, HMGA2, GTF2H1), retina function and eye development (SAMD7, CRYBB1, CRYBA4, SLC7A14), translational modifications (TPST2, DCAF4) as well as post-transcriptional gene regulation (TFIP11, RBM25), cell migration and signal transduction (SEC62, GPR160, CLDN11), pigmentation (HPS4, HPS5), circadian rhythm (SRRD), and ciliary function (DNAL1).

Candidate Genes Under Selection in the Lowlands

Similarly, we analyzed the top 10 candidate regions with signatures of selection private to the Lowlands (Table 3) and found that both FCS approaches revealed genes linked to analogous biological processes and pathways. However, as in the Highlands, these signals pinpointed distinct genes within the identified categories. The strongest signals were detected in genes related to immune response: FYB1, C9, and OASL identified through FCS_within, and HLA-DQB1, HLA-DQA2, SIRPA, ST3GAL1, BLNK, and DNTT identified when comparing Lowlands against Highlands with FCS_between. Other genes were also identified to be associated in neuronal functions and ion transport: NKAIN3, SCN3A, SCN2A (FCS_within), GPRIN3, and OPALIN (FCS_between), as well as genes involved in spermatogenesis and sperm function (C4orf45 obtained with FCS_within and ARRDC5 obtained with FCS_between). Other biological processes identified through both approaches include genes playing a role in histone modification and gene expression (HNF1A, FNIP2 with FCS_within and UHRF1, KDM4B with FCS_between). As well as cell adhesion and migration (IQSEC1 with FCS_within and MIR1246, LINC01116 with FCS_between).

When looking at genes identified only in the between population comparison of Lowlands against Highlands (FCS_between), we found MIR3135B and PDYN-AS1 to be top candidates, although there is currently no assigned function or potential relationship to any biological process in Uniprot and Genecards for these two genes.

Regarding genes revealed only through the FCS_within approach, we found candidates involved in folate metabolism (GGH), vitamin E transport (TTPA), solute transport and cellular metabolism (SLC38A11), skeletal development (TLL2), RNA metabolism and modification recognition (YTHDF3), and serum urate control (LINC01229).

One of the strongest signals detected in both ecoregions with the FCS_within approach is the region containing genes KAT6B, DUPD1, DUSP13, and SAMD8. In the Highlands, these four genes are part of the highest-ranking candidate region, resulting from the FCS_within approach, which comprises 11 candidate genes in total. These four genes are also present in 7th strongest candidate region found in the Lowlands using the same FCS approach.

Among the 10 strongest shared signals in overlapping candidates under selection in each ecoregion, we found genes involved in gene regulation and chromatin structure (KAT6B, ZBED9, ZNF311, ZSCAN12, ZSCAN23), cell signaling and post translational modifications (DUPD1, DUSP13, TRIM27, UBD). Some also constitute olfactory receptors (OR2B3, OR2J1, OR2J3), antioxidant enzymes (GPX5 and GPX6), and some others are related to the immune response (HCG14, HCG15), lipid metabolism (SAMD8), and RNA processing (SNORD32B). We also found shared uncharacterized long noncoding RNAs (LINC00533, LINC01623, LINC01556, LOC100129636, LOC101929006, LOC101928516, and LINC02147).

Haplotype Structure in Candidate Regions for Positive Selection in the Highlands and Lowlands

To better understand the haplotype structure of the top 10 candidate regions under selection when comparing Highlands and Lowlands (FCS_between), we constructed haplotype networks and determined the decay of site-specific extended haplotype homozygosity (EHHS) between ecoregions around each candidate region using the R packages pegas (Paradis 2010) and rehh (Gautier et al. 2017), respectively. In haplotype networks, although haplotypes in a candidate region for the highlands tend to be more frequent in highlanders than in lowlanders (Fig. 3a, supplementary fig. S3, Supplementary Material online), and vice versa (Fig. 3c, supplementary fig. S4, Supplementary Material online), we did not observe a clear set of haplotypes exclusive to highlanders in any Highlands candidate region, nor exclusive to lowlanders in any Lowlands candidate region. This suggests that populations in the Highlands and Lowlands have not reached complete differentiation in these haplotypes, potentially indicating their shared ancestral genetic background or ongoing gene flow. Furthermore, it is likely that the selection observed is acting on standing genetic variation, meaning that the advantageous haplotype was already present in the low-altitude population before becoming selected for in the high-altitude environment.

For some candidate regions, we observe a clear structure between several haplotypes, evidenced by a large number of mutations separating haplotypes that are more frequent in highlanders than in lowlanders, and vice versa (Fig. 3c, supplementary fig. S3, Supplementary Material online, supplementary fig. S4, Supplementary Material online). This could indicate that, for a candidate region identified in one ecoregion, selection is not acting upon, or at least not as strongly, in the other ecoregion at the same genomic location. Furthermore, we find the decay of EHHS to be slower in highlander haplotypes compared to lowlander ones at candidate regions of the Highlands (Fig. 3b, supplementary fig. S5, Supplementary Material online), and a slower decay in lowlander haplotypes compared to highlander ones at candidate regions of the Lowlands (Fig. 3d, supplementary fig. S6, Supplementary Material online). Additionally, at certain candidate regions, haplotypes in both ecoregions decay at similar rates, although highlander haplotypes in candidate regions of the Highlands are present at a higher frequency than in the Lowlands, and vice versa (supplementary fig. S7, Supplementary Material online, supplementary fig. S8, Supplementary Material online). This is consistent with our FCS_between approach, designed to capture regions with both EHHS and high-frequency differentiation between populations of each ecoregion.

Functional Annotation of Candidate Regions for Positive Selection

To functionally characterize the identified candidate regions to be evolving under selection, we conducted a gene-set over-representation analysis on the genes in candidate regions unique to each ecoregion, per FCS computation, using the functional annotation clustering tool from DAVID (Sherman et al. 2022). This allowed us to identify distinct trends of gene clustering whose functions shed light on which biological processes might underlie adaptive events in the Highlands compared to the Lowlands, and vice versa.

In the Highlands, genes obtained from our FCS_within analysis revealed clustering annotation of candidate genes associated with immune response (Allograft rejection), whereas candidate genes from FCS_between against the Lowlands revealed clustering annotation of terms associated with inflammatory response (Acute phase response), constituents of the eye lens (Eye lens Protein), oxidative processes (hydrogen peroxide catabolic process), lipid metabolism (acyl-CoA hydrolase activity), virus transcription elongation (HIV transcription elongation), and cellular signaling and development (Notch signaling pathway) (Fig. 4a and supplementary table S5, Supplementary Material online).

Conversely, in the Lowlands, clustering annotation of candidate genes was primarily observed in terms such as RNA metabolism (RNA metabolic process), obtained with FCS_within, as well as neurotransmitter receptor activity (trace-amine receptor activity) with FCS_between against the Highlands (Fig. 4b and supplementary table S6, Supplementary Material online).

Archaic Introgression in Candidate Genomic Regions Under Selection

Given the known history of archaic introgression in PNG, we looked for overlap between candidate genomic regions under selection and archaic segments identified in 25 whole genomes of PNG highlanders (Malaspinas et al. 2016) using IBDmix (Chen et al. 2020). As the method is only suited for whole-genome sequencing data, we used these genomes to identify regions of archaic introgression in PNG highlanders. Given that we include individuals from the same highland locations reported in Malaspinas et al. and have an overlap of 19 samples shared between data sets, we used the identified archaic regions as a proxy to compare with the selection results of highlanders obtained from the array data. Using the whole-genome data, for each candidate region under selection, we created bins along the genome of the same length and quantified the number of highlander individuals that had segments of archaic introgression in each bin. We then used the extreme of the distribution (top 0.05) of these values to define a threshold for the significance of archaic introgression. Interestingly, two of our candidate regions for selection passed the overlap threshold of archaic introgression. The first significant result is an FCS_within hit (individuals with overlaps = 17; 0.05 threshold = 10 individuals) found within the HLA region in chromosome 6 (Fig. 5a). It comprises several genes and pseudogenes, with functions related to immunity (HCG4, HLA-V, HLA-G, HLA-A, HCG9, HLA-J, TRIM31, TRIM31-AS1, TRIM10, IFITM4P, RNF39) and the less characterized ZNRD1ASP. The second significant hit we identified in highlanders was in a candidate region shared with lowlanders (individuals with overlaps = 11; 0.05 threshold = 10 individuals) in chromosome 10 (Fig. 5b). This region includes mostly protein-coding genes related to a wide variety of functions including protein transport (AP3M1), adenosine metabolism (ADK), chromatin remodeling (KAT6B), phosphatase activity (DUPD1, DUSP13), lipid metabolism (SAMD8), and energetic exchange (VDAC2).

Discussion

In this study, we identified different genomic regions as candidates to be evolving under positive selection in populations from the Highlands and Lowlands of PNG through the analysis of genome-wide genotype data from 62 highlanders and 43 lowlanders from the OGVP.

PNG Highlanders Exhibit Signals Related to High Altitude Adaptation

Given the difference in altitude between populations from the Highlands and the Lowlands, we initially hypothesized that selection could have targeted genomic regions in loci relevant for high altitude adaptation in populations from the Highlands. Previous studies on high-altitude adaptation in Tibetans, Andeans and Ethiopians have implicated genes like EPAS1, EGLN1, and DEC2 (Bigham et al. 2010; Huerta-Sánchez et al. 2013; Li et al. 2019) to be involved in the HIF pathway as the primary targets of selection. Although those genes were not among our candidates of selection in PNG highlanders, we found new candidates that can be playing a role under hypoxic conditions and in processes related to the body's response to high altitude.

The HIF pathway and hypoxia have several effects on different aspects of metabolism like glycolysis, anaerobic metabolism, oxidative stress, fatty acid metabolism and the inflammatory response (Eltzschig and Carmeliet 2011; Ge et al. 2015; McGarry et al. 2018; O’Brien et al. 2020; Taylor and Scholz 2022). One of the top candidate genes under selection found in the Highlands is LDHA (FCS_between). Its encoded enzyme, lactate dehydrogenase A, is integral in glycolysis and plays a crucial role in the metabolism of lactate, while also being particularly active in anaerobic metabolism. Moreover, the transcription factor HIF1 is known to enhance the expression of LDHA (Firth et al. 1995; Semenza et al. 1996), facilitating a transition from oxidative to glycolytic metabolism during periods of hypoxia (Taylor and Scholz 2022). LDHA is also a known biomarker for acute mountain sickness (AMS), showing increased levels of lactate dehydrogenase A in nonaffected individuals (Yang et al. 2022). In high-altitude populations from Peru, expression levels of LDHA have been suggested to predict the presence of CMS as well (Xing et al. 2008). A study on a high-altitude Tibetan population found that increased lactate concentration is positively correlated with high-altitude adaptation, while also suggesting a possible decrease in the activity of fatty acid oxidation at high altitude (Ge et al. 2012).

Hypoxic conditions also favor oxidation processes leading to an increase in reactive oxygen species (ROS), which in turn aids in promoting angiogenesis (Kim and Byzova 2014), the process of new blood vessel growth, which has been found to play a key adaptive role in other high altitude populations (Droma et al. 2022; Ferraretti et al. 2023). Within our top candidate regions under selection, we found the dual oxidase genes DUOX2, DUOX1, and DUOXA1 in chromosome 15 as part of the 13th strongest candidate region (supplementary table S2, Supplementary Material online) and the antioxidant enzyme PRDX1 in chromosome 1 as part of the strongest candidate region (FCS_between) (Table 1) forming an enriched cluster related to hydrogen peroxide catabolic process in our functional annotation clustering approach (Fig. 4a). PRDX1, primarily involved in peroxide elimination accumulated during metabolism, has been identified as a candidate gene under selection in an Andean population residing at intermediate altitudes (Eichstaedt et al. 2015). Moreover, DUOX2, which is expressed mainly in the lung and blood vessels and plays a role in the hydrogen peroxide pathway and ROS, has also been identified as a candidate gene under selection in highland Andeans (Jacovas et al. 2018; Borda et al. 2020), underscoring the importance of the interplay between ROS management and angiogenesis in populations adapted to high altitude.

Inflammation is interdependent with the production of ROS and hypoxia (McGarry et al. 2018). Within the top 10 candidate regions under selection in the Highlands we found the serum amyloid A (SAA) genes SAA1, SAA2-SAA4, SAA4, and SAA2 in chromosome 11 (FCS_between) to be also present in a functional cluster related to the acute phase response (Fig. 4a), a series of events that occur as a response to inflammation. SAA4 and SAA1 have been identified as biomarkers to predict AMS in humans (Julian et al. 2014; Lu et al. 2016). Additionally, SAA is also implicated in the cardiovascular system (Shridas and Tannock 2019), which can be affected by hypoxia and is targeted by selection in other high-altitude populations (Mallet et al. 2021). Moreover, we find the gene BRINP2 among our top 10 candidate regions under selection (FCS_between). BRINP2 is an homolog of BRINP3 (Kawano et al. 2004), a gene associated with cardiovascular phenotypes in humans and found to be under selection in high-altitude Andean populations (Crawford et al. 2017). The high homology between their coded proteins suggests they might perform a similar or interchangeable function (Berkowicz et al. 2016), therefore playing a possibly similar role in high altitude adaptation in PNG highlanders. Furthermore, we also found SANBR (aka KIAA1841) among our candidate genes under selection (FCS_between) (supplementary table S2, Supplementary Material online), a newly identified gene to be under selection in high-altitude Tibetan populations and suggested to play a role in the adaptation of cardio-pulmonary functions (Zheng et al. 2023).

Within the candidate regions under selection in the Highlands we identified four genes associated with the Notch signaling pathway (supplementary table S2, Supplementary Material online), three of them identified by FCS_between (NUMB, PSEN1, RBPJ), found also to form a cluster using DAVID (Fig. 4a), and one of them identified through FCS_within, MAML3 (SM table S2). Of these, NUMB and PSEN1 are notably among the top 10 candidate regions under selection (Table 1). The Notch signaling pathway is a key modulator in cell fate determination in the development of vertebrate and invertebrate species, a regulator of blood vessel structure and angiogenesis, and participates in a cross-talk with the HIF pathway in the cellular response to low oxygen conditions (Landor and Lendahl 2017). It has been proposed as a candidate pathway for high altitude adaptation in multiple organisms (O’brien et al. 2022), including human Tibetan and Andean populations (Bigham et al. 2009; Xin et al. 2020). The NUMB gene is involved in cell fate determination and has been shown to enhance Notch signaling by ubiquitination processes (Luo et al. 2020). NUMB has been identified as a key differentially expressed gene involved in high-altitude adaptation in Tibetan chickens when compared against lowland breeds (Zhang et al. 2020). Although PSEN1 has been primarily studied in the context of Alzheimer's disease, the protein encoded by PSEN1 is essential for the cleavage of Notch receptors during cell signaling and has also been suggested to play a part in the regulation of the oxygen sensing pathway in mice (Kaufmann et al. 2013; Newman et al. 2014). RBPJ is an essential transcription factor in the Notch signaling pathway that is also involved in angiogenesis (Wang et al. 2009; Lake et al. 2014; Ramasamy et al. 2014). It was identified as a candidate gene to be under selection in a study of a high-altitude macaque population and also as a potential novel regulator of hypoxia adaptation in Tibetan pigs (Jia et al. 2016; Szpiech et al. 2021). The gene MAML3 is also a transcriptional mediator of the Notch signaling pathway (Kitagawa 2015; Yamasaki et al. 2020). Although not much research has been done on this gene in the context of adaptation, it has been shown to be upregulated, along with RBPJ, by hypoxic conditions in pancreatic adenocarcinoma cells (Onishi et al. 2016). Moreover, MAML3 and ASPHD2, another candidate gene identified by FCS_between that is among our top 10 candidate regions under selection, are also upregulated together in glioblastoma cells under hypoxic conditions (Seidel et al. 2010).

The acyl-CoA thioesterase (ACOT) gene family encodes enzymes, located primarily in the mitochondria and peroxisomes, that play an important role in lipid metabolism by maintaining the balance of free and activated fatty acids in the body through oxidation processes (Brocker et al. 2010). We find the ACOT genes ACOT2, ACOT4, and ACOT6 to be within the genes in our top 10 candidate regions under selection obtained from FCS_between (Table 1), and as part of a functional annotation cluster (Fig. 4a). Although not much has been studied on these particular gene family in the context of human adaptation, hypoxic conditions have been found to attenuate fatty acid oxidation in adapted high-altitude populations through targeting other genes involved in lipid metabolism (O’Brien et al. 2020).

We also found MMACHC, involved in vitamin B₁₂ metabolism, in the candidate genomic region under selection in chromosome 1 showcasing the strongest signals (FCS_between) (Table 1). Under hypoxia, it has been shown to be significantly repressed by the transcription factor HIF1 (Kiessling et al. 2022).

Thus far our results suggest that, in contrast to other locations where human populations have adapted to high altitude, the HIF pathway has not been the primary target of natural selection related to increased high altitude in the Highlands of PNG. Instead, the observed molecular signatures suggest that these populations might be adapting to high altitude through pathways that mediate the response to inflammatory and oxidative stresses inducible by hypoxia as well as through the Notch signaling pathway, a regulator of the response to low oxygen conditions, with a potential link between these and angiogenesis. Furthermore, we note that ongoing or past selective pressures related to pathogens might have also given rise to signals in genes involved in inflammation and oxidative stress, as these too are processes affected by the systemic defense responses to pathogens (Tschopp 2011). Future research that can identify a causal distinction between the two is needed.

Signals in Genes Related to Eye Components in PNG Highlanders

In populations from the Highlands, we identify five candidate genes related to retina function and eye development: SAMD7, CRYBB1, CRYBA4, SLC7A14, and LIM2 (FCS_between) (supplementary table S2, Supplementary Material online). With SAMD7, CRYBB1, CRYBA4, and SLC7A14 being found within our top 10 candidate regions under selection (Table 1). In our functional annotation clustering approach, we find CRYBB1, CRYBA4, and LIM2 comprising an annotation cluster related to Eye lens proteins (Fig. 4a). SAMD7 is a crucial component in the development of rod photoreceptor cells in the retina (Omori et al. 2017). CRYBB1 and CRYBA4 are two adjacent genes that form part of the crystallin gene cluster on chromosome 22, coding for major proteins of the eye that maintain the transparency and refractive index of the lens (Siggs et al. 2017). Previous studies have found that mutations in any of these two genes can lead to congenital cataracts (Mackay et al. 2002; Billingsley et al. 2006). SLC7A14 is a cationic transporter that is involved in the development of the retina, with mutations associated with retinitis pigmentosa (Jin et al. 2014). Other candidate genes encoding solute carrier proteins (SLC) have also been previously identified as candidates under selection in indigenous populations of Malaysia (Deng et al. 2014). LIM2 codes for an integral membrane protein of the lens with a potential role in cell junction communication (Louis et al. 1989; Tenbroek et al. 1992), and mutations on this gene have also been associated with congenital cataracts (Irum et al. 2016). A 2017 study (Lee et al. 2019) found that PNG has a high global prevalence of blindness in individuals over 50, primarily due to cataracts and uncorrected refractive errors (Lee et al. 2019). The prevalence is notably higher in the Highlands, suggesting health disparities. Although trachoma, an eye disease caused by Chlamydia trachomatis, is endemic in PNG (Bourne et al. 2013; Handley et al. 2018), it is not a primary cause of blindness, implying a potential genetic resistance among the PNG population (Ko et al. 2016). Research in a Gambian population identified the glycoprotein-encoding gene TNC as upregulated during trachoma infection (Natividad et al. 2010), which is also among the candidate genes under selection in PNG highlanders (FCS_between) (supplementary table S2, Supplementary Material online). A recent study revealed that the glycolysis pathway plays a role in adapting human lens epithelial cells to hypoxic conditions (Huang et al. 2023) with LDHA, described above, identified among other genes as upregulated under hypoxia. However, glycolysis-related genes were not prominent in our DAVID annotation clustering analysis. Acknowledging the complexity of eye diseases, which involve both social and demographic factors as well as care access, we concur with other authors that the variations in eye health in PNG may reflect broader health disparities across different regions. Moreover, these findings enhance the understanding of eye diseases in PNG from a genetic standpoint.

Signals Related to the Immune Response in PNG Highlanders

We identified genes linked to the immune system, consistent with evidence of pathogens being a primary force driving evolutionary adaptation across human populations (Fumagalli et al. 2011). In PNG, where infectious diseases like TB, malaria, and pneumonia are prevalent (World Health Organization [WHO] 2018), we found the candidate genes ELL, GTF2H1, and ELOB (supplementary table S2, Supplementary Material online) (FCS_between), in our clustering results in DAVID related to HIV Transcription Elongation (Fig. 4a). It is important to note that while these genes are currently annotated as being involved in HIV infection, given the relatively recent appearance of the virus, we believe that these identified selection candidates might be involved with viral infection response in a wider sense and not necessarily related to HIV. These identified candidate genes are associated with RNA elongation and are crucial in viral infection (Schaeffer et al. 1993; Aso et al. 1995; Kong et al. 2005; Bergeron et al. 2010; Jain et al. 2016; Liu et al. 2020). However, there is currently no conclusive evidence that genetic variations in any of them provide adaptive advantages against any infectious disease. Notably, ELOB, which interacts with HIF1A (Greer et al. 2012), could also be selected for high-altitude adaptation.

Allograft rejection, a process related to the immune response, was an annotation cluster formed by our candidate genes obtained from FCS_within in the Highlands that passed our cluster enrichment thresholds in DAVID. The cluster is comprised of the genes HLA-A, HLA-G, and FASLG (Fig. 4a). FASLG, prominently expressed in B cells and involved in apoptosis (Hahne et al. 1996), has been recognized for its role in immunity against Mycobacterium tuberculosis (Van Rensburg et al. 2017; Loxton and Van Rensburg 2021). In PNG, the incidence of TB has been observed to be higher in Western Province (McBryde 2012), included in this study as part of the Lowlands ecoregion. Moreover, HLA-A and HLA-G genes, central in the immune response, are of note: HLA-G potentially serves as a TB disease state biomarker (Saurabh et al. 2016), and HLA-A alleles have been linked to TB in a Mali population (Kone et al. 2019).

Signals in PNG Lowlanders are Mainly Related to the Immune Response

To put into contrast the pattern of signatures of selection found in PNG highlanders, we also scanned for signatures distinguishable in lowlanders. Given the generally warmer and humid conditions in the Lowlands ecoregion, we initially hypothesized that selection could have targeted genomic regions mainly related to the immune response in populations from this area compared to the Highlands.

The candidate genomic region with the strongest signal of selection in the Lowlands when compared against the Highlands, FCS_between, encompassed genes from the human leukocyte antigen (HLA), a crucial component of the human immune system, HLA-DQB1 and HLA-DQA2 (Fig. 2b, Table 3). This suggests that pathogen-driven selection has likely acted particularly in the PNG Lowland populations as well. Moreover, we found the trace-amine-associated receptor (TAAR) genes TAAR6, TAAR9, and TAAR8 obtained with FCS_between (supplementary table S3, Supplementary Material online) to cluster together in terms related to the trace−amine receptor activity (Fig. 4b). TAARs belong to a family of transmembrane receptors and have been mainly studied for their involvement in functions of the central nervous system (Rutigliano and Zucchi 2020). TAARs have also been suggested to play a role in leukocyte function (Christian and Berry 2018). We find that functional annotation clustering of terms related to RNA metabolic processes comprises the genes RNASE4, RNASE6 and ANG obtained with FCS_within. RNASE4 and RNASE6 code for proteins belonging to the ribonuclease A superfamily, with roles in immune modulation and inflammation processes (Lu et al. 2018). RNase 6 shows high antimicrobial activity against Gram-positive and Gram-negative bacteria (Becknell et al. 2015; Pulido et al. 2016). Whereas RNase 4 contributes to innate defenses in the kidney and urinary tract by showing antimicrobial activity (Bender et al. 2021). The gene ANG also codes for a ribonuclease, RNase 5, that is the same superfamily as RNase 4 and RNase 6. ANG is a mediator of new blood vessel formation through angiogenesis that shows both antimicrobial and antiviral capabilities as well as being involved in the inflammatory response (Hooper et al. 2003; Cocchi et al. 2012; Sheng and Xu 2016).

Shared Signals of Selection Between Ecoregions

Candidate regions under selection shared between the Highlands and Lowlands, among the top 10 signals of any ecoregion (Table 3), include genes related to cell signaling, post-translational modifications, and to a lesser extent antioxidant enzymes and lipid metabolism. Shared selection signals in these genes could hint at pressures that populations from both regions are potentially subjected to, or at adaptive benefits in both environmental contexts. We note that such patterns could also be the result of gene flow across ecoregions, suggesting that gene flow between PNG highlanders and lowlanders has been more significant than previously thought. Moreover, it is also possible that we are detecting selection signatures that were already present in an ancestral population before divergence of the Highlands and Lowlands. Further research is warranted to elucidate this.

Candidate Regions Under Selection in PNG Overlap with Archaic Introgressed Segments

Previous studies have identified an important role of archaic introgression in human adaptation (Racimo et al. 2015; Gittelman et al. 2016). In some cases, this adaptation appears to be related to the immune response directly (Abi-Rached et al. 2011; Urnikyte et al. 2023). Here, we find a proportion of candidate genes under selection within archaic introgressed regions in PNG highlanders (35th strongest signal with FCS_within) to be part of a region in the HLA system (HLA-V, HLA-G, HLA-A, HLA-J, HCG4, HCG9, RNF39, TRIM31, TRIM31-AS1, TRIM10), crucial in the immune response, aligning with previous findings on the role of archaic introgression in immunity in PNG populations (Vespasiani et al. 2022; Yermakovich et al. 2024). This observation provides evidence that adaptive signals in PNG highlanders could be partially attributed to archaic genetic introgression, possibly as a response to local pathogenic challenges or other environmental stressors. Alternatively, the observed patterns could be the result of shared ancestral variation maintained by balancing selection, which is hard to distinguish from introgression in such a complex region like the HLA. We find a second candidate region under selection that is also shared with PNG lowlanders and located within an archaic introgressed segment identified in PNG highlanders. This shared region constitutes the strongest ranking signal in PNG highlanders (FCS_within) and the 7th strongest signal in PNG lowlanders (FCS_within) (Table 3), implying again that some of the adaptive traits in both highland and lowland populations may have deep evolutionary roots originating from archaic hominins. This candidate region comprises the genes VCL, AP3M1, ADK, KAT6B, DUPD1, DUSP13, SAMD8, VDAC2, COMTD1, and ZNF503−AS1, which appear to not be directly centered on the immune response and share no common function among all of them. The candidate region identified in PNG lowlanders is entirely contained within that of the one identified in PNG highlanders (Fig. 5b). It is interesting to note that a majority of archaic segments stack within the ADK gene, in a section of the shared signal that seems exclusive to the highlanders, since it does not directly overlap with the limits of the region identified in lowlanders. ADK codes for the enzyme adenosine kinase, an important regulator of homeostatic and metabolic networks that is mainly involved in energy metabolism, neurodevelopment, and angiogenesis. Given the diverse functions of these genes within the shared candidate region, further research is warranted to better understand their potential contributions to adaptative processes. We also note that the results presented here are restricted to a comparison of candidate regions under selection (identified as particular to highlanders, or shared between them and lowlanders) with regions of introgression identified only in highlander populations, due to limited data availability where whole genomes are required (see “Archaic introgression in candidate regions under selection” in Materials and Methods). Further analyses on whole-genome data from lowlander populations may reveal new patterns of introgression.

Conclusions

Our work focused on studying genetic adaptation in populations from the highlands and lowlands of PNG. The signals observed in the highlands clustered in annotations of diverse functional terms related to the inflammatory response, oxidative processes, constituents of the eye lens, lipid metabolism, virus transcription elongation, and Notch signaling. We propose that even if highlanders of PNG are exposed to moderate altitude, adaptation to lower oxygen conditions might be occurring in genes and processes that respond to the physiological stress caused by increased altitude or by playing a role in modulating the HIF pathway response. Additionally, we identified signatures at genes known to serve as biomarkers for mountain sickness, hinting at the possibility of complex adaptive responses in high-altitude settings. It is interesting to note that some of the candidate genes detected here have also been identified as under selection in other human and animal populations residing at intermediate and high altitudes (Ge et al. 2012; Eichstaedt et al. 2015; Jia et al. 2016; Jacovas et al. 2018; Borda et al. 2020; Zhang et al. 2020; Szpiech et al. 2021; Zheng et al. 2023) strongly indicating convergent evolution in response to living at increased altitude and highlighting the significance of these genes across species and geographic locations. Furthermore, we also identify one candidate region in the Highlands composed by the genes VAMP1, GAPDH, LPAR5, CHD4, IFFO1, and NCAPD2 that has been also identified by André et al in a highland PNG population, with an association to blood composition in the UK biobank André et al. (2024). We report this region as ranking 20th in our analyses and believe that discrepancies in our results stem mainly from differences in the surveyed populations, methodology, and type of genetic data. As for the Lowlands, we found signals mainly related to the immune response, likely reflecting a higher exposure to pathogens in this ecoregion.

While the observed selection signatures in PNG predominantly suggest distinct selective pressures across the two ecoregions, we also identified shared selection signals between the two. These shared signals, rather than pointing to a single set of functional implications, span a diverse array of biological processes. This complexity highlights the need for further investigation to fully understand the evolutionary dynamics and potential interplay of selective forces acting upon these regions. Finally, we identify selection signals within archaic introgressed segments to be mainly related to immune-associated genes. This finding indicates that archaic genetic contributions may have been pivotal in shaping the immune response of modern humans inhabiting islands in the Pacific. Further integration of past and recent demographic events experienced by populations of PNG into appropriate demography and selection models is required to enhance the understanding of adaptive events and their temporal dynamics. In addition, exploring analyses of polygenic selection could offer valuable insights into the genetics of adaptation of PNG populations. By characterizing signatures of positive selection in populations from PNG across different ecoregions, this work contributes to our understanding on the genetic basis of human adaptation in Pacific islanders.

Materials and Methods

Samples and Genotype Data

The genotype data analyzed here derives from the Oceanian Genome Variation Project (OGVP) in which 981 individuals from a broad range of diverse populations across Oceania were genotyped for 1.8 million SNPs with the Illumina Infinium Multi Ethnic Global Array (MEGA) and mapped to the GRCh37 (hg19) version of the human reference genome (Quinto-Cortés et al. 2024). This data was merged with reference continental populations from HGDP (Bergström et al. 2020): Yoruba from Nigeria in Africa, French from France in Europe, and Han from China in East Asia at overlapping sites with plink 1.9 (Chang et al. 2015), solving strand incompatibilities by flipping SNPs on the reverse strand to the forward using snpflip (binary retrieved January 4, 2022 from https://github.com/biocore-ntnu/snpflip).

Quality control on the merged data set was performed using Plink 1.9. We restricted analyses to keep individuals from continental references and the populations reported in this paper, as well as autosomal biallelic SNPs, removing all variants with a missing genotype rate higher than 5% and individual samples with missing genotype rates higher than 10%. We also filtered out SNPs with minor allele frequencies less than 1%. The final merged data set retained 780 individuals and 543,287 SNPs.

Population Structure and Subsetting of PNG Samples

PCA was performed on the OGVP-HGDP merged data set using Plink 1.9 (supplementary fig. S1, Supplementary Material online) and ancestry proportions were estimated using ADMIXTURE 1.3.0 (Alexander et al. 2009). In accordance with the scope of this study, we restricted our analyses to PNG. To reduce the potential confounding impact of admixture on the detection of signals of positive selection, which can introduce novel genetic variations and linkage disequilibrium patterns that may either distort or hide the hallmarks of a selective sweep (Souilmi et al. 2022), we limited our samples to individuals from the main island of PNG (Fig. 1a to c). This decision was informed by the substantial heterogeneity in admixture profiles observed among individuals from the offshore islands (Figure S3 in Quinto-Cortés et al. 2024). Therefore, individuals from the main island of PNG were subsetted, retaining a total of 105 individuals from 12 provinces across the main island and 84 individuals from continental references.

We determined elevation data for populations in the main island of PNG using the geographic coordinates of sampling locations, obtaining elevation data in meters relative to the local mean sea level from the Google Maps Elevation API (Google Maps Elevation API. Retrieved June 7, 2023, from https://developers.google.com/maps/documentation/elevation/start). With this information, we define the ecoregions “Highlands” and “Lowlands” as regions above and below 1000 m relative to the local mean sea level (mLMSL), respectively. Since two provinces with two different sampling locations displayed different altitudes, we further divided them indicating the location (south and central) of each one, leading to a total of 14 different locations (supplementary table S1, Supplementary Material online). With this, we categorized samples from the provinces of Enga, Southern Highlands, Eastern Highlands, Chimbu (Simbu), Hela, Western Highlands, Madang and Eastern Highlands into the Highlands ecoregion, and categorized samples from the provinces of Gulf, Oro, Madang, Central, Western and East Sepik into the Lowlands. We detected two sampling locations in Madang with different altitudes, so we assigned each one to its corresponding ecoregion. In total, the Highlands comprised 62 samples and the Lowlands 43 samples. We note that the terms “PNG highlander” and “PNG lowlander” used throughout the text do not intend to reflect specific cultural, historical, social, ethnic, or linguistic groupings. Instead, they are meant to merely denote human groups living at distinct altitudes.

Analyses of Genomic Signatures of Positive Selection

We investigated signatures of positive selection through population differentiation and haplotype-based methods, characterizing signals particular to populations in each ecoregion. In total, 62 individuals from the Highlands and 43 individuals from the Lowlands were analyzed for signatures of positive selection.

After phasing the merged data set with SHAPEIT4 (Delaneau et al. 2019), haplotype-based statistics were computed using selscan (Szpiech and Hernandez 2014). To look for signals within populations, we computed the iHS (Voight et al. 2006), a statistic based on haplotype homozygosity used to identify areas of the genome that have undergone recent positive selection, and the nSL (Ferrer-Admetlla et al. 2014), which measures the decay of the number of segregating sites with distance from a core allele, in a similar fashion to the incorporation of the extended haplotype homozygosity (EHH) in iHS. Ancestral and derived allele information necessary for these statistics was obtained from 1KGP (The 1000 Genomes Project Consortium et al. 2015). To look for signals between populations we computed the XP-nSL, an implementation of nSL into a cross-population statistic that uses a reference population to compare against (Szpiech et al. 2021), and the PBS, based on differences in allele frequency between two closely related population groups in comparison to an outgroup (Yi et al. 2010). Genome-wide PBS values were calculated using an in-house R script, implementing F_st values (Weir and Cockerham 1984) obtained with PLINK 1.9 from the pairwise comparisons between Highlands, Lowlands and French, using the latter as an outgroup.

To improve the robustness of signals obtained from selection statistics, we integrated the results of compatible statistics into two FCS. With the FCS score being equal to the sum, over two statistics, of –log10 (rank of the SNP in a given statistic/total number of SNPs) following (Lopez et al. 2019). We integrated the ranked results of iHS and nSL as FCS_within, comprising signals obtained from within population statistics (supplementary tables S8 and S10, Supplementary Material online). And integrated the ranked results of XP-nSL and PBS as FCS_between, comprising signals obtained from between population statistics (supplementary tables S9 and S11, Supplementary Material online). Additionally, the integration of these two statistics into a FCS leverages the identification of regions with both EHH and high frequency differentiation between populations of each ecoregion. In the Highlands, results of FCS_between represent signals obtained from comparing against the Lowlands and conversely, results of FCS_between in the Lowlands represent signals obtained from comparing against the Highlands.

To reduce false positives, strongest candidate regions (i.e. regions showcasing the highest FCS values) were identified through an empirical outlier approach as in Caro-Consuegra et al. 2022. For each ecoregion and FCS computation, candidate genomic regions of positive selection were defined as those harboring SNPs above the top 5% of the empirical distribution, only if they also contained at least three other SNPs with scores above the top 1% in a surrounding ±100 Kb vicinity. We kept the initial top 50 strongest candidate regions in each ecoregion, per FCS, for downstream analysis.

We considered candidate regions of positive selection specific to each ecoregion as those not overlapping in 50 Kb or more a candidate region in the top 50 of the other ecoregion, or if one candidate region is entirely contained within another from the other ecoregion. If a shared candidate region was found within the top 50, we assigned that region as “shared” and excluded it from the analysis of both ecoregions, keeping the remaining strong candidate regions that are unique in each ecoregion. Shared candidate regions were analyzed as its own group (supplementary tables S4 and S12 to S15, Supplementary Material online).

Furthermore, SNPs within the strongest candidate regions specific to each ecoregion and those found to be shared between them were annotated with ANNOVAR (Wang et al. 2010). For the top 10 strongest candidate regions unique to an ecoregion, we determined putative candidate genes driving the signal by considering the location of the highest-scoring SNP in each of them. Intergenic signals in the strongest candidate regions were assigned to surrounding genes only if they were at a distance below 50 Kb. Intergenic SNPs that were not assigned to any gene using this criterion were annotated with GeneHancer, a data base of human regulatory elements (enhancers and promoters) and their inferred target genes widely used to annotate noncoding variants (Fishilevich et al. 2017) (supplementary tables S16 and S17, Supplementary Material online). The gene with the highest association score was assigned to each intergenic SNP and the putative type of regulatory element (Promoter, Enhancer, Promoter/Enhancer) was annotated. The function of candidate genes was characterized using UniProt (The UniProt Consortium 2017) and GeneCards (Stelzer et al. 2016). Derived allele frequencies (DAF) were obtained for the top-scoring SNP in the top 10 candidate regions of each FCS approach.

Haplotype Networks and Decay of Site-specific EHH

Haplotype networks were constructed including highlanders and lowlanders for each of the top 10 candidate regions, in the Highlands and lowlands, obtained by FCS_between (supplementary figs. S3 and S4, Supplementary Material online). To limit the number of haplotypes to display, we restricted the networks to include only haplotypes with a minimum frequency of two chromosomes. We used the R package pegas (Paradis 2010) to build networks that connect haplotypes based on pairwise differences.

The site-specific EHHS measures the probability that any two randomly chosen chromosomes from a population are homozygous over a given surrounding chromosomal region of a focal marker (Sabeti et al. 2007; Tang et al. 2007). EHHS was calculated to compare the decay of homozygosity between Highlands and Lowlands around candidate regions obtained by FCS_between, using the R package rehh (Gautier et al. 2017). Both normalized (Tang et al. 2007) and unnormalized (Sabeti et al. 2007) EHHS comparisons between ecoregions were calculated for each of the top 10 candidate regions identified in the Highlands (supplementary figs. S5 and S7, Supplementary Material online), and for each of the top 10 candidate regions identified in the Lowlands (supplementary figs. S6 and S8, Supplementary Material online). Each EHHS calculation began at a focal marker set to the position of the top-scoring SNP in that candidate region. Normalized EHHS facilitates the direct comparison of the rate of decay of homozygosity between population haplotypes by standardizing the initial EHHS value at the focal marker to 1. This removes the influence of different baseline levels of homozygosity, which may vary due to differences in allele frequencies between populations. In contrast, the unnormalized version shows the decay in homozygosity while retaining differences in initial haplotype frequency at and around the focal marker. These differences can reflect genomic regions with EHH in one population that are also differentiated when compared to another, highlighting genomic regions that could have potentially been targeted by selection.

Functional Annotation Clustering

Analyses of functional annotation clustering were conducted separately on the candidate genes located within the remaining unique candidate regions under selection in the Highlands (41 FCS_within, 42 FCS_between), and the remaining unique candidate regions under selection in the Lowlands (40 FCS_within, 48 FCS_between) after filtering out shared signals between the two. For each FCS and ecoregion, these analyses utilized the functional annotation clustering tool provided by the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (Sherman et al. 2022), using the whole-genome as background and high classification stringency. We included OMIM (Amberger et al. 2015), GO (Harris et al. 2004), KEGG (Kanehisa et al. 2017), WikiPathways (Agrawal et al. 2023) and REACTOME (Fabregat et al. 2017) as query databases and only considered clusters with enrichment scores equal to and above 1.3, equivalent to a P-value of 0.05 (supplementary tables S5 and S6, Supplementary Material online).

Archaic Introgression in Candidate Regions Under Selection

Archaic introgression was extracted from whole-genome sequences of 25 Papuan highlander individuals (Malaspinas et al. 2016). For this purpose, we used IBDmix (Chen et al. 2020), a probabilistic method based on identity-by-descent (IBD) designed to identify segments of archaic introgression without need for an unadmixed modern human reference. The data was analyzed per chromosome, and only autosomes were used. The Papuan data was restricted only to bi-allelic single nucleotide variants (SNVs). For the Neanderthal and Denisovan individuals (Meyer et al. 2012; Prüfer et al. 2014), the data was restricted to SNVs, and multi allelic sites were removed. For each archaic hominin, IBDmix genotype data was created using the generate_gt command with the Papuan samples as the modern data. A mask was created for analysis with each of the Neanderthal and Denisova by combining the corresponding minimal filter mask (Prüfer et al. 2014) for each hominin with the 1000 Genomes Project accessibility mask using bedtools (Quinlan and Hall 2010). IBDmix was then run on the data using the ibdmix command, with a minor allele count threshold (-m) of 2 and applying the combined mask for each archaic hominin (-r), and all other values set as default. We created summary files of our results using the summary.sh script provided by IBDmix, filtering with a minimum length of 50 kb and a minimum LOD score of five (supplementary table S7, Supplementary Material online).

Papuans are believed to have introgression from both the Neanderthal and the Denisova, and several regions were called as both Neanderthal and Denisovan. For this reason, our calls were combined into a single “archaic” category. Using the GenomicRanges R package (Lawrence et al. 2013), we looked for overlap between candidate regions under selection and archaic segments. For each candidate region under selection, the number of Papuan samples with archaic segments that overlap in at least 50 kb with the region was quantified. For each region with at least three hits, the genome was divided into bins of the same length as the region, and the distribution of the number of individuals with archaic segments overlapping at least 50 kb with each bin was calculated. The top 5% of this distribution was used as the threshold for significance of a region (supplementary fig. S9, Supplementary Material online).

Supplementary Material

Supplementary material online is available at Genome Biology and Evolution online.

Acknowledgments

We thank all sample donors who contributed to this study; Michael Alpers, former head of the Papua New Guinea Institute of Medical Research (IMR), and George Koki for valuable support from the IMR; Sarah Kaewert, Victor Moreno-Mayar, Diego Ortega del Vecchyo, Alexander de Luna Fors, and Rafael Montiel Duarte for assistance on methods discussion; Santiago Medina Muñoz, Daniela Orozco, and Karla Benítez for helpful feedback. Data generation for this work was supported by the Oxford-Stanford Big Data Initiative and we thank the Mexican Council for Science, Humanities, and Technology (CONAHCyT) for fellowship support awarded to R.G-B., S.V-S, M.C-H, and C.B-J (CVU numbers: 918462, 1242720, 1234769, and 873785, respectively).

Data Availability

Data availability is described in the source publication of the Oceanian Genome Variation Project (OGVP) (https://arxiv.org/abs/2405.09216). Genotype data for the subset of samples analyzed here are available upon request by contacting the corresponding author.

Ethics Statement

Samples were collected under supervision of the Papua New Guinea Institute of Medical Research (IMR) and have been stored with approval from the IMR for future research. Use of the samples for the Oceanian Genome Variation Project (OGVP) was approved by the Oxford Tropical Research Ethics Committee at the University of Oxford (approval reference: OXTREC 537-14). All samples were anonymised and coded at the time of collection according to the region of origin of the participant. Further sampling details and community engagement procedures are described in (Bergström et al. 2017).

Literature Cited