African populations speaking Bantu languages (Bantu-speaking populations (BSP)) constitute about 30% of Africa’s total population, of which about 350 million people across 9 million km² speak more than 500 Bantu languages^1,11. Archaeological, linguistic, historical and anthropological sources attest to the complex history of the expansion of BSP across subequatorial Africa, which fundamentally reshaped the linguistic, cultural and biological landscape of the continent. There is a broad interdisciplinary consensus that the initial spread of Bantu languages was a demic expansion^{2,4,5,6,7,12,13,14,15} and ancestral BSP migrated first through the Congo rainforest and later to the savannas further east and south^{2,4,5,7,15,16,17,18,19,20,21,22,23}. However, debates persist on the pathways and modes of the expansion (Supplementary Fig. 1 and Supplementary Note 1) (see refs. ^18,19 for linguistic-based syntheses).

Whereas most recent human expansions involved latitudinal movements through regions with similar climatic conditions^24,25, the expansion of the BSP is notable for its primarily longitudinal trajectory, traversing regions with highly diverse climates and biomes, including the highlands of Cameroon, central African rainforests, African savannas and arid south-western Africa. Despite consensus on its demic nature, genetic studies of the BSP expansion have not revealed the typical serial-founder effect observed when small migrant groups settle in new areas, leading to a decrease in genetic diversity with increasing distance from the putative homeland^13,26,27. This might be a result of subsequent genetic diversity increases from admixture with local populations or long-distance interactions with later Bantu-speaking migrations, known as ‘spread-over-spread’ events^2,20,28,29. This underlying complexity, coupled with the different migration histories proposed by linguistics, archaeology and genetics, makes the expansion of BSP interesting for exploration with newer population genetic methods and modelling approaches that are spatiotemporally sensitive.

Although whole-genome studies of African populations have become available recently^30,31,32,33 and some localized genome-wide genotype studies exist^4,14,15, comprehensive genomic data for BSP across sub-Saharan Africa remain limited. To better understand their large-scale spread, we collected a genome-wide genotype dataset of 1,763 individuals (Supplementary Table 1), including 1,526 Bantu-speaking individuals from 147 BSP across 14 African countries and 237 other sub-Saharan African individuals (Supplementary Fig. 2 and Extended Data Fig. 1). This dataset includes 117 populations not represented in previous genomic studies and covers most main branches of the Bantu language family¹⁸: north-western 2 (2 NW-BSP 2), west-western (7 WW-BSP), south-western (13 SW-BSP plus the Damara, a Khoe-Kwadi speaking population from Namibia sharing a genetic profile with BSP) and eastern (44 E-BSP) (Fig. 1a and Supplementary Table 2). In addition, we generated genomic data for 12 ancient individuals from Late Iron Age sites of south-central and southern Africa (present-day Zambia and South Africa), spanning 97–688 years before present (bp). This comprehensive dataset allowed us to explore the demographic history of BSP using allele-frequency and haplotype-based methods, genetic diversity summary statistics and spatial modelling providing insights into African human history.

After quality control and merging with publicly available data from representative ethnolinguistic groups (Supplementary Information), we assembled a dataset of 4,950 individuals from 124 populations (111 sub-Saharan African and 13 Eurasian populations with at least 10 individuals per population), hereafter referred to as the ‘AfricanNeo’ dataset (Supplementary Figs. 3–5 and Supplementary Table 2). To visually represent genetic affinities between populations, we applied four dimensionality reduction methods (Supplementary Information) and found evidence for fine-scale population structure between sub-Saharan African populations with a clear geographical component and a broad correspondence with the main linguistic groups in Africa (Fig. 1b,c, Extended Data Fig. 2 and Supplementary Figs. 6–13). Population substructure within BSP (NW-BSP 2, WW-BSP, SW-BSP and E-BSP) can be distinguished (Supplementary Note 2).

Population substructure and suggestions of admixture are also apparent in model-based clustering analyses (Fig. 2a, Supplementary Figs. 14–27 and Supplementary Note 2) and show a finer representation of population ancestries with three main BSP-associated genetic components (Extended Data Figs. 3 and 4 and Supplementary Table 3): the dark-green component found in most BSP, the light-blue component shared between non-Bantu Niger-Congo and western BSP (NW-BSP 2, WW-BSP and SW-BSP) and the orange component mainly found in south-eastern BSP.

BSP also show differential genetic affinities with other populations (Fig. 1b,c, Extended Data Fig. 2 and Supplementary Table 4). This pattern may be the result of genetic admixture with local populations, during and/or after the expansion of BSP across subequatorial Africa^{7,13,34,35,36}. We have formally tested the hypothesis of admixture and its regional character using f3 and f4 statistics (Supplementary Information). The results confirm significant and differential contributions of Afro-Asiatic-related ancestry in eastern BSP from Kenya and Uganda, of western rainforest hunter-gatherer (wRHG)-related ancestry in western BSP from the Democratic Republic of Congo (DRC) and the Central African Republic (CAR), and of Khoe-San-related ancestry in southern BSP from South Africa, Botswana, Zambia (Fwe population) and Namibia (Extended Data Table 1, Supplementary Figs. 28–32 and Supplementary Tables 5 and 6). These findings underscore the intricate genetic history of BSP, characterized by distinct admixture patterns with diverse local groups in specific geographic regions of subequatorial Africa (Supplementary Note 2).

To assess if admixture with local groups is the main process driving spatial patterns of substructure in BSP (Figs. 1c and 2a), we masked out admixed genomic regions in BSP³⁷ keeping only west-central African (WCA) genomic components (Supplementary Information). This masked dataset allowed us to minimize the influence of non-Bantu-speaker ancestries in subsequent analyses (Supplementary Fig. 33). Principal component analysis (PCA) on the admixture-masked dataset (Supplementary Note 3) shows that BSP retain a clear genetic structure which aligns with geographic (Fig. 3b,d) and linguistic (Fig. 3a,c) features, suggesting that processes other than genetic admixture influence spatial patterns of BSP diversity. However, this structure could also be driven by outlier BSP with increased genetic drift, for example, Herero and Himba from Namibia largely influence PC2 (Figs 3a,b and Supplementary Fig. 34). We therefore repeated the PCA after excluding Himba and Herero from the analysis and still observed population substructure in the remaining BSP (Extended Data Fig. 5).

We then investigated whether BSP harbour signatures of genetic isolation and population size changes potentially driving the observed structure. Most BSP show similar patterns of genetic drift reflected in runs of homozygosity (ROH) (Supplementary Figs. 35–47, Supplementary Table 7 and Supplementary Note 4) and changes in the effective population sizes (N_e) that support population expansion signatures in the past 10–30 generations (Supplementary Figs. 49–51 and Supplementary Note 5). The Himba and Herero populations notably deviate from the general patterns of BSP, showing higher values of the genomic inbreeding coefficient (Supplementary Figs. 39 and 48), higher averages for the longest ROH length categories (Supplementary Table 7 and Supplementary Fig 41) and higher intensities of founder events (I_f = 1.6% and 1.2%, respectively; Supplementary Figs. 52 and 53 and Supplementary Table 8) than do other BSP (Supplementary Notes 4 and 5). These signatures can be the consequence of genetic isolation since their arrival in south-western Africa or shared cultural practices, for example, endogamic practices linked to cattle herding lifeways³⁸, as suggested by mitochondrial DNA data³⁹, genotype data⁴⁰ and exome sequencing data³⁸. The early twentieth century Herero genocide by imperial Germany is not expected to trigger an increased ROH signal, as this was a relatively recent event in the context of population demographic histories (Supplementary Note 4).

Exploring different models and analysing their fit to the observed genetic data can yield valuable insights into the population history that underlies the genetic patterns and geographic distribution of BSP. A strong correlation between genetic relatedness and geography suggests an isolation-by-distance (IBD) model, which assumes stepwise gene flow between neighbouring groups leading to gradients of genetic affinity across geographic space. Our dataset of BSP fits an IBD pattern (Supplementary Figs. 54–57 and Supplementary Information), including when admixture is removed (Supplementary Figs. 58–61 and Supplementary Note 6), consistent with previous findings based on fewer BSP and a smaller dataset². Alternative models could, however, explain these patterns. For instance, under a serial-founder model we also expect a strong correlation between shared genetic ancestry and geography. However, in contrast to IBD models, a serial-founder model would also show a decrease in genetic diversity from the putative region of origin. To distinguish between these two models, we investigated the spatial distribution of three genetic diversity summary statistics (haplotype richness, haplotype heterozygosity and linkage disequilibrium; Supplementary Information) suitable for array-based genotype data¹². The estimated statistics support a serial-founder model in which the highest genetic diversity is found in western BSP with a steady decline with distance towards eastern and southern BSP (Supplementary Figs. 62–65 and Supplementary Note 7). This pattern is stronger in the admixture-masked dataset (Fig. 4, Extended Data Fig. 6 and Supplementary Figs. 66–75). Further evidence supports serial-founder dynamics during the BSP expansion from west-central Africa; for example, significant demographic founder events have been inferred in 19 BSP (Supplementary Table 8 and Supplementary Note 5) and a maximum-likelihood tree of the admixture-masked BSP dataset shows north-western BSP 2 at the base of the tree and most eastern BSP forming a monophyletic group (Extended Data Fig. 7, Supplementary Fig. 76 and Supplementary Note 8). By contrast, admixture largely drove the shape of the maximum-likelihood trees for the unmasked datasets (Supplementary Figs. 77–80).

Overall, these analyses support the suggestion that the expansion of BSP started from west-central Africa and spread mostly through serial bottlenecks throughout subequatorial Africa. The negative correlation between genetic diversity and distance from the source, even in the unmasked dataset, suggests admixture had a small impact on genetic diversity of BSP, either because there was not that much gene flow with indigenous groups or some BSP moved on before they received substantial local gene flow. The fact that the admixture patterns are largely region-specific (Fig. 2a and Extended Data Fig. 3)—that is, in each population, we detect non-BSP ancestry from local groups and not from elsewhere—suggests the latter.

To gain a deeper insight into how the expansion of BSP unfolded, we investigated the spatial routes and timing of their movements. First, we used a climate-informed spatially explicit model^9,41 to infer the most likely initial expansion routes, with specific scenarios of population expansion that correspond to the ‘late-split’ and ‘early-split’ hypotheses proposed by linguistic studies (Supplementary Fig. 81). We ran one million Wright–Fisher simulations to test three expansion scenarios that differ in whether BSP were allowed to spread south through the Congo rainforest (that is, ‘southern route‘ or late-split; Supplementary Fig. 1b), north of the rainforest (‘northern route’ or early-split; Supplementary Fig. 1a) or both routes. In each simulation, gene genealogies were generated for selected African populations³⁰ (Supplementary Table 9). The scenario with only the northern route received substantially less statistical support from the data compared to scenarios for both routes or only the southern route (r² = 0.19, 0.32 and 0.34, respectively; Supplementary Information). Therefore, results support the late-split hypothesis, in agreement with recent linguistic, archaeological and genetic evidence^{4,15,18,19,20,42} and highlight the importance of the Congo rainforest in the initial expansion of BSP.

Previous studies proposed gene flow between the western and eastern branches of Bantu speakers^2,5. Even though populations speaking western and eastern Bantu languages are more separated in the PCA towards the terminal parts of the distribution, there is overlap toward the middle, particularly in BSP from current-day Zambia and the DRC (Fig. 3). These two countries thus represent interaction zones between different linguistic subgroupings, which is also reflected in their genetic composition (Supplementary Figs. 26 and 27). Our inference of BSP expansion routes by tracing nearest genetic distance (fixation index (F_ST)) values over the geographic landscape also indicates Zambia as a possible interaction nexus (Fig. 5a, Supplementary Figs. 82 and 83 and Supplementary Note 9). Specifically, the Lozi population represents the proxy population of the Bantu-speaking migrants from the western DRC to Zimbabwe, Mozambique, Eswatini (former Swaziland) and South Africa (Supplementary Fig. 83a,c). However, the Lozi language, widely used as a lingua franca in Zambia’s Western Province and adjacent areas, was only introduced into the region in the nineteenth century ce by Sotho-speaking immigrants from what is today South Africa^43,44. Removing the Lozi population from the analysis moves the connection point between eastern and south-eastern BSP with western BSP to the western DRC (Fig. 5b and Supplementary Fig. 83b).

Current-day Zambia as a point of divergence between expansion routes of BSP was previously proposed by ref. ³⁰ using a few BSP not representative of the whole Bantu-speaking area. Here, with notably better geographic representation of BSP, we identified Zambia and the western DRC as important nexus zones. However, spatially explicit analyses using EEMS, FEEMS, MAPS and F_ST estimates (Supplementary Information, Fig. 5c,d and Supplementary Figs. 84–91) and clustering methods (Extended Data Fig. 4) suggest barriers to gene flow and population structure even in these zones (Supplementary Note 10), possibly caused by the linguistic division between BSP. Future research and model-testing methods could further establish whether these are interaction zones between populations speaking the eastern and south-western branches of Bantu languages or splitting points in past expansion routes. Spatial methods further indicate high effective migration rates along the Indian Ocean coast from Kenya to eastern South Africa (blue areas in Fig. 5c and Supplementary Fig. 87), as was reported previously¹⁵, together with longitudinal zones of lower migration in the central parts of the continent (brown and dark red areas in Fig. 5c,d, respectively).

Dating of admixture events⁴⁵ strongly supports the main direction of BSP expansion across subequatorial Africa (Supplementary Note 11). Admixture dates (Fig. 2b) significantly correlate with geographic distance from the BSP homeland (r² = 0.20, P = 2.6 × 10⁻⁵; Fig. 2c), with earlier admixture dates in west-central Africa and more recent dates toward the extremes of the expansion (Supplementary Figs. 92 and 93 and Supplementary Table 10). These results also suggest that the rate of movement of BSP was more or less constant through time, despite the wide variety of environments and population interactions. Admixture dates seem to be older than expected in BSP from western regions (for example, between BSP and wRHG in western DRC) and younger in certain eastern regions (for example, between BSP and different eastern African groups in Uganda and Kenya) (Fig. 2b,c), suggesting that the rates of movement of BSP into these regions were either faster or slower than the average speed or that the admixture occurred earlier or later after arrival than in the other regions. Further investigation into the sociocultural aspects of the interactions with the linguistically and culturally diverse populations and the environmental challenges encountered by BSP during their expansion, particularly in adapting to diverse ecological zones and acquiring new subsistence practices, presents a promising avenue for future cross-disciplinary research.

The initial BSP expansion across subequatorial Africa may have been followed by subsequent migrations along similar routes, creating a pattern of spread-over-spread events⁴⁶. In some cases, these later migrations might have replaced earlier settlers and their languages^20,29,47. Consequently, certain branches of the Bantu language family tree may no longer accurately represent the initial BSP expansion³. This raises questions about the reliability of using only lexical and geographical data from modern Bantu languages for phylogeographic analyses to depict the ancestral BSP migration^18,19. Contact and admixture between incoming and previously settled Bantu-speaking groups could lead to genetic data reflecting a mixture of migration events, whereas linguistic data may represent only the latest spread event. As a result, both linguistics and genetics may correlate with geography but not necessarily with each other.

We tested this using Mantel tests (Supplementary Tables 11 and 12). Pairwise-population linguistic and geographic distances are significantly correlated (r-statistics = 0.6457; P = 0.0002), as are genetic and geographic distances (r-statistics = 0.1666; P = 0.0158). However, the correlation between genetic and linguistic distances is not significant (r-statistics = 0.0104; P = 0.4153). The correlation between genetics and geography increases after controlling for linguistic data, as well as between linguistics and geography after controlling for genetics (Supplementary Information). A marginally significant negative correlation between linguistic and genetic data is observed after controlling for geography (r-statistics = −0.1291; P = 0.0496). This overall weaker correlation between genetics and linguistics (whereas both correlate strongly with geography) could point to separate histories underlying the genetic and linguistic data that could involve secondary, and potentially more localized, spread waves. Other explanations are also possible, for example, admixture between linguistically distantly related BSP.

To further explore the possibility of spread-over-spread events, we compared the genetic diversity of present-day BSP and ancient (aDNA) individuals from Africa, including whole-genome sequencing data of 12 individuals from this study (97–688 years bp) and data from 83 individuals (150–8,895 years bp) from previous aDNA studies (Extended Data Fig. 1 and Supplementary Table 13). See refs. ^8,48 and Supplementary Table 14 for the archaeological and morphological descriptions and dating of the sequenced individuals (Supplementary Fig. 94 and Supplementary Information). Dimensionality reduction and clustering analyses represent genetic affinities between aDNA and modern-day individuals (Fig. 1d, Extended Data Fig. 8, Supplementary Figs. 95–98 and Supplementary Note 12). In South Africa, Late Iron Age aDNA individuals (since 688 bp) show homogeneity and genetic affinity with local modern BSP (Extended Data Table 2, Supplementary Figs. 99–101 and Supplementary Table 15), thus largely supporting a scenario of genetic continuity since the Late Iron Age. Our new Late Iron Age aDNA individuals from Zambia (since 311 BP), however, have a more heterogeneous genetic makeup showing genetic affinities with modern BSP from a wider geographical area (Supplementary Figs. 98 and 102–104). This supports the suggestion that Zambia might have been a crossroad for different movements of BSP.

Our dataset demonstrates its potential to provide an effective modern-day background genetic dataset to compare with aDNA individuals (Fig. 1d, Extended Data Fig. 8 and Supplementary Note 12). The underlying historical patterns in BSP are very difficult to distinguish on the basis of modern-day data only. Both IBD and serial-founder models can represent more complex underlying population histories among studied BSP, such as multiple overlapping expansions from the same location following similar routes. A clear manifestation of this pattern has been seen in the comparison of European history inferences based on modern DNA⁴⁹ and aDNA⁵⁰. Analyses such as our Mantel test correlations between linguistics, geography and genetics tentatively point to complex histories and possible spread-over-spread events (Supplementary Table 11), in agreement with recent archaeological studies²⁰ (Supplementary Note 1). Future aDNA studies on human remains from different archaeological contexts, associated with the Early, Middle and Late Iron Age in Africa, as well as different pottery traditions, will be necessary for assessing the affinity of the Bantu-speaker-related ancestry to each other and to current-day BSP. Therefore, the availability of our extensive genomic dataset, encompassing the full geographic expansion range of BSP, will enable further testing of these spread-over-spread proposals using aDNA.

Our study supports a large demic expansion of BSP with ancestry from western Africa spreading through the Congo rainforest to eastern and southern Africa in a serial-founder fashion. This finding is supported by patterns of decreasing genetic diversity and increasing F_ST from their point of origin, as well as admixture dates with local groups that decrease with distance from western Africa. Although our genetic findings provide less precision compared to existing linguistic models^18,19, they caution against relying solely on modern language data for tracing BSP dispersion because of potential spread-over-spread events and genetic admixture between linguistically distantly related BSP. Our genetic findings highlight the need for a comprehensive interdisciplinary study into how the demographic history of BSP influenced their language evolution. The significant correlation of admixture times with distance from the BSP source argues for a relatively constant rate of BSP expansion despite the extremely heterogeneous nature of the landscape. Although there were corridors of higher and lower effective migration rates across the African landscape, current-day Zambia and the DRC seem to be important crossroads or interaction points for the expansion of BSP. Future aDNA studies using our dataset as comparative data and new spatial modelling methods will refine our understanding of BSP expansion and their interactions with other African populations. The new findings and data will be useful not only to population geneticists, archaeologists, historical linguists, anthropologists and historians focusing on population history in Africa but also to the medical and health sector studying human genetic variation and human health in African and African-descendant populations.

Genotyping and assembled datasets

In total, 1,763 samples encompassing 163 African populations were collected with informed consent in 14 sub-Saharan African countries (Supplementary Figs. 2–4 and Supplementary Table 1). Ethical permits and sampling permission were obtained in African countries and the study as a whole was approved by the Swedish Ethical Review Board (DNR-2021-01448). DNA samples were genotyped at the SNP&SEQ Technology Platform, NGI/SciLifeLab Genomics (Sweden). We used seven genotyping batches on the Illumina Infinium H3Africa Consortium array (about 2.4 million single nucleotide polymorphisms (SNPs) and one batch on the Illumina HumanOmni2.5-Octo BeadChip. After merging all newly genotyped data and quality control steps using PLINK v.1.90b6.4 (ref. ⁵¹), genotype data consisted of 2,221,827 autosomal SNPs (Supplementary Fig. 5). After removing 67 samples because of low genotyping rate and 105 individuals because of their first- or second-degree kinship with other samples, we obtained 1,591 individuals and 2,221,827 SNPs for the ‘genotyped’ dataset. After merging the genotyped dataset with comparative data and performing quality control steps, we assembled the ‘Full-Genotyped’ dataset that contains 482,459 SNPs and 5,341 individuals from 227 populations (including 81 populations with sample sizes lower than 10 individuals) and three subdatasets with selected African populations. We included 4,950 individuals from 124 African and Eurasian populations in the AfricanNeo dataset (Supplementary Fig. 3a,b); 3,902 individuals from 111 sub-Saharan African populations in the ‘Only-African’ dataset (Supplementary Fig. 4a,c); and 2,108 individuals from 67 populations speaking Bantu languages (BSP) in the ‘Only-BSP’ dataset (Supplementary Fig. 4b,d). BSP with fewer than ten individuals were removed from specific analyses. To avoid sample-size biases, for some analyses (for example, local ancestry inference and analyses using the masked and imputed data) populations with large sample sizes were randomly downsampled to 30 individuals and we obtained 1,495 individuals from 124 populations in the downsampled admixture-masked Only-BSP dataset.

Dimensionality reduction and clustering methods

To visualize genetic variation and population structure in BSP, we applied four dimensionality reduction methods for genome-wide SNP data. We first used the uniform manifold approximation and projection (UMAP) approach⁵² directly on the genotype data. Second, we applied PCA using smartpca⁵³. To combine the information of the first ten principal components, we then used the PCA-UMAP approach⁵⁴. Fourth, we used the genotype convolutional autoencoder (GCAE) method⁵⁵. In addition, we applied an unsupervised clustering-based approach using ADMIXTURE software v.1.3.0 (ref. ⁵⁶) and cluster numbers ranging from K = 2 to K = 25.

Ancient DNA samples

To compare the genetic affinities of ancient and present-day BSP, we merged the AfricanNeo dataset with 12 aDNA individuals from southern and south-central Africa (at present Zambia and South Africa) and 83 aDNA individuals from previous studies^{23,57,58,59,60,61,62} (Extended Data Fig. 1, Supplementary Fig. 95 and Supplementary Table 13). We then projected the aDNA individuals onto a background of present-day populations using PCA. After merging haplodized modern samples and pseudohaplodized aDNA individuals and performing quality control and linkage disequilibrium-pruning steps, we used unsupervised ADMIXTURE analysis from K = 2 to K = 12.

Runs of homozygosity

We used PLINK to calculate five parameters of ROH in BSP and worldwide populations: mean ROH size, total length of ROH, sum of short ROH, sum of long ROH and ROH-based inbreeding coefficient (or F_ROH). For each studied population, we also calculated six ROH length classes. To estimate effective population sizes over the past 50 generations, we used IBDNe⁶³. To infer both the age and strength of demographic founder events in BSP, we used ASCEND v.10 (ref. ⁶⁴). To identify significant founder events, we followed the four criteria recommended by ref. ⁶⁴.

Admixture timing analysis and admixture masking

To estimate admixture dates, we applied haplotype-based admixture inference methods. We used MOSAIC v.1.4 (ref. ⁴⁵) with two- and three-way admixture models for each BSP included in the AfricanNeo dataset. For haplotype phasing of the AfricanNeo dataset, we used SHAPEIT v.2.r904 (ref. ⁶⁵). For local ancestry inference, we estimated haplotypic admixture from six reference panels in individuals from BSP using RFMix software v.1.5.4 (ref. ³⁷). To avoid the influence of admixture patterns in BSP in our ancestry-specific analyses, we removed haplotypes with non-WCA-related ancestry from each haploid genome of each Bantu-speaking individual using a masking approach (Supplementary Figs. 105–107). For each assembled dataset, we explored patterns of population structure between and within populations using smartPCA⁵³.

Phylogenetic analyses and correlations

To investigate phylogenetic relationships between all the BSP, we used TreeMix v.1.13 (ref. ⁶⁶). The likelihood of each proposed population-based maximum-likelihood TreeMix topology was assessed by bootstrapping blocks of 500 SNPs and assigning Ju/’hoansi from Namibia as the root of the population tree. To test the correlation between genetic, linguistic and geographical distances, we performed Mantel tests and partial Mantel tests using the R package ncf⁶⁷. As genetic distances, we computed pairwise F_ST between populations included in the ancestry-masked Only-BSP dataset using EIGENSOFT package v.6 (ref. ⁶⁸). Geographic distances were calculated as pairwise great circle distances between the studied populations using the R package geosphere⁶⁹. For linguistic distances, we used a linguistic dataset from the multistate matrix of cognate sets identified by ref. ¹⁸. In total, 38 BSP matched the genetic dataset and the linguistic dataset of 409 Bantu languages studied by ref. ¹⁸.

Patterns of genetic diversity

To investigate spatial patterns of genetic diversity of studied African populations, we calculated statistics based on haplotype diversity and linkage disequilibrium information. Haplotype heterozygosity and haplotype richness were computed following recommendations from ref. ¹² with homemade scripts implemented in Python. Each calculation was also repeated ten times. Values were calculated per chromosome and then averaged across the genome of each individual in each population. We characterized linkage disequilibrium patterns in each population with more than ten individuals by measuring the correlation coefficient (r²) between all pairs of SNPs within 500 kilobase pairs windows using PLINK. To assess whether haplotype heterozygosity, haplotype richness and linkage disequilibrium patterns in BSP were consistent with a history of expansion from the homeland of BSP, we performed a linear regression between the three summary statistics and geographical distances from Cameroon, assuming that the BSP expansion started in that region³.

Pairwise genetic distances

To reconstruct potential routes of expansion of BSP, pairwise F_ST values were calculated between one population from Cameroon (Nzime) and each of the studied BSP for the masked and imputed Only-BSP dataset. We also applied the GenGrad method from ref. ⁷⁰ but using F_ST as the genetic distance metric and with slightly adjusted parameters to better fit the smaller study area.

Effective migration rates

To further investigate spatial population structure in sub-Saharan African populations, we used EEMS software⁷¹ and its implementation FEEMS⁷² and MAPS⁷³. EEMS and FEEMS were performed on the Only-African and Only-BSP datasets before and after using the masking approach and MAPS was performed on the masked Only-BSP dataset. EEMS analysis was repeated three times and an average was taken as input for the visualization as recommended in the EEMS manual.

Testing isolation-by-distance models

To test four models of migration, we used SpaceMix v.0.13 (ref. ⁷⁴). The software generates geogenetic maps in which genetics rather than physical distances determine the distances between individuals/populations. The general underlying assumption evaluated with SpaceMix is that under an IBD pattern, geographic and geogenetic positions will be similar, which is a pattern of IBD. The best-fitting model was evaluated using Pearson correlations between the expected and observed data.

Testing models of migration routes

To test different demographic scenarios for the BSP expansion, we used a spatiotemporally explicit population genetic framework^9,41. Here, we adapted the extension of the model presented by ref. ⁴¹ to apply multiple local expansions for different scenarios of expansion. We considered three demographic scenarios in which the expansion of BSP proceeded north of the rainforest, south through a rainforest corridor or using both northern and southern routes. For each demographic scenario, we ran one million simulations with parameters drawn from an independent uniform distribution for parameters characterizing the BSP expansion and from a log-uniform distribution for parameters describing the initial global expansion of anatomically modern humans taken from ref. ⁴¹.

Ethics and inclusion

This study was conducted according to the Declaration of Helsinki (World Medical Association, 2013). DNA samples were collected with informed consent from participants. Ethical permits and sampling permission were obtained in African countries (Methods) and the study as a whole was approved by the Swedish Ethical Review Board (DNR-2021-01448). The sampling for this study emerged from population genetic and archaeology projects that involved local research institutions and the participation of local communities. Local institution involvement included research design, selection of archaeological material for analyses, modern-day DNA sample collection, community involvement, permit clearance, feedback on analyses, help with writing and feedback on the manuscript. Before submission of this study for publication, the corresponding author and first author participated in-person in the ICHG conference held in Cape Town, as well as held online presentations and meetings with local researchers.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.