The study was conducted in communities on the eastern slope of Mount Cameroon with varying malaria transmission profiles and topographical characteristics. The area has a forested equatorial climate, which is influenced by the ocean and mountains, with two distinct seasons: a brief dry season (from November to March) and a lengthy rainy season (from March to November) (Wanji et al, 2003). It is characterized by relatively steady temperatures that range from 18°C in August to 35°C in March, and the relative humidity (75%–80%), average annual rainfall (2,625 mm), and precipitation (2,000–10,000 mm) are relatively high (Wanji et al., 2003). Malaria transmission is widespread and perennial in the region, with parasitemia levels being highest during the rainy season and at lower altitudes (Achidi et al., 2008). Plasmodium falciparum is responsible for most infections, whereas Anopheles gambiae is the area’s dominant, most aggressive, and most active malarial vector. Although the Bakweri tribe, of the Bantu ethnic group, is indigenous to this area (Ngwa, 2000), the region’s fertile volcanic soils and large plantations have attracted individuals from other parts of the nation, mostly from the northwest’s Semi-Bantu ethnic group. There is a substantial level of human migration between localities, mainly for educational, recreational, and commercial purposes.

Individuals with asymptomatic or uncomplicated malaria were prescreened using light microscopy and enrolled through cross-sectional surveys in 10 randomly selected rural, semirural, and urban communities between May 2013 and March 2014, as described previously by Apinjoh et al. (2015; 2017). Communities below 200 m, between 385 m and 575 m, and greater than 636 m above sea level were considered to be at low, intermediate, and high altitudes, respectively. DNA extracted (Qiagen, UK) from CF11 leukocyte-depleted (wwarn.org) venous blood was then sequenced on the Illumina HiSeq platform (Illumina, San Diego, CA, USA), and, subsequently, genotyped using well-established methods, as described previously by Amambua-Ngwa et al. (2019).

A total of 239 P. falciparum samples were sequenced on the Illumina Genome HiSeq platform. Library preparation and sequencing were undertaken by the Wellcome Trust Sanger Institute core teams, as described by Auburn et al. (2012), and as part of the MalariaGEN P. falciparum Community Project (https://www.malariagen.net/projects/p-falciparum-community-project). Short-sequence reads with 76 bp were generated, and variant discovery and genotype calling were carried out as described in the P. falciparum Community Project (https://www.malariagen.net/projects/p-falciparum-community-project). Briefly, after removing sequences that mapped to the human reference genome, the sequence reads were aligned to Pf3D7 reference version 3 (ftp://ftp.sanger.ac.uk/pub/project/pathogens/gff3/2016–07/Pfalciparum.genome.fasta.gz) using Burrows–Wheeler Alignment (BWA) mem tools version 0.7.15 (https://github.com/lh3/bwa.git) with -M parameter to mark shorter split hits as secondary.

The Picard version 2.6.0 (http://picard.sourceforge.net) tool was used to clean and mark duplicates of the aligned sequences, and default parameters of GATK base quality score recalibration, version 4.2.6.0, were applied to improve the quality of aligned bases using the pass variants from the Pf Genetic Crosses project (http://www.malariagen.net/data/pf-crosses-1.0). The stats utility from samtools version 1.2 (Li et al., 2009) was used to generate standard alignment metrics for each sample. GATK’s CallableLoci (Depristo et al., 2011), with default parameters, was used to determine the genomic positions callable in each sample.

Variant calling and genotyping were carried out following the GATK best practice pipeline (Auwera et al., 2013) using utilities HaplotypeCaller and GenotypeGCVFs of GATK version 4.2.6.0, respectively. Each variant was assigned a quality score with GATK’s VariantRecalibrator and ApplyRecalibration using the pass variants from the Pf Genetic Crosses project. This step assigned a quality score named VQSLOD to each variant. High values of VQSLOD indicated a higher quality. Variants were annotated using snpEff version 4.1 (Cingolani et al., 2012), with gene annotations downloaded from GeneDB (Logan-Klumpler et al., 2012) at the following URL: ftp://ftp.sanger.ac.uk/pub/project/pathogens/gff3/2016–10/Pfalciparum.noseq.gff3.gz . Each sample was genotyped for polymorphic-coding SNPs across the genome, ensuring that there was a minimum of 5 × paired-end coverage across each variant per sample. Polymorphic sites within hypervariable, telomeric, and repetitive sequence regions were excluded. Vcftools version 0.1.16 (Danecek et al., 2011) and bcftools version 1.10.2 (https://github.com/vcftools/vcftools.git) were used to filter the genomic data. Biallelic high-quality SNPs (VQSLOD ≥ 3) in the coding region loci with a minor allele frequency of at least 2%, individual samples, and SNP sites missing less than 10% across the isolates were extracted and used for downstream analysis as described (Amambua-Ngwa et al., 2019; Abera et al., 2021).

Plasmodium falciparum samples with low-quality sequences and no altitude information were then discarded (19 samples), and 220 high-quality sequences from high (51), intermediate (139), and low altitudes (30) were retained for assessing the effect of altitude on the genetic diversity of parasites, and its contribution to drug resistance. The minimum of 30 samples from each population category in the study is sufficiently large to determine population genomic variation. Analyses of the optimum sample size for population genetic studies (Trask et al., 2011; Hale et al., 2012; Nazareno et al., 2017; Flesch et al., 2018) suggest that even eight samples could be sufficient to accurately report intrapopulation (within-population) and interpopulation (between-population) genetic diversity estimates.

The genotype mixedness of infection (within host/sample diversity) was estimated using the inbreeding coefficient within a population (F_WS) metric as described by Manske et al. (2012):

F_WS = 1 − HW/HS, (1)

where HW represents the within-isolate expected heterozygosity, calculated from the relative allele frequencies for all genic SNPs, and HS represents the heterozygosity of the local population. The F_WS values range from zero (high diversity of infection) to one (single infection within the sample). For this analysis, individual alleles with coverage < 5 reads and positions with total coverage < 20 reads were classified as missing data. Isolates with > 10% missing SNP data and SNPs with > 10% missing isolate data were discarded. Isolates with F_WS scores > 0.95 were classified as having a single predominant genotype due to a low genome-wide diversity.

Principal component analysis (PCA) was used to estimate the population structure. This method simplified the complexity of high-dimensional data by geometrically projecting them onto a limited number of lower dimensions called principal components (PCs). The first PC minimizes the total distance between the data, as it maximizes the variance between the projected points; PC1 thus has the greatest variance as compared with other PCs, followed by PC2, and so on. However, all the PCs are not typically useful because much of the variance in the patterns of the data was covered by the first few PCs (Lever et al., 2017). In addition, the glPCA function in the R version 3.6.3 software package (The R Foundation for Statistical Computing, Vienna, Austria) was used to calculate the first 10 principal component axes, and the first three PCA axes were retained to explain the genomic variation among P. falciparum populations. The data was thinned down by pruning SNPs with pairwise linkage disequilibrium (LD) R ² values greater than 0.05. An allele-sharing matrix was then calculated using the pruned SNP loci with the R glPCA function that uses between- and within-group variances to determine the population’s genetic structure. Discriminant analysis of principal components (DAPC), which determines ancestry proportions and membership probability, was carried out using the adegenet version 2.1.10 package in R, as described previously (Jombart et al., 2010). This method identified and described population genetic clusters and their relationships and modeled genetic variation/allele sharing across space (altitude) to determine admixture. The optimum PCs for DAPC analysis, which prevent overfitting and instability in the membership probabilities or exclude useful information from the analysis, were selected using the optim.a. score function (Jombart et al., 2010).

Analyses of allele frequency distributions between population F_ST values (Weir & Cockerham, 1984) were calculated using vcftools version 0.1.16 after excluding SNPs missing more than 10% data, as described by Abera et al. (2021), correcting for differences due to population size variation. The F_ST values range from zero to one, and are indicative of there being no and complete population differentiation, respectively.

Allele frequency distributions within the population were assessed by estimating Tajima’s D indices (Tajima, 1989) using vcftools, after removing missing data per SNP. Tajima’s D values were determined for each SNP, and the average value was calculated for each gene using the tapply function in R. Genes with at least five SNPs and Tajima’s D values greater than one were considered to be under positive balancing selection. Standardized integrated haplotype score (iHS) analysis was undertaken to identify positive directional selection signatures using SNP data with allele frequency > 5%. iHS was estimated using the rehh version 3.2.2 package in R with default parameters (Gautier & Vitalis, 2012), after imputing missing SNP data using beagle version 5.2 (Browning & Browning, 2007). An |iHS| > 2.5 (top 1% of the expected distribution) (Voight et al., 2006) was used as a cutoff value to report genes under recent directional selection as reported for genome analysis of West African P. falciparum (Oyebola et al., 2017).

To fully capture evidence of potential selection signature differences in P. falciparum populations at high, intermediate, and low altitudes, the cross-population Rsb statistics of rehh version 3.2.2 package in R were used, as iHS may not be suitable for detecting positive selection for those SNPs that have reached or are near to fixation in the local P. falciparum population (Voight et al., 2006). An Rsb score of < −5 or > 5 indicates that there is strong evidence of potential selection signature differences among P. falciparum populations (Duffy et al., 2015). The Rsb metrics (Tang et al., 2007) in the rehh package (Gautier & Vitalis, 2012) were used to detect the directional selection signature difference between populations. The scan for population-specific directional selection was performed using SNPs with a minimum minor allele frequency of 0.05, and to assess the relative haplotype lengths between two populations standardized against the genome-wide average, as used previously (Mobegi et al., 2014; Duffy et al., 2015). Calculation of haplotype breakdown during Rsb scores was terminated if the gaps between adjacent SNPs were > 20 kb. For determining Rsb, putative selection windows were defined by calculating the distance over which the extended haplotypes decayed to a level of 0.05 in both directions. For the analysis of between-population differences using the Rsb metric, |Rsb| > 5 was used to identify high-scoring SNPs.

Of the 239 P. falciparum isolates sequenced, 19 were removed due to low-quality sequences or a lack of altitude information. High-quality sequence data obtained from 220 samples enabled the identification of 682,462 (51 samples), 681,773 (139 samples), and 681,909 (30 samples) biallelic SNPs at high, intermediate, and low altitudes, respectively, with less than 10% missing SNPs and missing individual isolate data. Sequences from the intergenic regions had a lower read coverage than the coding regions, as expected, because of extreme A + T nucleotide contents in the intergenic regions. Thus, 78.23% (534,085), 78.4% (533,716), and 78.3% (533,953) of all SNPs identified from high, intermediate, and low altitudes, respectively, were located within the genes. Approximately 3,608, 3,644, and 3,847 genes analyzed at high, intermediate, and low altitudes, respectively, had at least one SNP ( Figure 1 ).

Approximately 24,214, 24,426, and 28,613 SNPs at high, intermediate, and low altitudes, respectively, were polymorphic in at least one sample of P. falciparum. Of these, 43.8%, 44.4%, and 45.4% were non-synonymous-coding SNPs; 20.9%, 20.7%, and 21.3% were synonymous-coding SNPs; 28.31%, 28.26%, and 27% were in intergenic regions; 3.3%, 3.04%, and 2.89% were in other intragenic regions; and 3.66%, 3.6%, and 3.46% were in intron regions at high, intermediate, and low altitudes, respectively ( Figure 2 ; Supplementary Table S1 ).

Single-nucleotide polymorphisms with a minor allele frequency (MAF) < 5% were common in all analyzed P. falciparum populations following the exclusion of monomorphic SNPs. Furthermore, SNPs with minor alleles frequency < 5% occurred more frequently in samples from low altitudes than in isolates from high and intermediate altitudes ( Figure 3 ).

The distribution of F_WS scores was similar across isolates in each population, ranging from 0.363 to 0.997 (mean = 0.829, median = 0.937), 0.275 to 0.998 (mean = 0.831, median = 0.941), and 0.333 to 0.997 (mean = 0.772, median = 0.820) for P. falciparum infections at high, intermediate, and low attitudes, respectively ( Figure 4 ; Supplementary Table S2 ). The proportion of isolates with F_WS > 0.95, indicative of single genotype-dominated infections, was 25 (49%), 67 (48.2%), and 9 (30%) at high, intermediate, and low altitudes, respectively. No significant difference was observed in the mean F_WS scores between P. falciparum populations from high and intermediate altitudes (p = 0.61), high and low altitudes (p = 0.93), and intermediate and low attitudes (p = 0.76). The majority of the samples had multiple infections or high levels of genetic diversity in all populations, which was consistent with high malaria prevalence and transmission rates reported throughout the years in the study area.

The principal component analysis did not clearly distinguish the parasite isolates into their respective geographical origins/altitudinal groups from where they were sampled ( Figure 5 ). Admixture analyses were consistent with the outcome of PCA, with no separation between isolates from the three altitudes, as there was gene flow among the populations.

This admixture analysis showed that none of the three components were differentiated according to the optimized cluster value of K = 3. Multiple parasite subpopulations were detected in parasites from all altitudes of Mount Cameroon ( Figure 6 ). There was bidirectional gene flow among the isolates from high, intermediate, and low altitudes of Mount Cameroon, which is consistent with other studies carried out in West Africa and across the continent (Duffy et al., 2015; Amambua-Ngwa et al., 2019). These studies reported there were high probabilities of gene flow within or between the P. falciparum populations of Middle and West African countries.

The clustering of all parasite isolates was consistent with the fixation index (F_ST) values with or without correcting for sample size. The F_ST values of isolates from high versus intermediate altitude were 0.0054 (p < 2.2 × 10^–16), 0.0087 (p < 2.2 × 10^–16) for high versus low altitude, and 0.0075 (p < 2.2 × 10^–16) for intermediate versus low altitude ( Table 1 ) This suggests a low population differentiation among P. falciparum isolates from different altitudes along the slopes of Mount Cameroon.

The average Tajima’s D values for high, intermediate, and low altitude isolates were 0.01 (p = 0.16), 0.12 (p < 2.2 × 10^–16), and −0.02 (p = 0.004), respectively, across the entire genome. Approximately 117, 134, and 98 genes in the P. falciparum population at high (Supplementary Table S3 ), intermediate (Supplementary Table S4 ), and low (Supplementary Table S5 ) altitudes, respectively, had at least five SNPs with Tajima’s D values > 1, suggesting that they were most likely under balancing selection. It is notable that 93 of the 98 genes with a Tajima’s D value > 1 at a low altitude also had a value > 1 at intermediate and high altitudes (Supplementary Table S6 ). Among these genes, crt, ama1, eba175, msp1, trap, dblmsp, clag2, and clag8 have been previously reported to be under balancing or mixed selections. Of these genes, clag8, ama1, and trap encode surface proteins that undergo balancing selection and have previously been identified as potential vaccine candidates. The Duffy-binding-like merozoite surface protein 2 gene (pfdblmsp2) has been implicated in parasite resistance to halofantrine, mefloquine, and lumefantrine.

Genome-wide analysis for evidence of recent positive directional selection due to drug pressure, the impact of the immune response, or other mechanisms identified 36, 49, and 39 genes at high (Supplementary Figure S1 , Supplementary Table S7 ), intermediate (Supplementary Figure S2 , Supplementary Table S8 ), and low (Supplementary Figure S3 , Supplementary Table S9 ) altitudes, respectively, with at least two SNPs that could be under significant positive selection. Interestingly, 22 consistent sweeps were present in populations at all altitudes in the study area ( Table 2 ). The genes surfin8.2 on chromosome 8 and trap on chromosome 13, previously reported as erythrocytic and pre-erythrocytic vaccine candidates, respectively, were detected with higher SNPs at elevated |iHSs|, indicating continuous selection on these antigenic loci. In addition, other vaccine candidates such as the celtos, clag8, and ama1 genes, were also detected, as was the chloroquine-resistance transporter gene crt, although only at low and intermediate altitudes. No signatures around the multidrug resistance mdr1 gene and the antifolate drug target genes dhfr and dhps were detected.

Of the genes under positive directional selection, 26 of the 36 at high (Supplementary Table S10 ), 30 of the 49 at intermediate (Supplementary Table S11 ), and 24 of the 39 at low altitudes (Supplementary Table S12 ), showed both positive directional and balancing selections.

Overall, 17 of the 22 genes under positive directional selection from the three altitudinal zones in the study area showed both positive balancing and directional selections ( Table 3 ), among which are SURF4.2 on chromosome 4, cytoadherence-linked asexual protein 8 (CLAG8) on chromosome 8, and thrombospondin-related anonymous protein (TRAP) on chromosome 13. There was no detectable signature around the multidrug resistance mdr1 gene, and the antifolate drug target genes dhfr and dhps showed both positive balancing and directional selections.

Chromosomal loci with population-specific evidence of directional selection between high and intermediate altitudes were identified: five at a high altitude and four at an intermediate altitude ( Figure 7 , Supplementary Table S13 ). Selection signatures were observed on chromosome 4 for both high and low altitudes. Similarly, six chromosomal loci with population-specific evidence of directional selection between high and low altitudes were identified: three each at high and low altitudes ( Figure 7 , Supplementary Table S14 ). In addition, 11 chromosomal loci with evidence of directional selection between intermediate and low altitudes were detected: four at intermediate and seven at low altitudes ( Figure 7 , Supplementary Table S15 )

Population-specific strong recent directional selections were also observed at genomic regions of chromosomes 1 and 4 in a P. falciparum population from high altitude, which includes two genes (PF3D7_0104100 and ACS6). Similarly, specific strong recent directional selections were found at genomic regions of chromosomes 4 (containing one gene, SMC3), 6 (containing one gene, ACS12), 7 (containing one gene, PF3D7_0709700), 12 (containing one gene, CelTOS) and chromosomes 1 (containing one gene, PF3D7_0114500), 6 (containing one gene, PF3D7_0616400), 7 (containing four genes, CG2, PF3D7_0710000, PF3D7_0711500, and PF3D7_0713200), 8 (SURF8.2), 11 (containing one gene, AMA1), 13 (containing one gene, TRAP) of P. falciparum populations from intermediate and low altitudes, respectively (Supplementary Table S16 ). These genes include immune-targeted genes such as AMA1, TRAP, SURF8.2, and CelTOS. This observation was consistent with reports of there being specific strong recent directional selections in immune-targeted genes of West African P. falciparum populations (Duffy et al., 2015).

Non-synonymous SNPs in resistance-associated genes were identified in isolates from high, intermediate, and low altitudes ( Table 4 ). Four CQ-resistant alleles (pfcrt-K76T, pfcrt-A220S, pfcrt-Q271E pfcrt-I356T, and Pfmdr1-N86Y) and SP resistance-associated mutations in pfdhfr and pfdhfs genes were detected in isolates at high, intermediate, and low attitudes. The prevalence of pfdhr-N51I and pfdhr-C59R, the major pyrimethamine resistance-conferring alleles, was remarkably high at all attitudes (above 98%). Moreover, there was a variation in the prevalence of the drug resistance-conferring alleles Pfdhps-K142, Pfdhps-I431V, pfdhs-G437A, and Pfdhps-A581G, with pfdhps-I431V and pfdhps-A581G being the most prevalent ( Table 4 ).

Pfk13 mutations associated with artemisinin resistance showed a variation in all analyzed P. falciparum populations. The Pfk13-K189T mutation, previously identified in African P. falciparum populations (Verity et al., 2020; Abera et al., 2021), was observed at higher frequencies (> 35%) more than its pfk13-K189N (<7%), pfk13-A578S (<3%), pfk13-A175T (<3%) and pfk13-L119 (<3%) counterparts. In addition, Pfk13-K189T was observed in 47%, 36%, and 40% of isolates at high, intermediate, and low altitudes, respectively.