{'fraction_non_alphanumeric': 0.03838942307692308, 'fraction_numerical': 0.02608173076923077, 'mean_word_length': 5.25673033538878, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genetic structure and sex‐biased gene flow in the history of southern African populations

Abstract Objectives We investigated the genetic history of southern African populations with a special focus on their paternal history. We reexamined previous claims that the Y‐chromosome haplogroup E1b1b (E‐M293) was brought to southern Africa by pastoralists from eastern Africa, and investigated patterns of sex‐biased gene flow in southern Africa. Materials and methods We analyzed previously published complete mtDNA genome sequences and ∼900 kb of NRY sequences from 23 populations from Namibia, Botswana, and Zambia, as well as haplogroup frequencies from a large sample of southern African populations and 23 newly genotyped Y‐linked STR loci for samples assigned to haplogroup E1b1b. Results Our results support an eastern African origin for Y‐chromosome haplogroup E1b1b (E‐M293); however, its current distribution in southern Africa is not strongly associated with pastoralism, suggesting more complex demographic events and/or changes in subsistence practices in this region. The Bantu expansion in southern Africa had a notable genetic impact and was probably a rapid, male‐dominated expansion. Our finding of a significant increase in the intensity of the sex‐biased gene flow from north to south may reflect changes in the social dynamics between Khoisan and Bantu groups over time. Conclusions Our study shows that the population history of southern Africa has been complex, with different immigrating groups mixing to different degrees with the autochthonous populations. The Bantu expansion led to heavily sex‐biased admixture as a result of interactions between Khoisan females and Bantu males, with a geographic gradient which may reflect changes in the social dynamics between Khoisan and Bantu groups over time.

2007; Wood et al., 2005) reflects a long history of population movements and interactions. The so-called Khoisan populations are the descendants of some of the earliest humans inhabiting the region; they are or used to be foragers and pastoralists who speak indigenous non-Bantu languages characterized by the heavy use of click consonants. We use the term "Khoisan" without any assumption about their genetic or linguistic unity (cf. Barnard, 1992). Three language families are recognized among Khoisan (Supporting Information Figure S1): Tuu and Kx'a, spoken by populations known to have practiced hunting and gathering until recently, and Khoe-Kwadi, spoken by a large number of different ethnolinguistic groups practicing diverse subsistence strategies (Güldemann, 2004(Güldemann, , 2005(Güldemann, , 2008(Güldemann, , 2014Güldemann & Elderkin, 2010;Heine & Honken, 2010). Genetic data revealed that the Khoisan populations harbor some of the earliest branching mtDNA and NRY lineages Barbieri et al., 2014Barbieri et al., , 2016Rosa & Brehm, 2011;Tishkoff et al., 2007). Additionally, autosomal genetic data indicate complex patterns of ancestry for most Khoisan groups, reflecting substantial admixture with other groups as well as between different Khoisan groups (Montinaro et al., 2017;Pickrell et al., 2012;Pickrell et al., 2014;Uren et al., 2016).
It has been shown that there are at least two different sets of related episodes of gene flow in the history of Khoisan populations that could have contributed to their current genetic ancestry. The earlier admixture event involves a migration from eastern Africa that occurred 900-1,800 years ago (Montinaro et al., 2017;Pickrell et al., 2014;Schlebusch et al., 2017) and is supported by several independent lines of evidence from different disciplines. Archeological data support an introduction of pastoralism from eastern to southern Africa (Mitchell, 2002;Pleurdeau et al., 2012), while based on linguistic data it has been hypothesized that Khoe-Kwadi-speaking populations are the descendants of these pastoralist migrants from eastern Africa (Güldemann, 2008), where livestock is present from 4,000 years ago (Deacon & Deacon, 1999;Phillipson, 2005). Genetic evidence of shared ancestry between the Khoe-Kwadi-speaking Nama pastoralists and the ǂKhomani and Karretjie (whose heritage languages belonged at least in part to the Tuu family), and East African groups, specifically the Maasai, was observed in autosomal data (Schlebusch et al., 2012). Recent studies of ancient DNA from skeletal remains from Africa demonstrated that all modern-day Khoisan groups for which there are genetic data have been influenced by 9%-22% genetic admixture from East African/Eurasian pastoralist groups (Schlebusch et al., 2017;Skoglund et al., 2017).
Further evidence for a migration from East Africa comes from elevated frequencies of an East African lactase persistence allele in southern African pastoralist groups and in Khoe-speaking groups, particularly the Nama (Breton et al., 2014;Macholdt et al., 2014;Macholdt, Slatkin, Pakendorf, & Stoneking, 2015;Schlebusch et al., 2012), and mtDNA haplogroup L4b2. This haplogroup is found in southern African Nama and ǂKhomani San, as well as in high frequency in the Hadza and Sandawe from Tanzania (Knight et al., 2003;Tishkoff et al., 2007), who also make use of click consonants in their languages.
Additional support for a demic diffusion from East Africa is based on the distribution of Y-chromosome haplogroup E1b1b (E-M293). It has been suggested that this haplogroup spread through Tanzania to southern-central Africa via a movement of people who brought pastoralism 2,000 years ago (Henn et al., 2008), independently of the migration of Bantu-speaking peoples. New studies of autosomal data suggest that pastoralism, after being brought to southern Africa from eastern Africa via at least some degree of demic migration, was subsequently spread within southern Africa mostly via cultural diffusion (Montinaro et al., 2017;Uren et al., 2016).
The more recent admixture event reconstructed with genomic data is a consequence of the Bantu expansion that started around 5,000 years ago from the current Cameroon-Nigeria border. This expansion is one of the most influential demographic events on the African continent (Grollemund et al., 2015), and led to a sex-biased pattern of admixture between Bantu speakers and the local groups already present in territories settled by Bantu-speaking groups, including forager populations such as the rain forest hunter-gatherers in Central Africa and the Khoisan groups of southern Africa (Destro-Bisol et al., 2004;Schlebusch et al., 2011;Tishkoff et al., 2007;Verdu et al., 2009Verdu et al., , 2013Wood et al., 2005). This sex-biased pattern is the result of mating practices that typically involve Bantu males and autochthonous females, but hardly ever involve autochthonous males and Bantu females (Destro-Bisol et al., 2004). This results in the flow of Bantu Y-chromosomes into autochthonous communities, or of autochthonous mtDNAs into Bantu communities, depending on where the children are raised.
In this study, we use previously published mtDNA and NRY sequences collected from a large and comprehensive sample of Khoisan and Bantu-speaking populations to investigate (1) the genetic history and structure of southern African populations, with a focus on their previously undescribed paternal genetic history as well as the distribution of haplogroups and specific lineages within different ethnolinguistic groups; (2) the link between the NRY haplogroup E1b1b (E-M293) and pastoralism; and (3) the intensity of the sex-biased gene flow between incoming non-autochthonous and autochthonous peoples.

| Samples
We collected data from two available datasets which analyzed the same population samples. The Y-chromosome dataset consists of 900 kb sequences from the non-recombining region from 547 individuals belonging to 24 different populations (Supporting Information Table S1 and Figure 1; Barbieri et al., 2016), and the mtDNA dataset comprises complete mtDNA genome sequences from 680 individuals belonging to 26 different populations . The NRY sequences from the neighboring Khwe-speaking ǁAni and Buga populations were merged together into a combined ǁAni/Buga population due to their low sample sizes; the mtDNA sequences were similarly merged to be directly comparable to the NRY dataset. We refer to the dataset of 23 populations (17 Khoisan and six Bantu) that overlap at the level of populations between the NRY and mtDNA datasets from Namibia, Botswana and Zambia, as the "NBZ dataset" (Supporting Information   Table S2). Approximately 93% of the individuals included in the NRY dataset are also included in the mtDNA dataset, while 73% of the individuals included in the mtDNA dataset are also included in the NRY dataset. For analyses of autochthonous genetic structure before nonautochthonous haplogroups arrived in the area, we excluded from the NBZ data sets (both mtDNA and the NRY) all Bantu populations, nonautochthonous haplogroups from Khoisan populations, and Khoisan populations with sample sizes <8 after removal of individuals with nonautochthonous haplogroups. This filtering resulted in 13 populations overlapping between mtDNA and the NRY, and we refer to this dataset as the "AU-NBZ dataset" (Supporting Information Table S2). The analysis of the intensity of the sex-biased gene flow (ISBGF) included additional data from southern African populations for which both mtDNA and NRY haplogroup frequencies were previously published (Barbieri, Butthof, Bostoen, & Pakendorf, 2013;Coelho, Sequeira, Luiselli, Beleza, & Rocha, 2009;de Filippo, Heyn, Barham, Stoneking, & Pakendorf, 2010;Henn et al., 2011;Marks et al., 2015;Schlebusch et al., 2011), and we refer to this as the "SA dataset" (Supporting Information Table S3). The highly admixed Karretjie and Colesberg Coloured populations from South Africa were treated as Khoisan populations in this analysis as they are likely to be partly descended from Khoisan populations (see Schlebusch et al., 2011;Traill, 1996).
Since the time of sample collection, additional linguistic research on the Kx'a family has revealed that the language formerly referred to as ǂHoan consists of three dialects: N!aqriaxe, ǂHoan, and Sasi. The language is nowadays referred to as ǂ'Amkoe (Gerlach, 2016;Güldemann, 2014). Although the samples included under the name ǂHoan mainly stem from N!aqriaxe speakers and include only a few ǂHoan speakers, for ease of comparison with previous studies of these samples we continue to refer to them as ǂHoan speakers.
Most of the individuals assigned to major haplogroup E1b1b (n = 59) were assigned to E1b1b1b2b (E-M293; n = 36), while due to lower coverage the rest were assigned to E1b1b1b2 (E-CTS5487, n = 23). Manual inspection confirmed that the individuals carrying the derived allele for M293 are randomly distributed in all sub-branches within the E1b1b sequence-based network, and thus we assumed that even those individuals where the position was not successfully genotyped are derived for marker E-M293 (the marker genotyped by Henn et al., 2008). Individuals assigned to E1b1b were genotyped for a set of 23 STRs using the PowerPlex ® Y23 System (Promega, Mannheim, Germany; Supporting Information Table S4A) as described previously (Barbieri et al., 2016). To place the southern African samples in a broader picture and search for possible connections with eastern Africa, they were subsequently merged with publicly available STR datasets for the E1b1b haplogroup from Africa (Berniell- Lee et al., 2009;Gomes, Sánchez-Diz, Amorim, Carracedo, & Gusmão, 2010;Henn et al., 2008;Tishkoff et al., 2007), resulting in the "E1b1b-STR dataset" (Supporting Information

| Data analysis
We used previous haplogroup assignments for both the mtDNA and the NRY data (Barbieri et al., , 2016Barbieri, Vicente, et al., 2013). A network for the NRY sequences was generated previously by Barbieri et al. (2016), but was here analyzed by population for the first time to investigate the distribution of specific lineages within different ethnolinguistic groups. Branches in the NRY network were numbered according to the tree in Figure 1 of Barbieri et al. (2016). Additionally, Network 5.0.0.1. (Fluxus Engineering, http://www.fluxus-engineering.com) was used to visualize the relationships between the STR haplotypes genotyped for haplogroup E1b1b and previously published for B2a (Barbieri et al., 2016). We excluded DYS385 a/b, as the PCR analysis co-amplifies these two loci (so alleles cannot be definitively assigned to a locus). The STR-based networks were constructed by either applying the reduced-median method first and then the median-joining method (network with 10 STRs; Bandelt, Forster, & Röhl, 1999), or just the median-joining method (networks with 21 STRs) with STRs weighted according to their mutation rate (Heinila, 2012 (Excoffier, Laval, & Schneider, 2005) for FIGURE 1 (a) Map of the approximate location of the populations and their NRY haplogroup composition. Population labels are color-coded according to linguistic affiliation as indicated in the lower left corner of the map. Green shades within the pie charts indicate haplogroups that are traditionally considered as Khoisan-related, purple shades indicate Bantu-related haplogroups, yellow indicates E1b1b (thought to be East African), blue are Eurasian haplogroups, and gray represents B2a, which is probably associated with both Khoisan and Bantu populations. The Bantu-related, Eurasian, and E1b1b haplogroups are here defined as non-autochthonous haplogroups, while Khoisan-related haplogroups are defined as autochthonous. The yellow area indicates the Kalahari Desert. (b) Schematic representation of the NRY phylogenetic tree, from Barbieri et al. (2016) [Color figure can be viewed at wileyonlinelibrary.com] both the NBZ and AU-NBZ dataset. For the NRY calculations, we used a Tamura and Nei (TrN) model (Tamura & Nei, 1993) with a gamma value of 155.8, and for the mtDNA, a TrN model with a gamma value of 0.26, as suggested by results obtained with jMo-delTest 2 (Darriba, Taboada, Doallo, & Posada, 2012). An AMOVA was performed in parallel for both the mtDNA and NRY to investigate previously proposed groupings and to explore possible factors that might influence genetic structure. Detailed information about the AMOVA groupings and populations included in each of the groups can be found in Supporting Information Table S5. Ecotype information was obtained from the WWF Terrestrial Ecoregions of the World (TEOW) map (Olson et al., 2001) and populations were grouped accordingly. Nonmetric multidimensional scaling (MDS) analyses were performed with the function "isoMDS" from the package MASS (Venables & Ripley, 2002). Neighbor-joining trees (NJ trees) depicting population relationships were generated from the matrix of Φ ST distances with the function "nj" of the package ape (Paradis, Claude, & Strimmer, 2004). Correspondence analysis (CA) was performed with the package ca (Nenadić & Greenacre, 2007), using haplogroup frequencies as input variables. A Mantel test was performed between genetic (Φ ST ) and geographic distances with the R package ade4 (Chessel, Dufour, & Thioulouse, 2004); geographic distances between populations were averaged over GPS data from the individual sampling locations with the function rdist.earth of the package fields (Nychka et al., 2017).  Table S3). For this, we sampled mtDNA haplogroups and Y-chromosome haplogroups that were present in the population independently with replacement, calculated for each marker the frequency of autochthonous haplogroups and subsequently the ISBGF. Population sizes were kept constant throughout the bootstrap analysis. It is not clear if NRY haplogroup B2a is autochthonous to southern Africa or was brought by the Bantu expansion (or possibly both; see Barbieri et al., 2016); we therefore treated haplogroup B2a as ambiguous and ran analyses in parallel treating it as either autochthonous (AUT) or non-autochthonous (NAUT).   (Barbieri et al., 2016;Batini et al., 2011;Marks et al., 2015;Underhill et al., 2000;Wood et al., 2005). Even though macrohaplogroups A2 and B2b are also found in Central African rain forest hunter-gatherers, different subclades are characteristic of Khoisan populations of southern Africa (Batini et al., 2011). In our dataset Khoisan-related haplogroups are, as expected, observed in higher frequencies in Khoisan than in Bantu populations (Supporting Information Figure S2). However, these haplogroups are also observed in relatively high frequency (>14%) in two Bantu-speaking populations (Herero and Kgalagadi; Supporting Information Table S6).
The A2 (A-V50) haplogroup has a narrow distribution in the central Kalahari (Figure 1 and Supporting Information Table S6), and it has been detected in low frequency in Baka foragers from Cameroon and Gabon (Batini et al., 2011). In the network of A2 sequences (Supporting Information Figure S3), all but one haplotype of the Tuuspeaking populations occur on branch 9 (see Table S1 of Barbieri et al., 2016 for more information about branches), while Kx'a and Khoe-Kwadi populations are present in most of the other branches.
All !Xuun, Ju|'hoan South and Ju|'hoan North haplotypes occur in the same branches as Naro haplotypes (for reasons of simplicity from here on, we refer to Ju|'hoan North and Ju|'hoan South together as Ju|'hoan; and to Taa West, Taa North, and Taa East together as simply Taa). The ǂHoan haplotypes are more closely related to haplotypes from neighboring Taa and G|ui populations than to haplotypes from other Kx'a speakers.
Haplogroup A3b1 (A-M51) shows strong regional and linguistic clustering. It represents a major autochthonous NRY haplogroup in most Khoe-Kwadi-speaking populations (ranging in frequency from 0% to 43%; Supporting Information Table S6), and is the only autochthonous haplogroup in the Nama. Interestingly, almost all haplotypes from Khoekhoe speakers (Nama, Haiǁom and Damara) and !Xuun are found in a single branch of the network (branch 13; Supporting Information Figure S4). Another branch (branch 11; Supporting Information Figure S4) almost exclusively contains haplotypes found in the Khoe-Kwadi-speaking Khwe (ǁXo and ǁAni/Buga) and Naro, which inhabit the Okavango delta and neighboring Ghanzi District, respectively. The Ju|'hoan and Taa populations make up the majority of the haplotypes found in branch 18, and they harbor similar yet distinct sublineages within this branch.
Present in almost all Khoisan populations (except ǁXo, Gǁana and ǂHoan), B2b (B-M192) has frequencies higher than 23% in Ju|'hoan and Tshwa, and frequencies lower than 15% in the other populations (Supporting Information Table S6). All but two haplotypes from the Kx'a-speaking Ju|'hoan and !Xuun belong to branch 26 in the network (Supporting Information Figure S5), which they share with all Naro haplotypes and one Haiǁom haplotype. All Taa haplotypes are found in two distinct sublineages within branch 26 (Supporting Information Figure S5).
Given its ambiguous origin, we ran analyses in parallel treating it as either autochthonous or non-autochthonous. The highest frequency of this haplogroup (80%) is found in the Gǁana (Figure 1 and Supporting Information Table S6). In addition, B2a is also found in three
The network for haplogroup E1b1a + L485 (E-L485) sequences (Supporting Information Figure S7) shows a star-like pattern, suggestive of population expansion, that harbors haplotypes from all language families. The Damara haplotypes are found within branch 38, and they are similar to the haplotypes found in neighboring West Bantu-speaking Himba, Herero and Owambo. Interestingly, a newly described sub-lineage (Barbieri et al., 2016) of this haplogroup (branch 35) is present exclusively in Khoe-Kwadi-speaking groups, namely Khwe (ǁAni/Buga and ǁXo) and Haiǁom. The similarity of Damara and West Bantu haplotypes is also noticeable in branch 39 of the haplogroup E1b1a1 (E-M58) network (Supporting Information Figure S7).
The network for haplogroup E1b1a8a (E-U175) sequences is similar to that of haplogroup E1b1a + L485 in exhibiting a star-like pattern with haplotypes from all of the language families in the core (Supporting Information Figure S8). Finally, haplogroup E2 (E-M75) is found in frequencies lower than 5% in Taa, Damara, Haiǁom, Owambo, and Kgalagadi, while in Mbukushu and Tswana it reaches almost 17% ( Figure 1 and Supporting Information Table S6 and Figure S9).

| E1b1b (E-M293) haplogroup
This haplogroup is considered to have an East African origin, and it has been associated with the spread of pastoralism from East Africa to southern Africa (Henn et al., 2008;Trombetta et al., 2015). The sequence-based network of this haplogroup shows a star-like pattern with all language families represented in the core of the network (Supporting Information Figure S10A). Interestingly, most of the Khoekhoe-speaking Nama haplotypes are found in the core. Haplotypes found in the ǂHoan individuals are on a branch shared with Taa and G|ui haplotypes (dotted circle in Supporting Information Figure S10A). The analysis of this haplogroup based on STR data and its possible link to the spread of pastoralism is discussed in detail below.

| AMOVA
Factors that might influence the genetic structure of southern African populations can be explored by grouping populations using various criteria and then examining the proportion of genetic variance shared among groups, among populations within groups, and within populations, using AMOVA (Excoffier et al., 1992). When all Khoisan populations are analyzed together as one group in the AMOVA, they show slightly higher differentiation between populations for mtDNA than for the NRY (17% vs 15%; Table 1), which may reflect geographically structured mtDNA lineages and recent expansion of Banturelated NRY haplogroups. However, this pattern varies when the Khoisan language families are analyzed separately; the Khoe-Kwadi harbor the biggest proportion of between-population variance for both uniparental markers. Overall, Bantu groups harbor levels of between-group variation that are comparable to Khoisan populations for mtDNA but lower for the NRY (Table 1), in keeping with other evidence that the Bantu speakers incorporated more female than male lineages from autochthonous populations, and that Khoisan populations incorporated more male than female Bantu lineages (Destro-Bisol et al., 2004;Marks et al., 2015;Schlebusch et al., 2011;Tishkoff et al., 2007;Verdu et al., 2013;Wood et al., 2005).
The highest between-group variance for the NRY, and one of the highest for mtDNA, is seen when populations are grouped into five groups defined by distinct, geographically organized autosomal ancestry components inferred from unsupervised population structure analysis (Uren et al., 2016). Based on the distribution of autosomal ancestries they defined, Uren et al. (2016) concluded that the autosomal structure in Khoisan populations reflects the role of geographic barriers and the ecology of the greater Kalahari Basin. In order to test if our dataset of uniparental markers is in agreement with this conclusion, we obtained ecotype information from the WWF Terrestrial Ecoregions of the World (TEOW) map (Olson et al., 2001) and grouped populations accordingly (WWF, Table 1). The much lower among-group variance for this grouping (for both uniparental markers) suggests that the ecoregions in which the populations currently live does not explain the genetic structure of the studied populations.
Grouping populations by the geo-linguistic categories defined by Barbieri et al. (2014) results in capturing 19% of the variance seen in mtDNA between eight geo-linguistic categories (Geo-linguistic, Table 1) and is thus one of the best groupings for the mtDNA, but it captures just 9.4% variation for the NRY. When only linguistic criteria are applied, grouping by the two major linguistic groups (Bantu and Khoisan) better explains the variance between groups for both mtDNA and the NRY than when grouping by the four language families (Bantu, Kx'a, Tuu, and Khoe-Kwadi).
The results of the AMOVA analysis for the AU-NBZ data set, that is, after removing the Bantu populations and non-autochthonous uniparental lineages (Supporting Information Table S7) reveals that there is larger genetic differentiation among groups for the NRY than for the mtDNA (23% vs 8%), in contrast to the full dataset (17% vs 21%), suggesting differences in male vs. female migration between Khoisan and Bantu groups, a hypothesis that is analyzed in more detail below. However, although the AU-NBZ dataset was constructed to try to investigate Khoisan population structure prior to the arrival of Bantu populations, the surviving autochthonous genetic structure was undoubtedly influenced by contact with Bantu populations, that is, probably many autochthonous lineages were subsequently lost.    Table S11). Population symbols and colors indicate linguistic affiliation, as shown in Figure 1 [Color figure can be viewed at wileyonlinelibrary.com] continuity correction p value = 0.16). Overall, our data thus do not provide support for a link between haplogroup E1b1b and current pastoralists or Khoe-Kwadi speakers.

| E1b1b STRs
To further investigate the relationships of southern and eastern African E1b1b Y-chromosomes, we genotyped our E1b1b samples for 23 STR loci and merged these with previously published data indicates that the DYS389I haplotype with 10 repeats is most likely ancestral (Figure 3b), contrary to previous suggestions (Henn et al., 2008). However, these findings should be taken with precaution as there is no absolute agreement between the STR-and sequence-   in explaining the variation among groups for the NRY and second-best for mtDNA for both the entire dataset (NBZ dataset,

| Fine-scaled population structure
There are several populations that exhibit substantial frequencies

| Damara, Himba, and Herero
The Damara speak a Khoe language closely related to Nama, yet their current genetic makeup (Montinaro et al., 2017;Pickrell et al., 2012Pickrell et al., , 2014Uren et al., 2016) appears to reflect shared ancestry with Bantuspeaking populations (such as Himba and Herero), making it plausible that they have undergone language shift from a Bantu to a Khoekhoe language (Oliveira et al., 2018). These findings are supported by our data, as they appear to be genetically more similar to Bantu groups (in particular, the Himba and Herero) than to other Khoisan groups in all analyses (Figures 1-2, and 4 and Supporting Information Figure S11).
Other Khoe-Kwadi populations reported to have high autosomal Banturelated ancestry (e.g., Shua, ǁXo, and to a lesser extent, Tshwa; Pickrell et al., 2014) also tend to be more similar to Bantu populations based on uniparental haplogroup composition and analyses based on Φ ST distances ( Figure 2 and Supporting Information Figure S11 and Table S9).
This high level of Bantu-related ancestry, reflecting extensive admixture and/or language shift, needs to be taken into account when considering the relationships of these Khoe-Kwadi populations.

| Haiǁom and Nama
The Haiǁom and the Nama speak closely related Khoe languages. Barnard (1992:12) suggested that the Haiǁom originated from !Xuun speakers who shifted to the Khoekhoe language. This hypothesis is supported by the close similarity of their NRY sequences (branch 13 in Supporting Information Figure S4), and the low NRY Φ ST differentiation between the Nama, Haiǁom, and !Xuun (Figure 2 and Supporting Information Figure S11 and Table S9)

| ǂHoan
The ǂHoan language (nowadays referred to as ǂ'Amkoe, Gerlach, 2016;Güldemann, 2014) is geographically isolated from the other Kx'a languages with which it is related. Even though ǂHoan are linguistically related to Ju|'hoan and !Xuun, autosomal data showed their strong relationship to neighboring Taa populations (Montinaro et al., 2017;Pickrell et al., 2012Pickrell et al., , 2014Uren et al., 2016). The putative genetic relationship with other Kx'a speakers might be preserved in uniparental markers. However, network analysis shows that NRY sequences found in the ǂHoan are related to sequences found in the Tuu-speaking Taa populations and in the G|ui (Supporting Information Figures S3 and S10). The closer relationship of the ǂHoan to neighboring Khoe-Kwadi and Tuu-speakers than to other geographically more distant Kx'a-speakers is additionally supported by mtDNA Φ ST distances, as they appear more similar to neighboring Khoe-Kwadispeaking Gǁana and Naro, and Tuu-speaking Taa North and Taa West, than to other Kx'a speakers (Supporting Information Table S9). The genetic data thus do not provide clear evidence for a genetic relationship with other Kx'a speakers; rather, they support a scenario of longterm extensive contacts between the ǂHoan and surrounding groups, in good accordance with linguistic data (Gerlach, 2016).

| Naro
Naro speak a Khoe-Kwadi language. Their kinship system has been described as a simplified Khoe kinship system with some Ju|'hoan features, and it has been hypothesized that they may have spoken a Kx'a language in the past and subsequently shifted to a Khoe-Kwadi language (Barnard, 2016;Güldemann, 2008). If the hypothetical language shift in Naro was not purely a cultural process, we would expect to see genetic relationships with both Kx'a and Khoe-Kwadi speakers.
Based on both uniparental lineages, the Naro appear to be closely related to Tuu and Kx'a speakers (Supporting Information Table S9, Figure S3, branch 26 in Figure S5, dotted circle in Figure S21, dotted circle and arrow in Figure S22, haplotypes within L0d1c haplogroup indicated with arrows 2-6 in Figure S23, and haplotypes within L0d2a1 in Figure S24), which is in good accordance with previous findings based on autosomal data (Pickrell et al., 2012). However, possible genetic evidence of contact between proto-Naro, who may have spoken a Kx'a-related language, with Khoe-Kwadi populations that could have contributed to the language shift in Naro, may be found in branch 11 of the network for NRY haplogroup A3b1 (Supporting Information Figure S4). This branch harbors almost exclusively haplotypes from current Khoe-Kwadi speakers (i.e., Naro, ǁAni/Buga, and ǁXo, who all belong to the West Kalahari Khoe sub-branch; Supporting Information Figure S1). Although most mtDNA lineages in the Naro are shared with Kx'a-and Tuu-speaking groups, there are some lineages that are shared with Khoe-Kwadi-speaking groups, that is, one L0k lineage shared with Khwe and Gǁana (arrow in Supporting Information Figure S21), one L0d1c haplotype predominantly shared with G|ui and East Kalahari Khoe speakers (arrow 1 in Supporting Information Figure S23), and for L0d2ab there are haplotypes shared or in close proximity to haplotypes from Haiǁom and G|ui (Supporting Information Figure S24). Both the mtDNA and the Y-chromosome evidence thus suggests that the putative language shift in the Naro was accompanied by some gene flow.

| Gǁana and Kgalagadi
The Gǁana are multilingual, like most Khoisan populations; they speak both their own language (in this case, a Khoe-Kwadi language) and a Bantu language (in this case, Kgalagadi). Barnard (1992) (Figures 1 and 4 and Supporting Information Figure S11 and Table S9). Furthermore, they are close to the Kgalagadi in the NRYbased CA plot (Supporting Information Figure S11B). The cultural belief of extreme sex-biased admixture is also supported by the high frequency of NRY haplogroup B2a in Gǁana (80%) and the fact that the lowest pairwise Φ ST values for the NRY between Gǁana and any other population is with the Bantu-speaking Kgalagadi (Supporting Information Table S9). Additionally, STR analyses show that there is less diversity for B2a in the Gǁana than in Bantu populations (Supporting Information Figure S6b). This mixed ancestry for the Gǁana is also reflected in the autosomal data, as they have an estimated 30%-41% Bantu-related ancestry (Pickrell et al., 2014;Uren et al., 2016). Botswana as recently as 350 years ago (Kiyaga-Mulindwa, 1993;Segobye, 1998), Barnard (1992)

| Admixture with East African migrants
Archaeological and linguistic evidence has suggested a pre-Bantu migration from eastern Africa that brought pastoralism and Khoe-Kwadi languages (Blench, 2006;Güldemann, 2008;Mitchell, 2002;Pleurdeau et al., 2012) and thus had a significant impact on southern Africa. Although studies of different genetic markers support the demic migration model from East Africa, they vary in their conclusions concerning the impact and importance of this migration (see Introduction). For instance, the mtDNA data support very limited gene flow Uren et al., 2016), while the lactase persistence and limited NRY data undeniably support the demic diffusion model with significant population movement (Breton et al., 2014;Henn et al., 2008;Macholdt et al., 2014Macholdt et al., , 2015Schlebusch et al., 2012).
One possibility is that the signal of East African ancestry in Khoisan populations was shaped by a heavily male-mediated migration from eastern Africa (as previously proposed in Barbieri et al., 2014).
Thus, it is crucial to investigate the paternal history of Khoisan populations in order to differentiate between the spread of pastoralism due to limited demic migration with more significant cultural diffusion versus a heavily male-biased demic migration from the east that brought pastoralism. We do not find a higher frequency of the E1b1b NRY haplogroup, previously associated with the spread of pastoralism (Henn et al., 2008), in the pastoralist Nama (Supporting Information  Figure 3 and Supporting Information Figure S10a,b). This is contrary to the expectation that diversity and frequency should be highest in the population that most probably brought pastoralism, or represents a direct descendant of such a population. If E1b1b was indeed associated with a migration of pastoralists, subsequent demographic events (e.g., migrations, drift) and/or changes in subsistence practices in southern Africa have diminished the association.
However, even though it cannot be directly associated with pastoralism, haplogroup E1b1b clearly has an eastern African origin. The close relationship of southern African E1b1b STR haplotypes to haplotypes from the two eastern African foraging populations, Hadza and Sandawe (Figure 3), indicate a common origin of haplotypes in southern and eastern Africa. This is in agreement with the linguistic hypothesis of a relationship between proto-Khoe-Kwadi and Sandawe (Güldemann & Elderkin, 2010

| Bantu expansion
Our results confirm and extend previous findings concerning the enormous impact that the spread of iron-using agro-pastoralist populations speaking Bantu languages had on the genetic landscape of southern Africa Beleza et al., 2005;Coelho et al., 2009;Marks et al., 2015;Oliveira et al., 2018;Pickrell et al., 2012;Skoglund et al., 2017;Wood et al., 2005). Arriving in southern Africa 2,000-1,200 years ago (Phillipson, 2005 (Lee, 1993). Conversely, the flow of maternal lineages from forager to Bantu groups occurs if the children are brought up in the Bantu communities. Strong sociocultural taboos inhibit unions between forager males and Bantu females (Cavalli-Sforza, 1986;Lee, 1993). Although this expected sex-biased signature is stronger in Khoisan than in Bantu populations, there is considerable variation in the intensity of the sex bias among different populations ( Figure 5c and Supporting Information Figure S12c), and so other factors must also play a role.
Residential practice is another cultural trait that is likely to influence the distribution of genetic variation, and thus the signal of sexbiased gene flow. Patrilocal populations are expected to show less population differentiation for the mtDNA than the NRY because of higher rates of female migration between local groups, while the reverse is expected for matrilocal populations. Indeed, it has been observed that mtDNA is more geographically structured than the NRY in matrilocal populations (Oota, Settheetham-Ishida, Tiwawech, Ishida, & Stoneking, 2001;Bolnick et al., 2006;Kumar et al., 2006), while patrilocal populations show contrasting patterns (Kumar et al., 2006;Langergraber et al., 2007;Oota, Settheetham-Ishida, Tiwawech, Ishida, & Stoneking, 2001;Wilder, Kingan, Mobasher, Pilkington, & Hammer, 2004). Even though our data set consists of patrilocal Bantu populations and Khoisan populations that preferentially practice patrilocal postmarital residence after an initial period of matrilocality (and to a lesser extent neolocality; Barnard, 1992), we observe larger differences among populations for mtDNA than for the NRY (Table 1), which is the pattern characteristic for matrilocality and not for patrilocality. This deviation could be explained by a rapid male-dominated Bantu expansion over huge geographic areas and incorporation of already geographically structured mtDNA lineages into expanding Bantu populations. Marks, Levy, Martinez-Cadenas, Montinaro, & Capelli, (2012) showed that even though female migration is more frequent among patrilocal populations, males migrate preferentially at longer distances than females, suggesting that patrilocal residence is expected to mostly impact geographically close groups. As geographic distances between populations in our data set are mostly >200 km, which would be considered long-range distances (Marks, et al., 2012), it is possible that the observed pattern is due to a higher migration rate of men at longer geographic distances rather than an overall higher migration rate. However, our data are insufficient for separating the effects of migration rate versus geographic distance; further studies are needed.
One of the most striking findings of our study is the increase of sex-biased admixture from north to south. With the assumption that the initial contact between Bantu and Khoisan populations occurred in the north, the increasing ISBGF toward the south (Figure 5a,b and Supporting Information Figure S12a,b) could be interpreted as an indication that the initial contact involved less sex bias. The north to south increase in autochthonous uniparental lineages in Bantu populations (Supporting Information Figure S16) could suggest either more intense intermarriage with local populations in the south, and/or a gradual accumulation of autochthonous uniparental markers in Bantu populations toward the south as a result of interactions between the southwards migrating Bantu populations and autochthonous populations that they encountered on their way. The former seems more likely, as autochthonous lineages incorporated in Bantu populations tend to be regionally specific, for example, L0k is found mainly in the north, both in Khoisan and Bantu groups, but not in the south Bantu groups (Barbieri, Vicente, et al., 2013;Schlebusch, Lombard, & Soodyall, 2013), in contrast to what would be expected if autochthonous lineages had gradually accumulated during the southward migration of Bantu-speaking groups. This is compatible with the Static and Moving frontiers model that implies limited assimilation of foragers into agro-pastoralist populations during the initial expansion into a new territory, and increased likelihood of gene flow between populations with the occurrence of a static frontier once the carrying capacity of the new territory has brought the expansion of the agropastoralists to a halt (Marks et al., 2015). It is also consistent with Destro-Bisol et al.'s (2004) model that proposes the gradual establishment of social inequalities between agro-pastoralists and foragers, resulting in asymmetric gene flow between them. During the early phase of contact between agro-pastoralists and foragers, the survival of the newly arrived food producing societies in an unfamiliar habitat is probably heavily dependent on the knowledge of autochthonous foragers, and as such their contact might be more egalitarian. It has been suggested that Bantu languages of the Kavango-Zambezi transfrontier region that have incorporated click consonants might have acquired them during this more egalitarian phase, when Khoisan populations had higher social status (Pakendorf, Gunnink, Sands, & Bostoen, 2017). Later, as Bantu groups became more established and less dependent on local knowledge, they became more socially dominant, and hence sex bias increased in intensity due to the establishment of sociocultural taboos (Destro-Bisol et al., 2004). Under this interpretation, the geographic pattern would reflect an increase in sex bias due to changes in social interactions over time.
In conclusion, we have carried out a comprehensive populationlevel study of matrilineal and patrilineal lineages in southern African populations, integrated within their historical and anthropological background. Discrepancies found between the linguistic and genetic relationships of Damara, Naro, ǂHoan, and Haiǁom suggest probable language shift and/or extensive contact between these and other, linguistically unrelated, populations. We find support for a migration from eastern Africa but do not find an association of NRY haplogroup E1b1b with pastoralism today, suggesting that the arrival of pastoralism was more complex than previously suspected. Our study indicates that the Bantu expansion was probably a rapid, male-dominated expansion, during which local Khoisan females were much more likely to be absorbed into Bantu populations than Khoisan males. We