2.1. Sample information

This study includes 149 newly genotyped individuals from two ethnolinguistic groups in Namibia and Botswana who passed quality control and relatedness filtering: the Khwe (114 individuals) and the Kwangali (35 individuals) (Figure S1; Table S1). The Khwe speak a Khoe-Kwadi language [38] and have traditionally subsisted on foraging in areas of the Okavango River Basin spanning across southeastern Angola, northeastern Namibia, and Botswana’s Okavango Delta [39]. Based on linguistic variation and livelihood patterns, the Khwe can be divided into two main subgroups, both of which are represented in our sample: the ǁAni-Khwe (30 individuals), primarily fishermen currently residing near the Okavango Delta in northwestern Botswana; and the non-ǁAni-Khwe (84 individuals), who speak closely related ǁXom, ǁXoo, and Buga varieties, and historically inhabited the drier regions to the east and north of the Okavango Delta [40] (Figures S1 and S2).

The Kwangali are a West Bantu-speaking group within the Kavango language cluster (K.33), who primarily practice agriculture along the Namibian and Angolan banks of the lower Okavango River [41].

All data were collected under the TwinLab collaborative network linking CIBIO/InBIO with the University of Namibia and the University of Botswana. Official research permits were granted by the Ministries of Health of Namibia (RPIV00722019) and Botswana (HPRD: 6/14/1). At each sampling site, the study was explained to participants in their native language by a bilingual translator. Written informed consent was obtained from all volunteers prior to sample collection. From each participant, we collected a saliva sample along with ethnolinguistic information, including ethnic affiliation, place of birth, and language spoken up to the grandparental generation.

2.2. Genotyping and quality control

We genotyped all sampled individuals using the Affymetrix Axiom Genome-Wide Human Origins Array [42], and performed quality control and filtering with PLINK v1.9 [43].

Non-biallelic single nucleotide polymorphisms (SNPs), SNPs with call error rates exceeding 10%, and non-autosomal SNPs were excluded. After filtering, a total of 553,841 SNPs remained. No individuals had missing rates greater than 10%. To filter out closely related individuals, we excluded one individual from each pair with estimated relatedness (PI_HAT) > 0.2 [43].

For comparative purposes, we merged our data with that of 130 individuals across 14 contextual African populations previously genotyped on the same array, also using PLINK v1.9 (Figure 1; Table S1). The intersection between the two datasets yielded a total of 535,116 SNPs.

To independently genotype SNPs rs2814778 and rs12075, which determine Duffy allelic variation, we used Sanger sequencing following polymerase chain reaction (PCR) amplification with previously described primers [7].

Figure 1. Population structure and admixture in the Khwe. (A) Map showing the approximate geographic locations of newly sampled and contextual populations. (B) PCA plots of the populations shown in the map. In the right panel, only Khwe individuals and parental populations are highlighted in colour. (C) Individual ancestry proportions inferred using RFMix v2. Vertical lines represent the estimated ancestry proportions for each individual.

2.3. Phasing

Phasing was performed with Beagle 5.4 [44], employing the standard HapMap genetic map available on the Beagle website (https://faculty.washington.edu/browning/beagle/beagle.html).

To ensure results were not dependent on a specific phasing method, we also conducted analyses with SHAPEIT v2.r904 [45], using the HapMap b37 genetic map available on GitHub (https://github.com/odelaneau/shapeit4/tree/master/maps).

2.4. Population structure analysis

To examine population structure, we performed principal component analysis (PCA) using a dataset pruned for linkage disequilibrium (LD) with PLINK v1.9 [43]. The pruned dataset comprises 343,959 SNPs and was generated by removing SNPs with an r² value greater than 0.4 within 200-kb windows, shifted at 25 SNP intervals. Principal component scores were computed with PLINK and visualized with a PlotPCAs.r script.

Pairwise F_ST values between populations were calculated with EIGENSOFT [46,47] and visualized with a heatmap and Neighbour-Joining clustering using the R-packages ape and pheatmap [48–50].

2.5. Ancestry inference

Ancestry inference at the population, individual and local-levels was performed using RFMix v2 [51] with phased samples and a three-way admixture model.

The analysis included 114 Khwe individuals and employed four reference panels representing ancestries related to Bantu-speakers, indigenous southern African hunter-gatherers and eastern Africans. To maintain balanced panels, the number of individuals included from each source population was constrained by the smallest available sample size. The first panel consisted of 13 Kwangali (Bantu-related ancestry), 13 Juǀ’hoan (southern African hunter-gatherer-related ancestry), and 13 Somali (Afro-Asiatic eastern African-related ancestry) individuals.

To assess the influence of alternative parental populations on the detected local ancestry deviations, we tested three additional reference panels using the Hadza and Yoruba populations as alternative eastern African and Bantu-related source populations, respectively. The panels consisted of: (i) 13 Nyaneka individuals, 13 Juǀ’hoan individuals, and 13 Somali individuals; (ii) 15 Kwangali individuals, 15 Juǀ’hoan individuals, and 15 Hadza individuals; and (iii) 14 Nyaneka individuals, 14 Juǀ’hoan individuals, and 14 Hadza individuals.

Although the Hadza are hunter-gatherers with genetic profiles that differ substantially from those of most other eastern African populations [52,53], they were included as an alternative reference population in order to assess the robustness of local ancestry deviations to the choice of parental groups, rather than to model the eastern African pastoralist contribution directly.

We used a two SNP window size, three Expectation-Maximization iterations, and the option “reanalyze-reference” to account for possible admixture in the reference groups. To test the robustness of our results, we assumed two admixture times: 39 generations corresponding to the average between proposed dates of 60 generations ago for the initial admixture between East Africans and indigenous southern Africans, and 17 generations ago for the later admixture with Bantu-speaking groups [36]; and 25 generations corresponding to the average of 33 and 17 generations ago for the later admixture event [36,37].

Individual ancestry proportions were calculated by averaging the chromosome-specific ancestry estimates provided in RFMix Q output files and were visualized using a Python script (VisualizeAdmixture.rfmix.py). Population-level ancestry was then determined by averaging the individual ancestry proportions within each group.

Local ancestry assignments were performed using the RFMix fb output files, which provide marginal probabilities for each genomic window originating from a given parental population. Windows with marginal probabilities less than 1 represented 0.7% of total assignments (0.35%, 0.12% and 0.23% for Bantu, eastern African and southern African ancestries, respectively) and were excluded to remove segments with ambiguous ancestry assignments. Results were unchanged when local ancestry inference was performed without filtering low-confidence windows. SNPs located in centromeric regions and within 2 Mb of chromosome telomeres were additionally filtered out [54–56].

2.6. Assessing deviations from expected admixture allele frequencies

To evaluate whether the frequency of FY*B^ES in the Khwe population was higher than expected from admixture alone, we employed a binomial resampling approach, where the probability of success was defined as the product of the FY*B^ES allele frequency in the ancestral populations and their estimated admixture proportions [57]. The resampling procedure is implemented in the resampling.CLI.plots.py script.

To compare the FY*B^ES frequency elevation at the Duffy locus with deviations observed at other SNPs genotyped on the Genome-Wide Human Origins Array, we used the F_adm statistic, which measures, for each SNP, a normalized difference between observed allele frequencies in an admixed population, and frequencies expected from admixture proportions under neutrality [58].

2.7. Simulations of genetic drift

We performed forward simulations under the Wright–Fisher model to generate expected distributions of FY*B^ES allele frequencies in the absence of selection, using the Wright_Fisher_simulator.py script. The simulated admixed population was modelled as initially formed through contributions from resident southern African hunter-gatherers and incoming East African herders, followed by a secondary admixture event involving Bantu-speaking populations. We investigated the outcomes of varying effective population sizes, parental allele frequencies, admixture dynamics, and number of generations since admixture on the resulting allele frequency distributions.

2.8. Estimation of selection coefficient

We applied a deterministic framework, using equation 10.25 from [59] to estimate selection coefficients (s) that could account for the observed FY*B^ES allele frequencies, given different initial allele frequencies and varying number of generations since admixture. Fitness (W) was parameterized relative to FY*B^ES homozygotes as: W_FY∗B^ES_/FY*B^ES=1; W _{FY∗ B}^ES _/FY*Other=1−hs, and W_{FY*Other /FY*Other}=1−s, where FY*Other designates non-FY*B^ES alleles. Under this parametrization, h=1 corresponds to a fully recessive advantageous effect of FY*B^ES, with selection favouring only FY*B^ES homozygotes.

Allele frequency trajectories were calculated across a range of s values from 0 to 0.4, in increments of 0.001, using dominance coefficients (h) set at 1 (recessive model) or 0.5 (additive model). The selection coefficient that yielded the FY*B^ES allele frequency closest to the observed value was taken as the best estimate. The procedure was implemented using the Selection_coefficient_estimator.py script.

3.1. The admixture structure of the Khwe population

We assessed the population structure of the Khwe in comparison to other Khoe-Kwadi–speaking groups, as well as selected populations representing resident southern African, eastern African, and Bantu-related ancestries (Figure 1). In PCA, the Khwe are positioned between Bantu-speaking groups and Kx’a- and Tuu-speaking southern African hunter-gatherers (PC1 in Figure 1B), while also showing slightly greater proximity to eastern Africans than the Bantu populations included in this study (PC2 in Figure 1B). Pairwise F_ST estimates are consistent with the patterns of genetic differentiation observed in PCA (Figure S3).

Previous studies have shown that the genetic composition of the Khwe and other Khoe-Kwadi–speaking populations can be modelled as the result of an initial admixture between East African and indigenous southern African–related ancestries, followed by a second admixture event with incoming Bantu-speaking groups [36]. Assuming a generation time of 30 years [60], inferred dates for these events range from approximately 2,100–1,800 years ago (~70–60 generations ago) for the admixture between southern and eastern African ancestries, and from 1,000–500 years ago (~33–17 generations ago) for the admixture between the pre-Bantu hybrid population and the Bantu migrants [36,37,61,62].

Using RFMixv2 [51] and selecting the Kwangali, Juǀ’hoan, and Somali as proxies for the three parental populations, we estimated that Khwe individuals carry, on average, 58.0% Bantu-related ancestry (standard deviation (SD) 9.8%), 30.5% (SD 8.6%) resident southern African–related ancestry, and 11.5% (SD 2.6%) eastern African–related ancestry (Figure 1C; Table S2).

Our choice of parental populations was based on geographic, historical, and genetic criteria (Figure 1A; Figure S1). The West Bantu–speaking Kwangali, who inhabit the middle reaches of the Okavango River Basin, and the Kx’a-speaking Juǀ’hoan, who live to the west of the Okavango Delta, both show strong genetic affinities with groups that have a documented history of contact and admixture with the Khwe [35,37,63,64]. In addition, the Afro-Asiatic–speaking Somali share approximately 62% of their genome with a ~3,100-year-old individual from the Luxmanda pastoralist site in Tanzania, whose ancestry is closely related to the eastern African component found in modern Khoe-Kwadi–speaking populations [37,65].

As Khwe-speakers fall into two broad ethnolinguistic subgroups with slightly different genetic profiles (Figure 1B; Figure S2), here referred to as ǁAni-Khwe and non-ǁAni-Khwe, we additionally report results for these two sample subsets. On average, ǁAni-Khwe individuals show higher Bantu-related ancestry (62.6% vs. 56.3%) and lower resident southern African–related ancestry (25.5% vs. 32.2%) compared to non-ǁAni-Khwe individuals (Table S2). Despite these differences, the overall ancestry profiles are consistent with previous estimates—obtained with different methods and smaller sample sizes—showing that the Khwe carry a higher proportion of Bantu than resident southern African–related ancestry [36,37] (Table S2).

3.2. Elevated FY*BES frequency in the Khwe

We next investigated whether the frequency of FY*B^ES in the Khwe was higher than expected from admixture alone [57,66].

Because the SNPs rs12075 and rs2814778, which determine Duffy allelic variation, are not included in the Affymetrix Axiom Genome-Wide Human Origins array, we directly genotyped these two polymorphisms in the 114 Khwe and 26 Kwangali (Table S3). For the Somali and Juǀ’hoan populations, we used previously published FY*B^ES allele frequencies [67,68] (Table S4).

Using a binomial resampling procedure with 10⁶ replicates, we found that the observed FY*B^ES frequency of 94% in the Khwe is significantly higher (p < 10⁻⁶) than the 72% expected on the basis of their estimated ancestry proportions (Table S2) and the current FY*B^ES allele frequencies in the Kwangali (100%), Juǀ’hoan (11%), and Somali (92%) source populations (Figure S4; Table S4). Previous estimates of Khwe ancestry proportions based on different methods and parental populations consistently yielded expected FY*B^ES frequencies below 72% (Table S5) [36,37,69].

We further note that the use of present-day FY*B^ES frequencies in both resident southern African and eastern African reference populations is likely to overestimate the allele frequencies in the ancestral populations that contributed to the Khwe. In particular, FY*B^ES was likely introduced into Juǀ’hoan and other resident southern African groups through relatively recent admixture with Bantu-speaking populations [70].

The situation in eastern Africa is more uncertain: high FY*B^ES frequencies have been documented not only in the Somali, but also in other eastern African pastoralist populations, including the Maasai from Kenya (0.93; [71]) and populations across Sudan (average 0.92), where selection signals have been reported in several groups, including Afro-Asiatic speakers [24]. These observations suggest that FY*B^ES was already widespread in parts of eastern Africa prior to admixture between Bantu-speaking populations and resident eastern Africans around ~1,000 years ago [62,65,72]. However, admixture dates between Nilo-Saharan and Afro-Asiatic groups (~1,000–1,500 years ago [24]) still postdate the proposed arrival of eastern African pastoralists in southern Africa ~2,000 years ago, making it likely that ancestral eastern African populations carried FY*B^ES at lower frequencies than those observed today.

Accordingly, the expected FY*B^ES frequency in the Khwe under neutrality likely falls between two extremes: a lower bound in which the pre-Bantu eastern and southern African ancestral populations entirely lacked FY*B^ES, yielding an expected frequency close to the Bantu ancestry proportion (~58%), and an upper bound based on present-day Somali frequencies, yielding the more conservative estimate of ~72%. Under either scenario, the discrepancy with the observed FY*B^ES frequency of 94% remains substantial.

To further evaluate this deviation, we assessed the genome-wide distribution of the per-SNP F_adm statistic [58] to determine whether the divergence at the Duffy locus was unusually extreme compared to other loci. Notably, only 0.03% of all SNPs exhibited F_adm values equal to or greater than that of rs2814778, the SNP defining the FY*B^ES allele—even under the most conservative assumptions for parental allele frequencies: Kwangali 100%, Juǀ’hoan 11%, and Somali 92% (Figure S5).

In a less conservative scenario, assuming that the FY*B^ES allele was absent in the pre-Bantu hybrid population, rs2814778 displayed the highest genome-wide F_adm value (Figure S5).

While the empirical genome-wide distribution of F_adm values may be influenced by SNP ascertainment, the outlier status of the Duffy locus suggests that the elevated frequency of FY*B^ES is unlikely to be explained by neutral evolutionary processes alone.

3.3. Effect of genetic drift

To further assess whether the elevated frequency of FY*B^ES in the Khwe could be explained by genetic drift, we conducted simulations under a Wright–Fisher model. We explored a general framework in which a pre-Bantu population is first formed through a single pulse of admixture between southern African foragers and eastern African herders, followed by subsequent genetic contributions from Bantu-speaking migrants (Figure 2; Table S6).

Bantu ancestry was modelled as a second discrete admixture pulse (Figure 2A) or as continuous gene flow into the pre-Bantu population (Figure 2B), consistent with the spread of Khwe individuals along the PC1 axis in Figure 1B. Each admixture scenario was tested using two different sets of FY*B^ES allele frequencies in the parental populations: a more conservative set based on current frequencies in Juǀ’hoan (11%) and Somali (92%) (Figure 2C-D and G-H), and a less conservative set assuming that the FY*B^ES allele was introduced solely via Bantu-speaking groups (Figure 2E-F and I-J). All simulations used a previously estimated date of 60 generations for admixture between southern and eastern African populations [36], while the timing of Bantu admixture was allowed to vary between 60 and 25 generations ago, spanning the upper range of estimates for admixture between Bantu speakers and pre-Bantu groups in southern Africa [36, 37, 61,62,73].

Each admixture scenario and set of parental allele frequencies was further combined with varying effective population sizes (Figure 2D,F,H,J) to assess the probability of reaching the observed FY*B^ES frequency of 94% by drift alone. We conducted 100,000 replicates for each of 44 demographic configurations (Table S6).

Except for a single highly conservative condition involving early Bantu admixture, small effective population sizes (N=200), and high initial FY*B^ES frequencies in the pre-Bantu hybrid population (Figure 2C; scenario 1 in Table S6), all remaining 43 tested scenarios resulted in fewer than 5% of simulations reaching or exceeding the observed 94% frequency (Figure 2; Table S6).

Using a more conservative 99% quantile threshold, only five additional scenarios were compatible with genetic drift under conditions of small effective population size (scenarios 2, 3, 12, 23, and 24 in Table S6).

These findings suggest that the high FY*B^ES frequency observed in the Khwe is unlikely to be fully explained by genetic drift alone and are consistent with positive selection acting on the allele after admixture.

Figure 2. Simulations of admixture and genetic drift. Admixture models with two admixture pulses (A) and with one pulse followed by continuous gene flow (B). M1=East-African-related ancestry proportion; M2=Bantu-related ancestry proportion; m=migration rate from the Bantu donor population per generation; T1=number of generations between the first admixture pulse and the onset of Bantu admixture; T2=number of generations from Bantu admixture to present. (C-J) Distributions of simulated FY*BES allele frequencies in samples of 114 individuals, under different admixture models, with varying time intervals since admixture (T1 and T2), initial allele frequencies and population sizes (N), corresponding to scenarios outlined in Table S6: panels C-F refer to scenarios 1 to 22, and panels G-J refer to scenarios 23-44. In each panel, the dashed grey line indicates the observed allele frequency in the Khwe; the solid black line indicates the expected allele frequency; and the grey shaded area represents the 95% quantile of simulated allele frequencies.

3.4. Ancestry deviations at the Duffy locus

After admixture, the ancestry of a source population carrying a beneficial allele is expected to become overrepresented in the genomic regions surrounding the advantageous variant [58,74].

To assess whether genomic segments containing the Duffy locus showed deviations from genome-wide ancestry proportions in the Khwe genomes, we inferred local ancestry while excluding the SNPs that determine Duffy allelic variation from the analysis.

Although local ancestry inference based on genotype array data may have reduced power to detect shorter ancestry tracts generated by older admixture events, we identified a 390 kb region on chromosome 1, encompassing the Duffy gene, where Bantu-related ancestry was significantly elevated (75%) compared to the genome-wide average (58%). This difference corresponds to 2.68 SD above the mean, placing the region in the 99th percentile of the genome-wide distribution of Bantu ancestry (Figure 3A and D; Table S7).

An opposite and stronger deviation was observed in the same genomic region for resident southern African-related ancestry (9,3 vs 30, 5%), corresponding to 3.79 SD below the mean (Figure 3B and E; Table S7).

The asymmetry between the two ancestry deviations appears to be driven by a non-significant local elevation of East African ancestry (Figure 3C and F; Table S7), which could reflect random noise or the presence of FY*B^ES linked markers in the Somali population.

Figure 3. Local ancestry proportions along chromosome 1 and genome-wide distributions of local ancestry proportions in the Khwe. (A-C) Average local ancestry proportions in 114 Khwe individuals along chromosome 1 for Bantu (A), southern African (B), and eastern African-related (C) ancestries. Solid horizontal lines indicate genome-wide average ancestry proportions, while dashed horizontal lines correspond to the 99th and 1st percentiles of local ancestry distributions. Dashed vertical lines mark the position of the Duffy locus. (D–F) Genome-wide distributions of average Bantu (D), southern African (E), and eastern African-related (F) ancestry proportions. Dashed vertical lines indicate the local ancestry estimate for the RFMix window containing the Duffy locus. Kwangali, Juǀ’hoan and Somali were used as proxies for Bantu, southern African and eastern African-related ancestries, respectively. The Chromosome 1 region from position 119429575 to position 144929697 was excluded, due to the absence of markers on the Affymetrix Axiom Genome-Wide Human Origins Array. Coordinates are based on the GRCh37/hg19 version of the H. sapiens genome.

A similar excess of Bantu-related ancestry, accompanied by a reciprocal decrease in southern African-related ancestry, was observed when Hadza and Nyaneka populations were used as alternative eastern African and Bantu-related source populations, respectively (Figures S6-S8; Table S7). Moreover, both signals remained robust across different admixture times and phasing methods (Table S7).

Taken together, these results support the hypothesis that the FY*B^ES allele and its linked haplotypes underwent post-admixture positive selection in the Khwe, following the arrival of Bantu-speaking populations in the lower Okavango River Basin.

Additional genomic regions showing significant local ancestry deviations in at least three sets of source populations are listed in Table S8. Two regions displaying excess southern African–related ancestry contained genes associated with metabolic and lipid-related functions (e.g. ARL15 and LDLR), while one region with excess eastern African–related ancestry encompassed TMTC1, a gene involved in calcium homeostasis.

3.5. Strength of selection

Finally, we applied a deterministic model assuming an infinite population size to estimate the selection coefficient (s) required to reach the observed FY*B^ES allele frequency. This analysis was based on the admixture scenario shown in Figure 2A, using different sets of initial allele frequencies and admixture times (Figure S9; Table S9).

Assuming that selection acts only on individuals homozygous for FY*B^ES (h = 1), the estimated selection coefficients in the Khwe ranged from 0.035 to 0.114. Under an additive model (h = 0.5), the estimated s values ranged from 0.059 to 0.179 (Figure S9; Table S9). We additionally estimated selection coefficients separately for ǁAni-Khwe and non-ǁAni-Khwe populations (Table S9). Inferred selection was consistently stronger in ǁAni-Khwe than in non-ǁAni-Khwe across the tested scenarios (h = 0.5: 0.11–0.30 vs. 0.05–0.16; h = 1: 0.06–0.17 vs. 0.03–0.11). This pattern may be consistent with ecological differences between the two groups, including the occupation of comparatively drier environments by non-ǁAni-Khwe populations.

Overall, our estimates are comparable to those reported in previous admixture studies in Madagascar and Cape Verde [22,23], suggesting that strong recent selection favours Duffy negativity in different epidemiological and ecological contexts.

Ethics Statement

The study involved the collection of saliva samples and ethnolinguistic information from adult volunteers. The study was approved and authorized by the Ministry of Health and Social Services of Namibia (Permit No. RPIV00722019) and Botswana (Permit No. HPRD: 6/14/1) in accordance with local regulatory requirements. All participants were informed about the study in their native language through a bilingual translator, and written informed consent was obtained prior to sample collection.

Consent for Publication

Not applicable.

Availability of Data and Materials

The genotype data from 114 Khwe individuals and 35 Kwangali individuals may be obtained from the corresponding author on reasonable request. The code is available on GitHub (https://github.com/ArtemKasianov/DuffyPaper).

Funding

This work was supported by the Portuguese Foundation for Science and Technology (FCT) [grant numbers PTDC/BIA-GEN/29273/2017, 2022.06259.PTDC, CEECIND/02765/2017].

Competing Interests

Jorge Rocha 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: A-M.F., M.B., G.M., E.C.F. and J.R.; Methodology: A.K. and J.R.; Software: A.K.; Formal Analysis: A.K.; Investigation: A-M.F., M.B., G.M., E.C.F. and J.R.; Writing – Original Draft: J.R.; Visualization: A.K. and A-M.F.; Project Administration: J.R.; Funding Acquisition: J.R. and A-M.F.

Acknowledgements

This study would not have been possible without the administrative and logistic support of the TwinLabs established between BIOPOLIS-CIBIO and the Universities of Namibia and Botswana. We further extend our gratitude to the Ministries of Health of the governments of Namibia and Botswana for granting us research permits, as well as to the local authorities and institutions that supported our fieldwork. We specifically wish to acknowledge the Teemashani Trust representing the ǁAni-Khwe community at Mohembo (North-West District, Botswana), as well as Thaddeus Chedau, Mmoloki Mogomotsi, and Kelatlhilwe Gamaxo Moses who assisted us during the data collection. Finally, we thank Mark Stoneking for comments and suggestions on a previous version of this manuscript.