A catalogue of early diverged contemporary human genome variation reveals distinct Khoe-San populations

Domain: neuroscience
Published: 2026-02-10
Editorial outcome: major_revision

Abstract
Creating a catalogue of early diverged genome variation is critical to determine the true extent of human diversity and associated medical impact. Generating deep whole genome data for 150 Khoe-San (12 groups, 1 unclassified), and 40 regionally comparative Southern Africans (3 groups), we identify ~30 million small-to-large variants - over 1.3 million unknown single nucleotide variants. Representing shared traditionally forager lifestyles and click-speaking languages, we identify San and Damara as separate phylogenetic lineages, contributing two admixture waves to Nama. While San represented modern humans’ deep divergence (~115 thousand years ago), Damara divergence is recent, with both showing high effective population sizes between 45–150 thousand years ago. Developing an assembly-based test we report 1,376 genes under positive selection (dN/dS=19.46) of which 479 are significantly associated with forager peoples and, therefore, maintained ancestral alleles that differ from derived genetic variation observed in non-African biomedical resources.

Introduction
Together with the current global references of human populations, ethnically well-defined large-scale high-depth whole genome sequencing (WGS) efforts are required to make sense of population-specific medically relevant variation. The urgency to move away from a singular all-African model has more been recognised with the Nigerian 100K Genome Project. Here we focus on the Namibian population isolates clustered through their shared click languages (Khoi-San) and assimilated forager-herder lifestyles (Khoe-San), representing significant genetic distinction from the widely distributed Bantu and other global populations.

Linguistically, Khoi-San represents a broad group of click speakers, including the indigenous Kw’a (Northern Khoi-San), Tuu (Southern Khoi-San) and Khoe-Kwadi (Central Khoi-San) languages. While Kwadi is extinct, the Khoi-languages are clustered further as Kalahari, including in this study, the Kxoe (specifically Khwe) and Naro, and Khoekhoe (Nama-Damara language cluster and Hai//om). Recruited within the borders of Namibia, the 150 Khoe-San participants in this study represent the following language groups (1 unknown): Kw’a sub-branches, Ju (Ju|’hoan, !Xun and Makaukau, n=38) and ǂ’Amkoe (ǂHoan, n=2); Tuu sub-branch, Taa (!Xoon and ’Nǀoha, n=2); Khoe sub-branches, Kalahari Khoe and Naro (Khwe and Naro, n=18) and Khoekhoe (Hai//om and Nama-Damara, n=89).

Focusing on the cultural Khoe-San identifier, San stems from the Nama-derived ethnonym Saan used to distinguish traditionally forager from herder click-speakers, with largely unfavourable connotations. The anthropologist Alan Barnard describes how San interchanged with Bushman has come in and out of academic/political fashion over the decades. Those living in remote communities within the Nyae Nyae Conservancy of north-eastern Namibia, which include the Ju|’hoan, some !Xun and a few displaced Hai//om, have all elected to continue the traditional and legal practices of hunting and gathering. After 15 years of community engagement, Hayes and co-authors have observed a steady decline in natural resources, and in turn an increase in small-scale subsistence farming and reliance on tourism and government assistance. Concerned that their culture is being lost to the next generation by westernisation, the Nyae Nyae participants have expressed pride in their forager ways.

Outside the Ju|’hoan protected lands, the older generation of !Xun, Makaukau, ǂHoan, !Xoon, ’Nǀoha, Naro, Khwe and Hai//om (the latter sharing the Khoekhoegowab language with Nama-Damara) recall the times when hunting was common practice (now forbidden), while gathering of local plants remains a common source of nutrition and as such share the Bushman (or San) collective. Culturally, the Nama-Damara speakers distinguish themselves as more southerly pastoral Nama, while the northerly distributed Damara refer to themselves as ǂNūkhoen (meaning black) due to their darker complexions, yet mixed hunter-gatherer and herding practices. The Damara participants and community leaders in north-western Namibia speak of an ancestry that pre-dates the arrival of the Bushmen, staking claim to the roughly 5000-year-old hunter-gatherer depicting rock paintings scattered across the region, which pre-dates archaeological evidence for early pastoralism into the region roughly 2000 years ago, and reference locally as Black Bushmen. As such, the Damara distinguish themselves from the Nama-Pastoralist, who in turn distinguish themselves from the Bantu-Pastoralist.

To date, WGS data have been generated for 144 Southern African Khoe-San (82 at low and 62 at high coverage) representing seven linguistic groups from four countries, including South Africa (5 Karretjie, 3 Khomani and 82 Nama), Angola (5 !Xun), Botswana (5 |Gui & ||Gana, 15 Ju|’hoan and 15 !Xoon) and Namibia (5 Nama and 9 Ju|’hoan). To improve on the catalogue of early diverged genomic variation and its associations with forager population structure and positive selection signature of local adaptation, we double the number of available Khoe-San genomes, including 150 individuals (7 related) representing eleven ethno-linguistic groups scattered across the expanse of Namibia. Providing further contextual relevance, an additional 40 study participants were included to represent three distinct cultural, linguistic and/or genetic Khoe-San-related population identifiers. The remote semi-nomadic pastoral Himba peoples of the North-western Kunene region of Namibia (n=10), sharing a language with the Herero, were included for their geographical and reported genetic relationship to the Damara. Representing more traditional lifestyles, the participants live in remote family groups and self-identified as culturally distinct Himba, Tjimba and Twa. In contrast, archaeological evidence and linguistic phylogenies suggest the Southern Bantu peoples (n=12) followed an easterly Bantu migratory route arriving in modern-day South Africa between 1500–1800 years ago. Supported by contemporary and ancient DNA, non-Bantu hunter-gatherer genetic ancestry first appears in the Southern Bantu gene pool roughly 1,300 years ago. Lastly, the arrival of European settlers and East-Indian slaves at the southern tip of Africa in the mid 1600s gave rise to South African Coloured (SAC). In this study, SAC and European Afrikaans (n=18) were included, with the latter tracing back to the 1st generation of European settlement at the Cape. Appreciating that the objective is not a study of Khoe-San prehistory, the inclusion of African data from the Simons Genome Diversity Project (SGDP), representing nine populations (22 genomes), allows for a snapshot of comparative analyses across Sub-Saharan Africa.

Results
Whole genome sequence quality assessment
We generated an average of 41-fold sequencing coverage (range 33–55 fold; Supplementary Data 1). A total of 2.7×10^13 bases of DNA sequencing data were generated, consisting of an average of 139.2×10^9 bases of sequence per individual. For a typical individual, 97.16% of the bases in the reference genome were covered with at least 5-time sequencing coverage. All alignment data tested for sample mix-ups based on the average Euclidean distance of 1,484,044 single nucleotide variant (SNV) markers for each pair of samples showed no mix-ups at most (score 0.2–0.3; Supplementary Data 2 and Supplementary Fig. 1b). Only 11 pairs of individuals had mix-up scores less than 0.2 (range 0.139–0.199); together with pedigree information, pairwise kinship coefficients and identical by descent, consensus results removed seven samples with 1st-degree relationships (full-sib or parent-offspring; kinship coefficient, 0.148–0.257; Supplementary Fig. 1a).

Stringent consensus datasets for SNVs and indels (insertions and deletions smaller than 50 bases in length) were generated using two callers (see “Methods”), and subsequent joint-called variants were passed through machine learning variant quality score recalibration and genotype refinement with the 1000 Genomes Project to filter out low-quality calls (Supplementary Information 2.4). Then, resulting call sets were evaluated in several ways (Supplementary Fig. 1c, d and Supplementary Data 3). First, we compared variants for samples genotyped in duplicate using WGS, together with either the Human Omni1-Quad genotyping array (mean alternative single-base allele concordance 0.9388) or Complete Genomics Inc (concordance 0.95). Second, we compared consensus genotypes to those obtained using the Genome Analysis ToolKit (GATK) informatics tools (concordance 0.939; range 0.937–0.945). These reproducibility estimates indicate the high quality of the genotype calls and effectiveness of machine-learning-based quality filters.

For genomic rearrangements, we joint-genotyped using lobSTR across our custom hg38 repeat panel, with motif lengths ranging from 1–6 bp (Supplementary Fig. 2a), to generate the high quality of short tandem repeats (STRs) after removing three outliers with call rates outside three times the standard deviation of their mean (mean±SD=246,101±27,364 calls) and choosing STRs with reference alleles smaller than 105 bp due to short read lengths (Supplementary Fig. 2b). Furthermore, we used MELT to call ~12 thousand mobile element insertions (MEIs) after filtering based on evidence of discordant and split reads at the breakpoint, erroneous calls and low complexity regions. The joint-called MEIs across all samples encompassed 12,245 Alu, 1,772 LINE1, and 772 SVA elements, with the highest breakpoint resolution accuracy observed among them at 96.26, 87.02, and 60.49%, respectively (Supplementary Fig. 2c). Three main families of Alu elements (AluJ, AluS and AluY) were detected, with the single AluJ (Jo), five AluS (Sc, Sg, Sp, Sq and Sz) and 11 AluY subfamilies (Ya, Yb, Yc, Yd, Ye, Yf, Yg, Yh, Yi, Yj and Yk; Supplementary Fig. 2d).

A unique catalogue of Khoe-San genome variants
Focusing on the 150 contemporary Khoe-San genomes, after excluding for relatedness, 143 remained (87 female, 56 male), with further lifestyle classification as forager (92.3%), including 69 San (10 groups; 1 unknown) and 63 Damara, or pastoralist (7.7%), 11 Nama, while acknowledging the San Hai//om (n=9), Damara (n=65) and Nama (n=11) share a singular Khoi-language. Using the GATK workflow, we identified 30,411,977 single nucleotide variants (SNVs), insertions and deletions (26,427,613 SNVs and 3,984,364 indels), as well as 153,712 short tandem repeats and mobile element insertions (141,254 STRs and 12,458 MEIs) (Table 1 and Supplementary Data 4–6). Overall, 4.81% and 12.04% of these SNVs had yet to be described in dbSNP build 155 and TOPMed data, respectively with significantly more unknown SNVs identified in San (median 33,911 versus 13,120 variants per genome; t=−21.657, P-value<2.2×10^−16; Table 1). San participants revealed more small- and large-sized genome variants than those of Damara (median 7,641,488 versus 7,244,602 variants per individual; t=−20.597, P-value<2.2×10^−16), which was greater for most representative San populations compared to the Damara and Nama (Fig. 1b and Supplementary Data 4–6) and other regionally relevant neighbours sequenced in this study (Supplementary Fig. 3a). However, San showed an increase in homozygosity compared to Damara (median 1,558,203 versus 1,372,128 homozygous SNVs, t>−11.74, P-value<2.2×10^−16; Run of homozygosity, median 49,240 versus 34,714 kb). This would be the result of isolated populations, such as hunter-gatherer societies with traditionally controlled family sizes and more recent population declines. Notably, San maintained 21 ancient and inactive AluJo mobile elements, more than that observed for the Damara and West African Niger-Congo representative genomes (Supplementary Fig. 3b and Supplementary Data 6). These highly diverse elements are thought to have integrated into primate genomes ∼65 million years ago, supporting the preservation of ancestral sequences within peoples having preserved a largely forager subsistence.

Degrees of accumulated or derived (recently arisen) variants were estimated within Khoe-San using known ancestral sequences inferred from human and nonhuman primate references (Supplementary Information 4.2). The total number of derived variants with neutral and deleterious mutational effects based on Genomic Evolutionary Rate Profiling (GERP) scores was significantly higher in San over Damara (median 1,891,580 versus 1,770,072 variants; t=−16.585, P-value<2.2×10^−16) and other regionally relevant populations (t>−19.006; P-value<5.381×10^−12), Furthermore, significant differences were observed for those coding and noncoding variants that were derived and assessed via GERP > 4 or Prediction score = 1 as being deleterious (Fig. 1c, Supplementary Fig. 3c and Supplementary Data 4).

Khoe-San deep structure and phylogenetic analysis
Given that we have introduced 11 Khoe-San groups, seven of which were unreported in the previous African WGS efforts, we performed stringent DNA structure and principal component analyses to contextualise understudied populations. There is a wide consensus that the most divergent group of all modern people are found among the San based on deep population splits using whole genomes and deepest-rooting maternal L0 mitogenomes. Population structure analyses of 7,003,530 consensus SNVs across 316 individuals from three genome sequencing studies (this study, SGDP, 1000GP) optimally separated five major ancestral groups, representing Niger-Congo African, European, Asian, San and Damara (Fig. 2 and Supplementary Fig. 4). The latter shared the ancestry with Herero-Himba at K=4 and showed a gradient of San ancestry. On average, the Nama individuals carried 40% of the ancestry unique to the Damara and Herero-Himba groups, which was also found in relatively high levels in !Xun (7–34%), Khwe (46–59%) and Hai||om individuals (14–99%). The division of the San and Damara forager groups was in line with ADMIXTURE clustering, with highly significant correlations of admixture fractions between methods (Pearson correlation; R=0.76 and 1.00, P-value=1.21×10^−81 and 3.59×10^−13; Supplementary Fig. 5a, b).

Focusing on the top 100 heterozygous loci for each variant type, SNV, indel, STR, and MEI (Alu, LINE1 and SVA), while all clearly defined San substructure, only the highly diverse SNVs and STRs confirmed separation between the three Southern African population identifiers (Supplementary Fig. 6a). Note that each of the SNV and STR markers analysed had mean heterozygosity ranging from 0.32–0.50. Indels had the mean heterozygosity of 0.27, while MEI markers (Alu, SVA and LINE1) ranged from 0.07–0.08.

Furthermore, we used uniform manifold approximation and projection (UMAP) and principal components analysis (PCA) to characterise genetic relationships between Southern Africans (this study) and published ancient and present-day Africans (3,318 individuals in total). The UMAP project of 15 principal components among the merged dataset of WGS and genotyping arrays reflected only a major separation between African and Eurasian geographic regions (Fig. 3a), whereas its subset of the merged WGS data (2783 whole genomes) showed a distinct illustration of the complex structure and population relationships in a large and multi-ethnic dataset. Fine-scale structure was observed between geographic regions (Africa, Europe, Central South Asia and East Asia) and even at the smaller scale between African populations, including Central African, West African, Eastern Bantu, Damara/Herero-Himba and San (Fig. 3b). To better represent genetic distance within Khoe-San, the PCA of all variant types (small consensus variants and high-confidence genomic rearrangements; see “Methods”) also revealed similar fine-tuned relationships for two forager groups that were mostly differentiated from neighbouring populations outside of Namibia (Fig. 3c, d). We recommend the importance of WGS resources to further fine-tune genetic relationships between Khoe-San groups and global populations.

Combining with the above-mentioned DNA structure and PCA analyses, Maximum Likelihood phylogenetic analysis of 137,927 coding exons concatenated across 13 Khoe-San and 12 non-Khoe-San genomes was performed, with individual inclusion based on low admixture fractions (median 0%, range 0–14%, except for 22% in one Southern Bantu). Addressing the phylogenetic analysis for recombination effects between samples, robust DNA substitution models (GTR+ Optimisation of substitution rates+GAMMA model of rate heterogeneity) were independently estimated in each step (105 bootstraps), varying with sample-specific recombination rates. To optimally reconstruct a true tree, we pruned for rogue samples with less confident tree positioning that rendered insufficient phylogenetic signals caused by complex human histories. The reconstructed tree showed two monophyletic clades with 100% bootstrap value (BV) support, using Altai Neanderthal as an outgroup (Fig. 4). The first clade clustered all ten San representative peoples representing five diverse populations, supporting genetic distinction of this forager group as the early diverged human group. The second clade consisted of four subclades with highly supported branches: Twa (BV=100%), Tjimba-Himba (BV=89%), Southern Bantu (BV=89%), and Eurasian (BV=100%). The Southern Bantu clustered close to Damara foragers and Herero-Himba pastoralists with Eurasians falling within them (BV=100%). Damara was subdivided into multiple subclades (range, BV=53–94%), confirming population separations between San and Damara despite their broad Khoe-San classification. Additional evidence of the results using the Maximum Likelihood tree of genome-wide SNVs (8,491,748 markers) across 199 Khoe-San and regionally relevant genomes is also provided in Supplementary Fig. 6b.

Admixture history and dating in Khoe-San
Rare doubleton variants likely to inform recently shared ancestry across 149 Khoe-San genomes (143 our study and 6 SGDP) showed greater between-population allele sharing among San groups over Damara and Nama (Fig. 5a). Compared with regionally relevant genomes for possible admixture between populations, statistics of the form D (Chimp, X; A, B) that quantify the extent to which the allele frequencies of X are closer to those of A or B confirmed historically and culturally related populations of Ju|’hoan, !Xun and Damara, while adjusting various allele frequencies from West African Niger-Congo (Esan, Mandenka, Mende and Yoruba; Fig. 5b, c). In contrast, frequency data based on WGS suggested an excess sharing of the Damara with Herero-Himba speakers relative to Ju|’hoan and !Xun at any derived allele frequencies inferred from the ancestral sequences (Ancestor in D-statistics; Fig. 5d, e).

Haplotype-based admixture dating further suggests that contact occurred between the ancestors of the Damara and San populations within the past 1860 years, while the Nama reveals best-guess contributions from both Khoe-San groups during two different migration waves (Fig. 6a, b). Using 35 San and 21 Damara genomes with the identified single ancestral structures as reference datasets, the two ancestries were the major contributing sources of their respective populations, with ancestral proportions of 0.71–0.83 reported during the two waves. The estimated dates of European admixture in Damara and Nama (30–150 years ago) is consistent with the recent settlement of European populations in the region.

Early and late divergence timing of Khoe-San peoples
When San Ju|’hoansi were compared with other Africans and non-Africans, the divergence time from a common ancestor of all present-day humans was noted approximately at 115 kya using MSMC2 (mutation rate 1.25×10^−8 per base pair per generation, 30 years per generation; Fig. 7a), consistent with other studies. Having included most participants with no admixture (median 0%, range 0–9%, except for a !Xun genome admixed by 14%), single midpoint estimates (RCCR=0.5) are still inadequate for divergence timing, so our results would be to observe a gradual split over time of the curves between the populations using the same method and those obtained by simulating instant split scenarios of 20–500 kya without subsequent gene flow. Interestingly, the divergence of Damara (culturally, ǂNūkhoen Forager) and Yoruba (West African Niger-Congo) dated around 12.7 kya, with their cross-coalescence rate close to zero, which places this event post Eurasian ancestral genetic separation estimated at >40 kya for recent studies. A similar trend observed for Twa (Herero-Himba languages) and Yoruba at around 12.8 kya. Note that all these MSMC2 curves are inconsistent with simulating splits (20–500 kya), suggesting a picture where African divergence times were shaped by ongoing gene flow over tens of thousands of years, as evident until the time as recently as the method can infer. We therefore applied Relate, which incorporates phylogenetic tree-based information to enable inference of effective population sizes in the recent past of a larger number of genomes similar to our dN/dS ratio between populations described above.

Using the Medical Genome Reference Bank of healthy Australians as a control with selectively depleted pathogenic variants (MGRB, n=1144 deep sequenced genomes), we investigated positively selected outlier regions of Khoe-San genomes, using Fisher’s exact test, and their functional potential in human health. We identified 633 positively selected variants to be significantly different from those expected to be randomly selected between the populations compared (P-value<5×10^−5, -log P=7.307–246; Supplementary Fig. 10a, b and Supplementary Data 11). Inferred from the ancestral sequence of 12 primate genomes, the accumulation of derived variants among these variants or a higher derived allele frequency was shifted towards MGRB (also Eurasian) compared to Khoe-San groups. However, mild proportions of derived variants were different among their genomes using neutral markers outside genic and CpG regions (Fig. 8b). The neutral genetic markers are less biased for demographic modelling of complex evolutionary scenarios without selection. Taken together, forager populations in Southern Africa still maintained many ancestral alleles. This finding agrees with a higher accumulation rate of variation in non-Africans (Oceania, America and West Eurasia), supporting the out-of-Africa human migration that would have caused local adaption due to region-specific environmental challenges in the last 60,000 years.

Searching for trait-informative annotation, our significant selected variants contained 508 variants (80.25%) annotated to have functional potential based on sequence and physical properties of amino acids, conserved protein domains and trait association (Fig. 8c). Around 42% (266 variants) were significantly associated with drug responses and traits reported in phenotypic databases, PharmGKB, GPCRdb, NHGRI-EBI Catalogue, and GTEx. Particularly, 235 were correlated with differential gene transcription in particular human tissues, supporting our genome-wide amino acid-based selection scans developed to capture medically relevant diversity. Of all the 2840 positively selected variants, 74.81% (2124 variants) were annotated in the public domains, leaving approximately 25% of the adaptive nonsynonymous variants understudied in the literature. Those significant selected variants representing 479 annotated genes were found to be enriched within an olfactory signalling pathway (17 relevant genes; P adjusted=2.928×10^−5; Supplementary Fig. 10c), with 13/17 genes displaying 21 nonsynonymous variants within coding regions of G protein-coupled receptors (GPCRs; Supplementary Fig. 10d) that play a role in interacting with external stimuli to mediate cellular responses. Interestingly, PPFIBP2 and SYNE2 genes reported to have risk alleles for prostate cancer, with their frequencies higher in Africans than in Europeans, have significant nonsynonymous variants in our study, with one and three amino acid positions under strong positive selection (Fig. 8d). In contrast, an exemplary gene, ACE2 that is critical for SARS-CoV-2 transmission had DNA sequences very conserved with the mean ω value as little as 0.348 (negative selection) across the gene length. All 1,144 MGRB participants had no nonsynonymous changes in this gene, while three nonsynonymous variants were identified among the Khoe-San participants (10 times smaller number; n=111).

Methods
Population inclusion, community engagement and ethics
Irrespective of gender, 190 study participants (genetically determined as 99 female; 91 male) greater than 18 years, with a preference for older community members belonging to a Khoe-San language group (n=150) were recruited within the borders of Namibia (Supplementary Information 1 and https://www.hayeslab.net/videos) as part of the Diversity inclusive Genome Study for Southern Africa (DiGS-SA), a subsidiary of the Southern African Prostate Cancer Study (SAPCS) Namibia, under informed consent following extensive community discussions and continued engagement. Community over individual compensation, including addressing of community-specific needs, was provided post-initial sampling and under the guidance of the community members. Within South Africa, participants were recruited as part of the University of Limpopo Ubuntu Genomics Project (see https://www.hayeslab.net/videos), a subsidiary of the SAPCS and aimed to represent country-specific ethnolinguistic diversity broadly defined as Bantu, SAC and European (n=30). As participants were sourced from volunteered university staff and students, no participant compensation was received. We have previously reported on ancient Khoe-San heritage within South African Bantu and SAC communities while addressing cultural considerations based on population identifier terminologies. In brief, while Khoe-San (or Khoi-San) languages are grouped through their use of click consonants, Bantu includes a large group of non-click languages with Southern African exceptions such as isiXhosa and isiZulu having historically incorporated click consonants. All study participants self-identified back to two generations; final population labels indicated where the person lived, and mostly self-identified, although recognising the mixed heritage.

The research presented was performed in accordance with the ethical standards of the overseeing human research ethics committees (HREC) and local governance, as per the 1964 Helsinki Declaration. In Namibia, the study was reviewed and approved by the Ministry of Health and Social Services (MoHSS 17/3/3 HEAF003 2024 latest update; previously 17/3/3 HAYES 2008 and 17/3/3 HEYNS 2014) and in South Africa by the University of Pretoria HREC (43/2010 2022 latest update, merged with 280/2017 and including US Federal wide assurance FWA00002567 and IRB00002235 IORG0001762) and the University of Limpopo Faculties of Health Sciences (MREC/HS/214/2011). Samples were shipped to Australia in accordance with institutional and governmental Material Transfer Agreements (MTAs) and including the Republic of South Africa Department of Health Export Permit (National Health Act 2003, J1/2/4/2 No 1/12). WGS and analysis was performed in accordance with the approval granted by St. Vincent’s Hospital Human Research Ethics Committee (SVH/15/227).

WGS, data processing and variant calling
DNA extracted from blood underwent 2×150 bp sequencing on the Illumina HiSeq X instrument (Kinghorn Centre for Clinical Genomics, Garvan Institute of Medical Research). Raw sequencing reads were aligned against human reference hg38+alternative contigs using bwa (v.0.7.15). The Genome Analysis Toolkit (GATK, v.3.5-0) was used for duplicate marking, indel realignment and base quality recalibration. Germline small variants were identified using two callers (HaplotypeCaller and FreeBayes). We also processed whole-genome raw reads of 22 African genomes and two Asian genomes reported in a previous publication using our unified pipeline. Full details of the unified analysis are provided in Supplementary Information 2. The discovery of genomic rearrangements (short tandem repeats and mobile element insertions) using WGS data is detailed in Supplementary Information 5.

Familial relatives and quality control
We searched for close family relatives by computing, for each pair of individuals, sample mix-ups (using 1,484,044 SNVs and discordance cut-off at 0.20; Supplementary Data 2), kinship coefficients and identical by descent segments (using 15,603,438 SNVs) and exclusion of 1st-degree relationships (full-sib and parent-offspring). Joint-genotyping variation data suitable for downstream analyses were assessed through different analysis pipelines and genotyping/sequencing technologies for their concordance (Supplementary Data 3). The detailed explanation is provided in Supplementary Information 3 and Supplementary Fig. 1c, d.

Merged data of published present-day human and ancient genomes
We retrieved 380 relevant individuals from 13 present-day African populations genotyped on the Human Origins array and the Illumina Omni bead chip, and included published ancient South African genomes and high-coverage genomes from the Human Diversity Genome Project (HGDP) and 1000 Genomes Project (1000GP), including 2,778 individuals from 67 worldwide populations defined by gnomAD v3.1. After merging with our whole genomes and filtering with a minor allele frequency >0.01 and a missing genotyping rate <0.1, we obtained 93,976 autosomal SNVs across the worldwide populations. Kinship analysis (KING; Supplementary Information 3) was performed to remove 1st-degree relations, with a threshold of 0.177. More details are provided in Supplementary Data 7.

Population structure
In total, 143 Khoe-San and 33 regionally relevant genomes (this study), 28 previously published African genomes (SAPCS and SGDP) and 179 reference individuals from the 1000GP (Africans, African Americans and Eurasians) contained 7,003,530 filtered variants computed for robust population structure analyses (Supplementary Information 6). For ADMIXTURE v1.3.0 analyses of K=2–15 (10 replicates each) using the default five-fold cross validation, the above-mentioned merged dataset of high-coverage genomes was filtered for a minor allele frequency >0.1, individual missing genotype data <0.1 and variant call rate >0.9 and then pruned for linkage disequilibrium in PLINK v2.0 using a window size of 50 bp, a step size of 10 bp, and a r^2 threshold of 0.1, leaving 289,025 SNVs among 2964 individuals. Pong v1.5 was used to check for concordance between replicates for each K and for visualisation.

The non-linear dimension reduction UMAP method provided fine-scale relationships between geography and ethnolinguistic identities in Khoe-San populations, analysing our merged dataset with or without samples with genotyping arrays. Principal component analyses of six groups deeply sequenced (San, Damara, SAC, South African European, Herero-Himba and West African Niger-Congo; Supplementary Data 1) were computed across different variant types (SNV, indel, STR and MEI) as described in Supplementary Information 4, 5.

Phylogenomic analysis
The large-scale Maximum Likelihood tree of 16,083 genes from 33 ancestrally distinguishable genomes identified in this study and the Altai Neanderthal published previously was reconstructed using RAxML, as revealed in Supplementary Information 7. In total, whole-genome diploid consensus sequences encompassing 137,927 exons or 26,039,475 bp were converted from analysis-ready GATK data. Local alignment of those sequences was performed, and an extended majority-rule consensus phylogeny was built by reconciling the best tree with 100,000 bootstrap searches, using Altai Neanderthal as an outgroup. Eight individuals with insufficient phylogenetic signals were removed in a final phylogeny using the RogueNaRok programme.

Patterns of derived variation and doubleton allele sharing
Joint-called consensus SNV data using the two callers described above were used to infer derived variation, with ancestral state inferred from ancestral hg38 genome (Ensembl compara 93 database). Coding and noncoding regions were functionally annotated, with further details provided in Supplementary Information 4. Rare doubleton variants observed twice across the autosomal DNA dataset were identified across the 143 Khoe-San genomes using the biallelic SNV data. Loci that had missing genotypes greater than 10% were removed. We calculated f4 and D-statistics of those derived variants at different frequency intervals (≤0.1, 0.1–0.2, 0.2–0.3, 0.3–0.4, 0.4–0.5, 0.5–0.6, 0.6–0.7, 0.8–0.9 and >0.9) using the ADMIXTOOLS2 programme and reported D values and |Z|>3 as significant. We compared f-statistics between forager populations and Niger-Congo/Herero-Himba speakers to address a relationship between forager populations and admixture from neighbouring populations.

Admixture proportion and dating
We reconstructed the timeframe of admixture waves between major ancestral components in Khoe-San populations using the fast GLOBETROTTER programme. We followed the programme instruction for the analysis of 69 San, 63 Damara and 11 Nama genomes in this study, together with the merged WGS data described above and phased by SHAPEIT2. In brief, ChromoPainter v2 was run to estimate the effective population size and mutation rate of the data (407.590 and 0.000441, respectively). It painted haplotypically phased chromosomes of targeted admixed individuals with the haplotypes of all other individuals from other reference populations (113 Italy, 108 LWK, 18 Southern Bantu and 127 YRI genomes).

Population size and separation
We used the Pairwise Sequentially Markovian Coalescent (PSMC) to infer population sizes using a single diploid genome per population that was scaled with a generation time of 30 years and a mutation rate of 1.25×10^−8 mutations per nucleotide per generation. We used MSMC2 to estimate changes in population split times. This analysis used phased genome datasets and eight haplotypes (two genomes from each of a pair of populations). By expected site-frequency spectrum models using momi2, we also simulated all possible pairwise population divergence times under simple, clean split scenarios (20–500 kya). Diploid effective size and recombination rate parameters were set to 12,000 and 1.25×10^−8, respectively. More information is shown in Supplementary Information 8.

To complement the above-mentioned coalescent-based methods, Relate v1.2.2 was employed using 21 Damara and 30 Yoruba, with no suggestive admixture signals defined in this study to avoid deep coalescence time of San contribution in this tree-based analysis as shown in Supplementary Fig. 6b. Biallelic SNV consensus data (2 callers) were phased using SHAPEIT2 and filtered by the 1000GP accessible genome pilot mask to keep only passed positions. Ancestral/derived alleles for each SNV were defined by the human ancestor GRCh37 e59 before within-population and cross-population coalescence rates were computed using the EstimatePopulationSize.sh script. Using IBDNe, the ancestry-specific effective population size was estimated. We followed the guidelines provided by the authors, using our merged dataset of WGS and genotype data. Local ancestry information was inferred using 20 San, 20 Damara, 38 Southern Bantu and 113 Italian genomes from this study and previous publications as reference source populations.

Genome assembly-based analysis of positive selection
To assess the impact of positive selection on shaping human diversity and its relevance to population fitness, we have developed amino acid-based selection tests using de novo genome assemblies (see the Code availability section). The main steps of this analysis include i) generating an individual’s genome draft; ii) extracting sample-specific target sequences; iii) generating cohort-wise multiple sequence alignment per gene; iv) performing local and amino acid-guided alignment methods; and v) positive selection tests. Full details of analysis steps, parameters and relevant tests are provided in Supplementary Information 9, 10.

Alternative methods to identifying potential SNVs under positive selection in forager populations (admixture fraction <30%; 55 San and 56 Damara individuals) were computed using PBS and XP-EHH statistics. PBS utilises pairwise Fst values of exonic regions among three populations to quantify genetic differentiation along each branch of their corresponding three-population tree using 60 CEU genomes from 1000GP as a control to find outlier PBS scores. Normalised XP-EHH of phasing variant data detected selected variants that were near fixation in one forager group but not in the CEU population using the rehh package in R. Only SNVs with minor allele frequency >0.01 were considered while phased with a reference panel of haplotypes following SHAPEIT2 guidelines for the XP-EHH analysis. To identify candidate regions for selection, we then averaged the XP-EHH over windows of length 100,000 bp with a step length of one kb between windows; the windows chosen were required to have >50 SNVs and represented the top 1% with the highest percentage of outlier scores.

Khoe-San signatures for local adaptation
Medical Genome Reference Bank (MGRB) is a catalogue of deep whole-genome variation across thousands of healthy elderly individuals, mostly of European descent, with no reported history of cancer, cardiovascular disease, or dementia. Comparing this resource acts as a control population when linking to Khoe-San genomes and assists in variant filtering and assignment of pathogenicity for genetic traits (Supplementary Information 10). To identify outliers between the forager genomes (San and Damara, n=111) and the control population, Manhattan plots of positively selected variants were visualised with -log P>4.301, with the analysis of their potential in functional biology via ANNOVAR and significance in trait association through multiple databases (i.e., PharmGKB, GPCRdb, GWAS Catalogue, and GTEx). We generated homology protein models of G protein-coupled receptors, and the modelled protein structure was recorded in the Protein Data Bank format and visualised using the PyMOL Molecular Graphics System. Comparing the derived allele frequency of the positively selected variants with that of neutral variants outside genic regions and CpG islands for each dataset was performed to indicate selective signature for local adaptation.

Ethics and inclusion statement
As per the Health Equity Research Outcomes and improvement Consortium (HEROIC) Prostate Cancer Precision Health (PCaPH) Africa 1K inter-institutional Collaborative Research Agreement (CRA), which includes generating regionally representative Panel of Normal (PoN) genomes, including Namibia and South Africa, and the Global Code of Conduct for research in resource-poor settings, locals have been included in all aspects of the research including study design, local primary ethics approvals and stewardship, study implementation, analysis and authorship, to intellectual property and data ownership. Capacity building across the region includes (i) awarded and self-managed budget allocation to South Africa, which has led to numerous employments, (ii) sourcing infrastructure, resourcing and providing clinical training, (iii) co-supervision and exchanges for postgraduate students to genomic intensive partner laboratories, (iv) providing access to off-site high performance computational infrastructure, while (v) holding local training workshops in genomics. Through engagement and inclusion of local policymakers, local conservancy and public health leaders, the team is committed to the dissemination of scientific data back to communities and local government.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "We commend Jaratlerdsiri et al. on an interesting manuscript, which describes unprecedented diversity in Namibian genomes, and clarifies questions of population structure among Khoe-San-descent populations and neighboring groups. The addition of Damara, Khwe and Hai//om genomes from the Khoe-language family are especially exciting. We have several minor comments regarding nomenclature and clarification (there is some awkward phrasing throughout). We recommend de-emphasizing the estimates of population divergence as derived from MSMC analysis, as these estimates have been shown to be biased by gene flow, sequence coverage, cryptic structure and, as also highlighted here, phasing. Despite the appealing simplicity of MSMC analysis, fully parameterized demographic models are more interpretable (and the authors are missing key genomes from eastern Africa which represents a migration into southern populations). There is also a bit of confusion about using the term \" ǂNu ̄ khoen Forager\" as applied to populations other than the Damara and the ancestry component that is clearly shared among the Himba/Herero/Hai\\om/Nama/Tjimba etc.\nMinor comments:\n— We recommend hyphenating “Khoe-San” throughout the paper and figures when referring to the southern African Khoe-San populations collectively and where southern African ancestry is not referred to as “southern African Foragers”. In the past, members of the South African San Council requested to remain distinct from their Khoe-speaking neighbors, and to be referred to by their specific ethnic identities or generally as the San (Schlebusch, C. (2010). Issues raised by use of ethnic-group names in genome study. Nature, 464(7288), 487–487. https://doi.org/10.1038/464487a). Previously, we have acknowledged their difference in identity through hyphenating and capitalizing ‘Khoe’ and ‘San’ while referring to the collective ancestral component. See also Schlebusch 2020, MBE.\n— Our second recommendation regards the nomenclature around language. We’ve found it easier to refer to the population / ancestry as “Khoe-San”, given the above, and use “Khoisan” to reference the language phylum, consistent with how linguists have previously referred to the larger category.\n-- Line 60. There are 25 Khoe-San genomes in Schlebusch et al (2020, MBE); 84 Nama in Ragsdale et al (2023, Nature); 30 !Xoo and Juu in Fan et al. (2022, Cell); 6 in HGDP. This is a rough total of 145 genomes. --> \"To date, in-depth WGS has included only about 150 Khoe-San individuals representing eight linguistic groups from four countries. We double the number of available genomes, including previously unrepresented linguistic groups.\"\n-- Line 83. Awkward phrasing. --> \"Culturally, the Nama-Damara speakers are distinguished by sheep/cattle/goat herding by the Nama who are presently concentrated in central and southern Namibia, while the northerly distributed Damara refer to themselves\"\n-- Line 98. Remove \"with an association with behavioural modernity\". All human populations are exhibit modern human behavior. Also extended data Fig 9.\n--Line 175. \"genome variations\" --> \"genome variants\" and throughout the manuscript\n--Figure 1B (and in the text). Does \"small\" mean SNV? Are you summing SNVs and indels? Just keep these categories separate for comparison to other papers.\n--Line 178. Awkward phrasing, remove 'colors'. --> \"The ten Southern African Forager populations tend to have greater genomic variation, for all types, as compared to the Nama and Damara populations (Fig. 1b; Supplementary Tables 4-6)\"\n— Layering language families on top of the STRUCTURE plot in Figure 3 is a great idea, and the population column referencing the ancestral component each color represents across the K values will be useful for the audience. Half of the plot though is taken up by standard reference populations. It would be easier to visualize the structure if there were fewer populations. We would recommend removing the Mandinka, Esan, Mende, and African-American samples in the main figure. (No need to rerun things, just exclude these populations and you can still refer to the full suite of ref groups in the Supp.)\n--Line 241. \"Given that Southern African Forager, ǂNu ̄ khoen Forager and West African Niger-Congo were chosen based on little admixture (< 30%) observed in the DNA ancestral analysis described above, only principal components (PC1-3) of small variation showed the latter two groups closely related to each other.\" No idea what this sentence means, or why anyone would exclude putatively admixed individuals from a PCA, which reifies artificial human clusters.\n— Beginning at line 973 of the supplement, “In this study, MSMC estimates of effective population size for a San (!Xun) individual were vastly different around 50 kya between trio-phased and reference-phased data (Ne = 25,906 and 16,616, respectively), but the same for a mild bottleneck in the population history compared to Europeans.” Could you clarify what data was used for the trio-phased SHAPEIT run?\nMajor comments:\n--Abstract. Remove sentence. \"This study suggested the former separated ~115 thousand years ago, while the latter recently diverged 13 thousand years ago from a human ancestor involved in the ‘Out-of-Africa’ event.\" Not a complete logical thought. The Southern Foragers separated from whom 115KYA? The \"Nukhoen Forager\" 13KYA estimate in the paper is divergence from the Yoruba, not from an ancestor involved in the OOA (which long predates 13KYA).\n-- Figure 3. I'm not really understanding the motivation to describe the ancestry frequent in the Damara, Himba/Herero, Hai||om, and Nama as \"Nukhoen Forager\". It's clearly not a forager ancestry as it is shared among both pastoralist and forager populations. There is no evidence from the STRUCTURE analysis that the ancestry dates to a pre-Bantu population (because this is an atemporal clustering analysis), while such an idea is a hypothesis which can be explicitly tested later on in the paper. As the ancestry is most frequent in the Himba and Herero, it cannot be described as a \"phylogenetically distinct click-speaking groups\". It has been previously shown that the Himba/Herero/Damara share a similar genetic ancestry (Uren 2016, Genetics; Oliveiera et al 2018 AJPA) consistent with divergence from a recent common ancestor <2000ya. Please provide a different descriptive name for the Figure and in the text.\n--FIgure 3. In the legend it mentions 7M markers used for analysis, but in the Supplemental Methods description it is a little odd that the authors describe only 79,000 SNPs for the fastSTRUCTURE, which is a pruned version of the 7M SNPs extracted from the 1000G, present, SGDP merged dataset. This substantially fewer SNPs than are used from most array-based structure analyses. The authors have likely thrown away significant information by LD-thinning so extensively? Can the authors clarify these different datasets and which analysis goes with which? What is a \"logistic prior model\" as applied to an unsupervised analysis with K>2 clusters? K=8 is clearly overfitting, as it is pulling out small numbers of individuals assigned to high level ancestries.\n--Figure 3. Is inadequately described in the main text. The authors should note that the Damara sampled here have a gradient of Southern Forager ancestry. The Nama individuals carry ~50% of the green Himba/Herero/Damara/Hai||om ancestry - in contrast to multiple other studies where they have little admixture from this component. This suggests that this pattern is specific to the sampling locale in this study. There is also quite high levels of this 'HHDH' ancestry in the !Xun and Khwe.\n--MSMC/PSMC analyses have been performed quite a few times on both Khoe-San and other African genomes (see Schlebusch et al 2020, MBE among others). The use of the midpoint RCCR was originally intended to estimate divergence in a simple 2-pop isolation model. It is not meant for geographically proximate populations, or those which may have experienced long-term gene flow. The authors appear to acknowledge this: \"single midpoint estimates are still inadequate for divergence timing\" but then go on to use the method. For example, Europeans and Asians split long before 20KYA, and this is a perfect illustration of the misspecification associated with MSMC two-pop divergence approach. Unfortunately, this whole section is misleading. Trying to spin this is a more positive way, what is the main question here?\nOne hypothesis is that the Damara share recent common ancestry with neighboring Bantu-speaking cattle pastoralists, but lost their cattle and language generations ago, adopting cultural practices from the Nama or a related Khoekhoe group. The alternative is that the Damara/Himba/Herero etc belong to a central western African population which was distinct from the ancestors of the Bantu. Adoption of agriculture and pastoralism occurred???\n--One could rephrase the MSMC results around these central questions, such as below:\n[Under a two population, isolation without migration model, the MSMC2 analysis indicates that the Damara genomes have moderately deep coalescent times with the Tjimba, Himba, Herero, Twa (broadly southern Bantu-speaking cattle pastoralists) with an RCCR=0.5 varying between 5-15KYA. Older than one would expect if the group of populations descended from a simple western Bantu expansion circa 4,000 years ago. However, one possibility is recent flow from Southern Forager populations into the Damara (Figure 3) that introduces loci with older coalescent dates. Comparison between the Damara and the Ju|Hoansi genomes illustrates the much deeper coalescent times [120KY], supporting a long period of isolation between the ancestors of the Damara and Khoe-San groups.]\n-- It's not clear to me that a description of the MSMC2 curves can resolve the competing scenarios for the Damara 'origins'. One would need to build explicit models (ABC, moments or other parametric approach).\n--Line 295. (and Ex Data Fig 7). Other PSMC curves show that the Yoruba 'grow' in size after a local minima Ne ~50KYA, the contrast between the Damara and Yoruba in 6D is likely an artifact of the data merging. Further, population structure will result in dips and peaks (Mazet 2016, Heredity). We cannot conclude anything about past 'expansions alongside each other' based on PSMC plots. This whole paragraph should be removed, or revised to merely describe the curves without arguing about why there are relative differences.\n— Beginning at line 337, “To identify positively selected non-synonymous variants that may play a role in local adaptation to lifestyle-related behaviors, we compared two different extremes of lifestyles observed in our forager populations in Namibia versus the Australian Medical Genome Reference Bank (MGRB; n = 1,144 deep sequenced genomes) that represents healthy and well-adapted individuals to Western practices”. It might help to clearly state the hypothesis and expected outcomes of this analysis in the body of the paper. Are you expecting convergence in positively selected loci? At which level? What are the driving environmental features that would impact selection?", "reviewer_2": "In this manuscript W.Jaratlerdsiri et al. (i) reported novel genetic variations including SNPs and CNVs (ii) investigated population structure and (iii) genetic adaption of KhoeSan-related populations using 150 novel high-coverage WGS combined with publicly available sequences of different African groups.\nAmong all results the main novelty of this manuscript results in (i) the dataset itself with new WGS of KhoeSan-related groups who are underrepresented in public databases but very important when working with African populations (e.g. imputation, phasing ...) and (ii) the detection of new signatures of natural selection. However, I suggest the authors to restructure some parts of the manuscript (Introduction and results sections). I am also more skeptical and frustrated about the demographic history section of the manuscript. This part of the paper needs also more context.\nHere are my detailed comments:\nI. Introduction\n• Paragraph 3 of the introduction: from line 73 ‘Observing a steady decline in natural resources, we have seen …’ : Why using « we » ? Is this sentence a quote? Is it an observation that was made by some authors of the manuscript? I suggest reformulating this whole paragraph.\n• Line 98: “behavioral modernity”: I suggest the authors to clarify what they mean by modernity (agropastoralism?)\n• I have difficulties pinpointing clear objectives and aims in the introduction.\nII. Results\n• Whole genome sequence quality assessment: It appears that this section closely resembles the method section, and all figures cited are supplementary figures/tables. For clarity and readability, I recommend restructuring this section to differentiate it from the method section.\n• L144: “Concordance rate of 0.9388 …”. Are the individuals genotyped with the HumanOmni1-Quad a subset of individuals sequenced in the manuscript? If so, to assess the quality of the genotype calling I suggest the authors to calculate the per sample concordance of WGS genotypes versus genotyping array (of the same individual) at intersected sites (Malaspinas et al. Nature. 2016). To assess the quality of the newly generated dataset the authors can also calculate the transition-transversion ratio (TiTv) per sample (that should be ~ 2.1 for WGS).\n• Line 188: ‘again supporting their prehistoric lifestyles’: I respectfully disagree with the authors. I assume that behind ‘prehistoric lifestyles’ the authors refer to foraging/hunting and gathering mode of subsistence. First, I suggest the authors to avoid the use of ‘prehistoric’ as South African Khoe/KhoeSan peoples still practice hunting and gathering and are not the only modern human groups (for example, Hadza and Sandawe from Tanzania, rainforest hunter-gatherers from central Africa or from the Philippines). Second, the fact that Southern African foragers ‘maintained 21 ancient and inactive AluJo mobile elements’ is not related to their lifestyle but to their unique demographic history: (i) Oldest divergence of all modern populations, (ii) isolation as mentioned by the authors, which is associated with (iii) the limited admixture signals with Bantu farmers of some South African foragers (e.g <10% of Bantu ancestry in Taa or Naro Khoisan groups as reported by Pickrell et al. Nature Communication. 2012)\n• From line 203: I suggest moving the paragraph about the D statistics in another part than the current section entitled ‘A unique catalogue of KhoeSan genome variation’.\n• Admixture detection: The authors used STRUCTURE-like approach to quantify admixture in South African foragers. This approach does not formally test for admixture (Lawson, van Drop and Falush. Nature communication, 2018). The authors should benefit from this amazing dataset they generated to address questions related to admixture history (timing, proportions, sources) between South foragers and Bantu-related groups but also to test the presence (and proportions) of other ancestries, such as pastoralist-related ancestry as suggested by wang, Goldstein, Bleasdale et al. Science. 2020. To do so, the authors should use accurate methods: allele frequency-based methods (qpADM,DATES ...) or haplotype-based methods (Chromopainter/Finestructure/Globetrotter or other similar approaches)\n• Line 242 ‘[…] were chosen based on little admixture (<30%)’: 30% of admixture is huge and it is more than enough to bias analyses such as PSMC.\n• Line 280-284 ‘[…] divergence of Damara and Yoruba was very recent at around 12.7kya …. Yet preceding the Southern African Bantu migration event at 4-6 kya’: I am not convinced by this date. MSMC2 can consider migrations but only between the ‘haplotypes’ used as input. As mentioned above, I suggest the authors to deeply investigate admixture history in the region. In addition, although the authors filtered out samples based on admixture proportion, they still retained samples with up to 22% of admixture (median 0%, range 0-22%). Same apply for Twa Herero-Himba.\n• Line 283 ‘[…] genetic separation observed at 20kya consistent with the observation of basal Eurasian …’: In this paper the authors give a lower bound based on the time of the ancient sample they used (>24 kya). I think the separation is >35kya: Fu, Meyer, Stenzel et al. PNAS.2012; Malaspinas et al. Nature.2016; Choin et al.Nature. 2016 (the two last estimated divergence between Sardinians and Han ~40kya).\nA previous study (Pickrell, Patterson, Barbieri et al. Nature Communication.2012) using genotyping array suggests “Damara were a non-Khoisan population who acquired their language from their Khoisan neighbours (the Nama) with little Khoisan gene flow”.\nSchlebusch, Skoglund, Sjodin et al. Science. 2012 also found an East African connexion for the Nama population and ‘suggest that they originate from a Southern San group that adopted pastoralism with some introgression from an East African group that potentially brought pastoralist practices.‘\nOptional: The author can also try IBDNe tool which can infer ancestry-specific Ne (using non-Bantu ancestry masked data) : https://faculty.washington.edu/browning/ibdne.html\n• Positive natural selection: this part of the manuscript is not within my expertise especially the dN/dS statistics. Here is a general comment:\nThe authors found signatures of natural selection in 1,376 genes out of 18,680 tested: Compared to a previous study(Schlebusch, Skoglund, Sjodin et al. Science. 2012, using iHS and Fst-based statistics to detect natural selection), I am surprised by the number of genes under selection detected in the manuscript (high number). Did the authors try other standard approaches to detect natural selection in humans such as iHS/XP-EHH (for more recent period) and Fst-like approaches (for more ancient time frame)?", "reviewer_3": "Jaratlerdsiri and colleagues report about 150 new high-coverage Southern African genomes from Namibia. The coverage of several different groups makes this a highly valuable dataset. Unfortunately, the analysis presented in the paper does not do justice to this unique dataset. I have outlined my major concerns about the data and the analyses below.\n1. Ethics, consent, and sharing: The manuscript does not provide sufficient details on the data collection. Given that these genomes are from highly marginalized populations that have previously been aggrieved by genomics research, it is pivotal that the ethical aspects of the study are evaluated critically.\ni. Firstly, did the Namibian San council approve the study? Proof of approval needs to be provided.\nii. A full section on the recruitment of the participants and the ethics are provided. This should answer questions such as:\na. When were the Kho-San participants recruited? Who performed the recruitment?\nb. What did the participant information sheet/Case report form say?\nc. How were the concepts of whole genome sequencing, and secondary data sharing explained?\nd. Which languages the informed consents were translated into? Could the participants sign? If not, how it was captured?\ne. Did the informed consent have a provision for individual-level data sharing?\nf. How will the access to the data be regulated? Who will oversee that?\niii. Could the authors clarify how the \"extensive and continued engagement\" as mentioned in the manuscript was performed?\na. How many times the participants/communities were visited before the study and how many times they have been contacted after the sample collection?\nb. What was the mechanism and frequency of community feedback\na. The author list includes only two authors from Namibia and none of them seem to have expertise in genetics or history – so how was it ensured that the participants/communities understood the key messages and nuances?\n2. Population structure: The authors had the unique opportunity to investigate the population structure and correlation between genetics and ethnolinguistic identities of individual San populations. Surprisingly the authors did not compare these genomes to any of the public genomics datasets (Schlebusch et al. 2017, Choudhury et al. 2020, Sengupta et al. 2022, Ragsdale et al. 2023, Tallman et al. 2023, Fortes lima et al. 2023) from Southern Africa that could have a provided a framework to interpret the variations and the structure reported in these novel genomes.\nEven in terms of presentation - the 3-dimensional PCA plots with only a few populations are difficult to interpret. Phylogenetic trees with only a subset of genomes were not very informative either. The authors should present their genomes in the backdrop of at least some of the previous datasets and if possible, employ some of the more recent approaches such as UMAP and GCAE to visualize the population structure in a more informative way.\n3. Population divergence and demographic history: A 13,000-year split between ǂNūkhoen and Bantu speakers is one of the highlights of the results. However, this does not align with the current understanding and models of African population migration and history. Moreover, this claim is based only on a single method and concerns have been raised about the accuracy of PSMC/MSMC-based population divergence estimates. Therefore, the authors need to employ other approaches such as SMC++ and RELATE that can analyze a larger number of sequences per group to support their claim. Also, it would be useful to check whether the same results are obtained if other proxy Bantu-speaker populations are considered.\nAs a lot of work needs to be done to support this claim, in case the authors find that it might take a substantial time to do this, they can take this section out.\n4. The comparison of European ancestry Australians to Khoe-San to identify lifestyle-associated selection is perhaps the weakest section of the manuscript.\nThere are several layers of differences – including strikingly different demographic histories (out-of-Africa related and probably tens of strong bottleneck events thereafter in one group), genetic drift, non-random/assortative mating, introgression from groups such as Neanderthal and Denisovans, the geographically stratified impact of pathogens, the impact of different climatic conditions that could have caused these populations to harbor different alleles at different frequencies.\nThe attempt to relate the differences in the populations with such old divergences and extremely different histories to be just due to “forager\" and \"Western\" lifestyle is not meaningful. Had the comparison been performed to related African populations with different lifestyles it might have still made sense.\nAlthough the dN/dS ratio-based selection scan has been in use for a long time it is rarely employed these days for identifying signatures of selection in humans. In case, the authors believe that this is superior to other well-established (allele frequency-based. LD-based, SFS-based) methods then they need to demonstrate it using specific cases and in-depth comparisons with results from other methods. I am not sure if that is within the scope of the current manuscript. Given that the premise they are applying the selection scan to is in itself problematic. I believe this section should be taken out completely.\nMinor comments\n1. The scheme of nomenclature is not uniform – the abstract starts with Khoe-San and the introduction say this is a Khoe-San genomic resource, but then the rest of the manuscript goes on to describe the populations as “Southern African forgers”.\n2. Several statements that are not well-founded and need to be removed or reframed. For example, the paper states that - \"Our forager populations also maintained 21 ancient and inactive AluJo mobile elements (thought to integrate into primate genomes ∼65 million years ago41), again supporting their prehistoric lifestyles.\"\nHow do the authors know that the other populations did not lose these elements due to strong bottleneck or genetic drift but that the Khoe-San retained these due to prehistoric lifestyle? Did they check whether these elements are fully absent in all other humans, especially African populations? Is there a reason that the prehistoric lifestyle is benefited by this?\nIn summary, I would recommend the authors demonstrate that the data was collected, and stored with proper ethics, permissions (especially the San council), and community engagement. They also need to show that they have the necessary permissions from communities/authorities for cross-border sharing of individual-level data.\nThey need to expand their population structure section and bring in other relevant public datasets. I would recommend the authors take out sections on selection. Also, they need to either provide additional data for supporting the PSMC/MSMC-based claims or take out the section too.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Point-by-point response to the reviewers’ comments\nREVIEWER COMMENTS\nResponses/answers (A, blue) to questions (Q, black) have been fully addressed.\nChanges/updates to the ‘main text’ or ‘Supplement’ have been tracked in red.\nReviewer #1 (Remarks to the Author):\nQ 1: We commend Jaratlerdsiri et al. on an interesting manuscript, which describes\nunprecedented diversity in Namibian genomes, and clarifies questions of population structure\namong Khoe-San-descent populations and neighboring groups. The addition of Damara, Khwe\nand Hai//om genomes from the Khoe-language family are especially exciting.\nA: We are very grateful of your comments and insights, while addressing all of them as follows.\nQ 2: We have several minor comments regarding nomenclature and clarification (there is some\nawkward phrasing throughout). We recommend de-emphasizing the estimates of population\ndivergence as derived from MSMC analysis, as these estimates have been shown to be biased\nby gene flow, sequence coverage, cryptic structure and, as also highlighted here, phasing.\nDespite the appealing simplicity of MSMC analysis, fully parameterized demographic models\nare more interpretable (and the authors are missing key genomes from eastern Africa which\nrepresents a migration into southern populations).\nA: The manuscript includes minor changes throughout to ensure we adequately address any\nawkward phrasing, which includes clarity of phrasing within the context of the Abstract, the\nIntroduction has improved phrasing enabled a reduction in word count and textual flow, while\nagain this is addressed for clarity in the Discussion.\nAcknowledging MSMC limitations, we performed additional analyses including\nfastGLOBETROTTER, Relate and IBDNe. See Results and new Figures 6, 7 and\nSupplementary Fig 7.\nWhile we appreciate our study does not sequence east African populations, this does not form\nthe scope of this study but will certainly form part of future Pan-African collaborative efforts\nunderway.\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I would like to start by thanking the authors of the manuscript for thoroughly addressing all my questions and comments.\nBelow, I provide some minor comments for consideration:\n• The statement, “These highly diverse elements are thought to have integrated into primate genomes ∼65 million years ago38, supporting the preservation of ancestral sequences within peoples having preserved a largely forager subsistence”, could benefit from an additional analysis: I suggest exploring whether there is a correlation between the percentage of Khoe-San ancestry across all South African populations (including non-forager groups) and the number of AluJo elements.\n• I recommend replacing terms like “our San” and “our South Africans” with alternatives such as “San genomes generated for this study”, “San genomes generated here”, or “San genomes newly generated”\n• “The division of the ‘San’ and ‘Damara’ forager groups was in line with ADMIXTURE clustering, with highly significant correlations of admixture fractions between methods (Pearson correlation; R = 0.76 and 1.00, P-value = 1.21 × 10⁻⁸¹ and 3.59 × 10⁻¹³; Supplementary Fig. 5a, b). Both ADMIXTURE and STRUCTURE rely on the choice of K to estimate admixture proportions and although there is a high correlation between the two algorithms, they also both suffer from the same issues: (i) a population can appear as a mixture of two components without necessarily reflecting admixture and (ii) there are descriptive analyses like PCA.\nTo address these limitations, I recommend using qpADM to quantify admixture for all individuals. Alternatively, consider SourceFind (https://people.maths.bris.ac.uk/~madjl/finestructure/sourcefind.html) which leverages the output of ChromoPainter to estimate contributions from different sources to a target. SourceFind can be applied at both the population and individual scales.\n• XP-EHH and PBS scans:\nIf possible, I suggest highlighting genes and specifying the number of genes (in a table or supplementary figure) that overlap between these two scans and also with dN/dS.", "reviewer_3": "The authors have addressed several of my concerns in their revised manuscript and response to reviewers. However, the following concerns still need to be addressed:\n1. It is unclear from the response provided, whether the informed consent included provision for sharing data on a public database as well as secondary analysis of this data by external researchers. I can understand that the framework of consent/ethics might not have been as well-defined at the time of initial sample collection as it is now. Also, explaining concepts such as public databases to specific communities might be challenging. However, as this study is based on extremely marginalized communities it needs to be verified that the participants/communities were adequately informed about all the potential risks involved in such data sharing before the consent was sought. An alignment of the sample collection and related documentation with the current ethical standards would also ensure that this valuable dataset becomes and remains accessible to the broader scientific community. Therefore, a paragraph describing the exact information provided to the participant, especially concerning sharing data on public databases, secondary use, and the risks involved, is critically important.\n2. The addition of PBS and XPEHH has improved the section on positive selection. However, I am still worried about the use and interpretation of dN/dS.\nThe use of dN/dS at the genome-wide level for assessing positive selection in human populations is extremely rare in literature and over the years several studies have highlighted the nuances involved in the interpretation of this estimate when derived from a single population. Kryazhimskiy and Plotkin (2008) have pointed out that the interpretation of selection from dN/dS, when applied to a single population and short time scales could be inaccurate as they note - “the hallmark signature of positive selection over divergent lineages, dN/dS>1, is violated within a population. For population samples, the relationship between selection and dN/dS does not follow a monotonic function, and so it may be impossible to infer selection pressures from dN/dS”. This was supported by simulations performed by Mugal and colleagues (Mugal et al. 2014) that show dN/dS to be biased by polymorphisms at short time scales and that “it can take substantial time for the expected value to settle at its time limit where only fixed differences are considered\". This study further demonstrated that \"intrinsic fluctuations of a ratio of stochastic variables, can even under neutrality yield extreme values of dN/dS at short time scales or in regions of low mutation rate\". The review by Del Amparo and colleagues (Del Amparo 2021) notes other possible factors, such as genetic drift in a small population, which could amplify the biases. Moreover, this study reports about 8% of the genes to be under positive selection in a population, which is way higher than previous estimates of the number of positively selected loci in the human genome.\nTherefore, the authors need to convince the reader that they have considered and corrected for possible errors and biases, and their interpretation of all SNPs with dN/dS>1 are due of positive selection is correct. Also, it would be useful if the authors could provide a comparison, showing how many of the signals detected in their assembly-based workflow would have also been found had the authors used the more traditional approach for estimating dN/dS. This could provide an estimate of the extent to which the approach and pipeline employed in this study are contributing to the observation of such a high number of positively selected genes.\n3. In the statement - \"variants of a higher derived allele frequency was shifted towards the westernized MGRB compared to the forager groups\", the authors have ascribed derived allele frequency differences between the San and the MGRB dataset to westernized vs forager lifestyle.\nAssuming modern humans evolved 250 K years ago, all human populations have been hunter-gatherers until 10,000 years ago (96% of human history), and most retained this lifestyle until about 200 generations (~98% of human history). Whether the timeframe of about 200-300 generations would be sufficient to lead to such strong and widespread derived allele frequency differences needs to be supported with data. Moreover, several previous studies have noted such differences between African and European populations, irrespective of lifestyle, and have suggested these might be explained by the OAA bottleneck and adaptation to new environments since then. Unless the authors plan to add further analysis to support this, associating allele frequency differences with lifestyle needs to be avoided.\n4. Unfortunately, the GWAS-like comparison of Khoe-San with MGRB does not add much value. Firstly, the sample sizes of these datasets (23 vs 1144) are not comparable. Secondly, given that these are extremely distant populations unless an adjustment for population structure is made using genome-wide data it is difficult to separate functional differences from background differences. It would be important to know whether the proportion of “positively selected” SNPs that show P-value>10-5 is significantly different from what would be seen if an equal number of SNPs were randomly selected. Instead, of this analysis, the authors could simply present a heatmap showing the frequency of select variants in different gNOMAD/1000 Genomes populations - which could allow readers to not only see how these SNPs vary between non-African and Khoe-San but also between Khoe-San and African groups as well as between different non-African populations.", "reviewer_5": "Although substantial portion of results in the previous version, along with the abstract has been amended, some of the major concerns persist and new ones appear in the revised version.\n1. The section “Early and late divergence timing of Khoe-San peoples” raises a lot of questions. The authors use different three tools to analyse these genomes. The results fromMSMC2 suggesting a much older divergence of Damara and Yoruba doesn’t seem to be correct. As a complementary approach the authors have also run RELATE and the results seem much more reliable (divergence time ~ 2000 years). These two time points don’t corroborate at all. So, what is the message coming out from this analysis, the authors can’t just throw in results from different tools that suggest different results and not discuss a line other than “This recent population separation could, in part, be explained by the larger sample size analysed.” – for the RELATE results. Even though the RELATE divergence time fits better with the understanding of the history of these populations, they are buried deep in the supplementary. I think instead of relying on the MSMC2 results, the authors should run RELATE for all population pairs (Twa and Yuroba etc) which should then be a part of the main Figure7.\n2. I have some concerns about the ancestry specific IBDNe analyses interpretation as well. To understand the results, it is very important we see the TrueNe, without which, it is not possible to assess which ancestry is closer to the original population size trajectory. Also, the sentence “The pre-admixture effective sizes were driven by San, Damara and Bantu contributions, with the latter being the highest among the three studied groups 20-100 generations ago (Figure 7e-g).” is not very clear, could the authors reword the sentence. Sengupta et al., 2021 has shown the analyses line IBDNe can be impacted by cryptic relatedness. So, it would be good to check if the populations in this study have cryptic relatedness. If yes, it will be good to perform a sanity check to see if the results change after using a more stringent cutoff.\n3. In Line 623, “lternative methods to identifying potential SNVs under positive selection in forager populations (admixture fraction < 30%; 55 ‘San’ and ‘56’ Damara individuals) were computed using PBS and XP-EHH statistics.”. I still think ~30% admixture is quite a lot, and both XP-EHH and PBS results can be impacted by this. Signal of selection might be inflated/weakened due to the differential admixture, and hence a careful interpretation of the results is necessary. Could the authors use a more stringent cutoff for these selection scans? Or at least present a caveat it in the discussion section.\nOther comments.\n1. Though the paper presents many results, the discussion section after first few lines is entirely focussed on their newly developed method to identify adaptive genes under positive selection. Other results like admixture dating and population divergence time also need to be discussed and put in context with other studies.\n2. A number of questions raised by the reviewers related to the Damara (and Nama) ‘ancestry’, that could have made the results more reliable/robust were not answered clearly in the rebuttal. Moreover, responses like “This will be our future collaboration with a larger access of Bantu WGS data” do not address the core concerns.\nReviewer1 comments\n⁠16. I agree with the previous reviewer that there is a need to label Damara as a new ancestry just because of the study focus. These population and ancestries have been studied and a standard nomenclature is established. Unless there is a strong data based rationale for it, it would be confusing to add new ancestries.\n18.⁠ ⁠The authors in their response miss to incorporate the critical observation that Nama from this study to show a completely different ancestral composition compared to that observed in other studies. It might be useful to comment on whether this could be related to sampling geography.\n22.⁠ ⁠The authors do not answer the basic questions about the hypothesis they are trying to test and the their expected outcomes for the positive selection analysis. Given different genetic/ancestral backgrounds and sample-size differences the case-control analysis is very difficult to justify. Moreover, there is no attempt to unpack whether the dN/dS differences observed are due to the extremely old split or has anything to do with forager lifestyle."}, "author_response": "Response to reviewers’ comments in Blue\nReviewer #2 (Remarks to the Author):\nI would like to start by thanking the authors of the manuscript for thoroughly addressing all my\nquestions and comments. Below, I provide some minor comments for consideration:\nQ: The statement, “These highly diverse elements are thought to have integrated into primate\ngenomes ∼65 million years ago38, supporting the preservation of ancestral sequences within\npeoples having preserved a largely forager subsistence”, could benefit from an additional\nanalysis: I suggest exploring whether there is a correlation between the percentage of Khoe-\nSan ancestry across all South African populations (including non-forager groups) and the\nnumber of AluJo elements.\nA: We couldn’t afford more main/supplementary figures in this manuscript. We have AluJo\ninformation in Supplementary Table 6. Also, we perform similar analysis of ancestral/derived\nvariants across this study and GNOMAD data in new Supplementary Figure 10a.\n[...]"}, {"round": 3, "reviewer_comments": {"reviewer_2": "All previous major issues have been satisfactorily addressed, and I have no further concerns. Thank you for your efforts in clarifying and strengthening the manuscript.", "reviewer_3": "The authors have addressed my major concerns.", "reviewer_5": "The authors have addressed my concerns."}, "author_response": ""}]

Decision letter
[Editorial Note: See attachment]