African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana

To examine the relationship between African individuals and other worldwide samples, we used three complementary clustering approaches. Distance-based clustering by neighbor-joining reveals a clear split between Eurasian and African samples, indicating deep divergence between the continents (Fig. 2A and SI Appendix, Fig. S1). The majority of Eurasian samples form a nearly star-shaped phylogeny, consistent with recent expansion of these lineages. Conversely, longer more bifurcated branches separate African subclusters and previously identified Eurasian relicts from each other and from the nonrelict clade. In general, the Eurasian clades cluster consistently with the nine groups defined previously (4). Exceptions are the Central Europe clade, which separates into two clusters, and the Iberian relicts, which cluster with the Moroccan Rif-Zin population. Moroccan samples separate into four clades, reflecting their geographic distribution and South Africa and Tanzania cluster together in a single clade. The results for South Africa and Tanzania are striking because A. thaliana populations outside of Eurasia and North Africa were previously thought to be recently introduced by humans (27).

Similarly, principal component analysis (PCA) separates African populations from each other and from Eurasian populations (Fig. 1B). The first PC distinguishes sub-Saharan Africa, and the second separates the four Moroccan clusters from Eurasians. Subsequent PCs mainly discriminate populations within Africa, whereas Eurasian populations remain tightly clustered (SI Appendix, Fig. S3). Results from ADMIXTURE (28) reinforce this finding. Moroccan populations separate into three clusters and are distinguished from a single cluster of Eurasian samples (Fig. 1C, and SI Appendix, Figs. S4 and S5, and Table S2). PCA and ADMIXTURE results also suggest a Moroccan origin of the relicts in Iberia, which are spread between Rif-Zin and Iberian nonrelicts in PCA, and sizable portions of the Iberian relict genomes match the Moroccan clusters in ADMIXTURE. This finding is consistent with previous work (29) and with the accepted phylogeographical history of Mediterranean and North African flora characterized by a complex history of expansions and contractions driven by important climatic changes experienced in this vast region, particularly since the Pliocene (30, 31).

Furthermore, from pairwise differences, we recover the previously reported difference between Eurasian relicts and nonrelicts (4) and find that all African accessions are at least as divergent as samples previously classified as relicts (Fig. 3A and SI Appendix, Fig. S6). Therefore, in contrast to Eurasia, where most samples represent a single recently spread clade, all African individuals represent relict samples. The distribution of pairwise differences within Africa (Fig. 3B) further demonstrates the high diversity in these samples.

If populations in Africa truly are more ancient than the Eurasian clusters, we should also expect higher numbers of private variants in Africa. Indeed, we find that the Moroccan clusters and Iberian relicts, which appear likely derived from Morocco, harbor the highest numbers of private SNPs (Fig. 3C and SI Appendix, Fig. S8 and Table S3). This signal intensifies when we exclude recently arisen variants, which we find constitute the majority of the private variation in Eurasia. First, we excluded singletons, the class of SNPs most influenced by recent population growth, and found a 7.0- to 23.2-fold enrichment in Morocco compared with the top nonrelict cluster (Fig. 3C). Next, we considered the spatial distribution of private variants across the genome. Because novel private variants are unlikely to be tightly linked, clustering of several, contiguous private SNPs indicates old haplotypes. At the haplotype level, we again found extremely high enrichment (3.9- to 20.9-fold) for Moroccan clusters and Iberian relicts (4.5-fold) relative to the top Eurasian cluster (Fig. 3C). Notably, the Eastern Mediterranean and Caucasus regions, which were previously favored as points of origin and major centers of diversity of the species (3), do not exhibit a striking pattern for any of the metrics examined (see IBC in Fig. 3C). These findings evoke a model in which polymorphism in Africa is because of ancient variation and Eurasian polymorphism is mainly because of recent expansion.

To better understand the relevance of variation in Africa for the early history of the A. thaliana lineage, we examined variation at the locus that confers self-compatibility. S-locus variation in Africa differs in several ways from what is found in Eurasia (Fig. 4 and SI Appendix, Table S4). We found all three S-haplogroups together in a single geographic region, with S-haplogroup B private to Africa. In addition, the A-C recombinant is also present at low frequency in Morocco, in contrast to what is found in Eurasia (32). Finally, we discovered deletion haplotypes in haplogroups A and C in Morocco (SI Appendix, Fig. S9). The finding that all S-locus haplogroups are present together implies that selfing evolved in a single geographic region.

Taken together, the patterns in population structure, and levels of variation across the genome and at the S-locus, specifically, imply a deep history in Africa. To clarify the details of past demographic events, we inferred historical effective population sizes (N_e) and split times among populations based on cross-coalescent rates (CCR) using a multiple sequentially Markovian coalescent approach (MSMC) (33).

In ancient times, we find the highest N_e is in Africa, peaking at around 500–400 kya (Figs. 5 and SI Appendix, Fig. S10). All Eurasian populations (including Iberian relicts) show the same trajectories as the Africans, but with lower amplitudes. Given that these curves are in phase with one another and we do not see evidence for a population split until 120 kya (Fig. 6A), we interpret this as population structure in the ancestral population combined with bottlenecks as more derived populations migrated away from this ancestral population. This finding is consistent with our finding that variation in Eurasia is often a subset of variation present within Africa and is similar to the situation in humans (34). Notably, the IBC cluster, which includes previously hypothesized A. thaliana origins and refugia (3, 5, 7), exhibits a much lower ancient population size than Africa (SI Appendix, Fig. S10).

At 120–90 kya, there are bottlenecks in all populations (Fig. 5) and a split among the major clades (Moroccan, Tanzanian, and Levant) (Fig. 6 and SI Appendix, Fig. S11). This roughly corresponds to the Abbassia Pluvial, a period when migration corridors were open because of high precipitation and humidity in Africa (35) and also marks the last interglacial at Marine Isotope Stage 5e (130–116 kya), when temperatures were 1–2° warmer than present-day conditions, providing favorable conditions in Eurasia. As this interglacial period came to a close, there was a worldwide shift toward cooler, drier conditions as the most recent and severe Pleistocene glaciation phase began (36). Beginning at this time, CCR implies a progressive decline in population connectivity, consistent with decreasing temperature and increasing aridity (Fig. 6A). We checked for consistency using a complementary method that relies on the joint site frequency spectrum between populations (δaδi) (37) and found a slightly older estimate for the split and overlapping confidence intervals (141–116 kya) (SI Appendix, Table S5). We propose that the most likely scenario is that A. thaliana was colonizing broadly within Africa as well as in the Levant during the last interglacial (130–116 kya), but that connections between populations began to break down as populations spread and as climate became cooler and drier at the end of this period.

MSMC based on eight haplotypes detects several changes in N_e in more recent times (Fig. 5B). Since ∼120 kya, the population size changes in Europe and Asia are often out of phase with those in Africa, consistent with geographical separation and exposure to different climatic regimes. Maxima in African populations occur at around 60–40 kya and 11–5 kya, corresponding to orbital-scale climate shifts from arid to moist conditions (i.e., pluvial periods) that occurred at 11–5 kya and 59–47 kya in North Africa (38, 39). We find relatively recent split times between sub-Saharan African clades (South Africa and Tanzania) and between Western Europe and Central Asia, at around 40 kya (Fig. 6).

There are a few caveats to consider regarding the demographic inference. Because MSMC can spread instantaneous population size changes over time, maxima and minima are informative but the slope should be interpreted with caution (33). In addition, the precise timing associated with inferred population size changes and splits is dependent on parameters that are difficult to measure and may vary over space and time, including mutation rate, degree of purifying selection, and the possible input from a seed bank. We used the best available data for mutation rate [based on mutation accumulation experiments (40)] and made the usual simplifying assumptions for other parameters (one generation per year), but the timing we infer would need to be revised if these assumptions turned out to be incorrect.

Genomic studies thus far have amassed data for nearly 2,000 Eurasian A. thaliana accessions but were unable to provide insight into the early history of the species. Here, in a genome-scale sequencing effort focused on African accessions, we find clear evidence for a deep history of African A. thaliana populations, which harbor variation that was either lost or never present in Eurasia. Several specific results were unexpected based on current knowledge in this well-studied species. First, we discovered surprising and clear evidence that A. thaliana is native not only to North Africa but also to Afro-alpine regions of sub-Saharan Africa. Second, our results revealed that the deepest splits species-wide separate the African lineages from one another and that in ancient times, the effective population size was largest in Africa. Finally, we learned that variation at the S-locus is highest in Africa and that all three S-haplogroups are present there.

Based on our results, we can outline a model for the early history and transition to selfing in A. thaliana (detailed in SI Appendix, Fig. S12). In the first step, we infer that the population ancestral to A. thaliana became geographically separated from its parental out-crossing population. Our results suggest that this separation involved migration of the ancestral subpopulation into Africa by 1.2–0.8 Mya. This timing corresponds to the Middle Pleistocene Transition, a shift to drier more variable climate and more open habitats in Africa (i.e., grasslands versus woodlands), as evidenced by soil carbon analysis showing an increase in the ratio of C4 to C3 plants (41, 42).

Although the estimated divergence times between A. thaliana and Arabidopsis lyrata center around 5–7 Mya (9, 43), the origin of A. thaliana itself appears to be much younger. Our model predicts that there was an initial bottleneck as the subpopulation that led to A. thaliana split from a A. lyrata-like ancestral population [similar to that observed in Mimulus nasutus (44) and Capsella rubella (45–47)], followed by an expansion in N_e as the selfing population began to spread. In this case, we could interpret the MSMC results to suggest that the transition to selfing occurred between 1 Mya and 500 kya, before the most ancient maximum in N_e. This finding is in line with an estimate based on the depth of the A. thaliana genealogy (0.84% maximum divergence among individuals sampled here) under a simple model (T ∼ D/2 μ ∼ 598 kya). Our estimated timing is also consistent with previous estimates for the loss of self-incompatibility and origin of selfing (24–26).

Once selfing was established, traits associated with the “selfing syndrome” would have been favored, including reduced pollen number and petal size (48). Such phenotypic shifts are common in predominantly selfing species and have occurred in A. thaliana compared with its closest relatives (26). At the genomic level, A. thaliana exhibits major chromosomal rearrangements and a reduction in genome size and number of chromosomes (49). This genomic reduction is also likely a by-product of the shift to predominant selfing in A. thaliana (9–11), consistent with an observed link between reduced genome size and selfing in other plant species (11–13). These changes introduce a strong reproductive barrier as found in hybrids of A. thaliana and A. lyrata, which are infertile because of the chromosomal rearrangements that occurred in A. thaliana (14).

Given that all three S-haplogroups co-occur in Morocco, we hypothesize that the transition to predominant selfing occurred in a single region, best represented today in Morocco. This finding differs from previous assertions that these events likely happened separately in geographically distinct populations (15, 16). Moreover, it allows for the possibility that the transition to selfing was aided by a shared precursor mutation, a shared climate, and the bottleneck that occurred during the migration away from the ancestral population (22). Our proposed model parallels observations in partially selfing populations of A. lyrata (50). Here, self-compatibility is associated with two different S-haplogroups in Great Lakes populations and self-compatibility may have been favored because of the bottleneck that initially limited S-haplogroup diversity and thus mate availability.

After the origin and initial population size increase of A. thaliana, we infer several demographic changes that are congruent with known climatic shifts. At 120–90 kya, we find evidence from MSMC and δaδi for splitting among the major clades: Morocco, Levant and sub-Saharan Africa. This split corresponds to the Abbassia pluvial, which produced migration corridors within Africa (120–90 kya) (35, 39) as well as Marine Isotope Stage 5e (130–116 kya), the last interglacial period, when worldwide temperatures were 1–2° warmer than they are currently (51, 52). This is consistent with a model in which A. thaliana spread widely throughout Africa and into Eurasia when conditions were favorable (∼120 kya), with isolation as gene flow was reduced (SI Appendix, Fig. S12). More recent major demographic events include the split between European and Asian populations at around 40 kya and the increase in N_e within Africa during the most recent pluvial.

The patterns we observe and their concordance with climatic events suggest that the transition to selfing and speciation occurred within Africa, with subsequent migration out of Africa into Eurasia. However, it is also possible that the initial transition to selfing occurred within Eurasia followed by migration into Africa and concomitant loss of variation in Eurasia. This alternative would require that the ancient variation in the A. thaliana lineage was either lost or has not been sampled in Eurasia and the bottleneck into Africa was mild enough to preserve high levels of genetic variation.

Overall, the patterns in A. thaliana bear striking similarities to those observed for human populations, particularly in the larger effective population size in Africa (34), the exodus from Africa approximately 120 kya (39, 53–55), and the splitting of major human populations in Europe and Asia (approximately 45–35 kya) (53, 54). Analogous to what we propose here, demographic events in human populations have been attributed to major climate transitions (35, 39, 56).

Moreover, the timing and types of demographic events we infer during the history of A. thaliana are consistent with previous observations in a broad range of other plant species. Specifically, the shift to predominance of C4 plants across Africa at 1.2–0.8 Mya and the intensification of glacial cycles worldwide (57) correspond with our estimated timing of the evolution of selfing in A. thaliana and a clustering of speciation events more generally (58). The geographic expansion approximately 120 kya corresponds to an African pluvial and worldwide interglacial, which resulted in expansion of forests across Africa (59) and Eurasia (51, 52). Finally, we see evidence of an increase in effective population size overlapping with the most recent and well-described African pluvial at 11–5 kya, when the Sahara was heavily vegetated and filled with lakes (38, 60, 61). The concordance between inferred population size changes, climate, and reports for other species implies that the patterns we observe in A. thaliana may be representative of climate-mediated population dynamics across diverse taxa.

For full materials and methods, please see SI Appendix, Supplementary Text.

We sequenced the genomes of 79 A. thaliana individuals, including 70 fresh samples and 9 herbarium samples (SI Appendix, Table S1). For fresh leaf samples, sequencing libraries were prepared using Illumina TruSeq DNA sample prep kits (Illumina) and sequenced on Illumina Hi-Seq instruments. DNA from herbarium specimens was extracted, authenticated, and treated with uracil glycosylase to remove damaged nucleotides in a clean room facility at the University of Tübingen. To align the sequences to the TAIR10 reference genome and to call variants, we used two different pipelines: the MPI-SHORE pipeline (62) and a more conservative pipeline designed to reduce false positives resulting from indels.

For population structure analyses, we subsampled the complete dataset to match sample sizes across clusters as some Eurasian geographic regions are heavily oversampled, which could cause biases in some analyses, and we pruned SNPs based on LD to select a representative set. For ADMIXTURE, the number of clusters (K) was determined based on the outcome of cross-validation analyses.

To infer patterns of effective population size and population separations over historical time, we used a MSMC v2 (33). Because A. thaliana accessions are inbred, we created pseudodiploids by combining chromosomes from pairs of individuals from the same populations and ran MSMC in the two- and eight-haplotype configurations (Fig. 5). We assumed a mutation rate of 7.1 × 10⁻⁹ based on results of mutation accumulation experiments (40) and a generation time of 1 y. We confirmed inferences using δa δi (37) on joint site frequency spectra from pairs of populations.

We thank the South African National Biodiversity Institute, the University of Vienna, and Naturalis for access to herbarium samples; and V. Castric and M. Horton for helpful comments on the manuscript. Support for sequencing and plant growth were provided by the Vienna Biocenter Core Campus Science Support Facilities and the Max Planck Gesellschaft (MPG). Computational resources were from the Vienna Scientific Cluster, Center for Integrative Bioinformatics Vienna, and MPG were used. The project was supported by European Union CIG 304301, European Research Council Grant CVI_ADAPT 638810, and MPG (to A.M.H.); Austrian Science Funds FWF W1225; Grants BIO2016-75754-P (to C.A.-B.) and CGL2016-77720-P (to F.X.P.) from the Agencia Estatal de Investigación of Spain and the Fondo Europeo de Desarrollo Regional; Japan Society for the Promotion of Science KAKENHI 15K18583 (to T.T.); and the Presidential Innovation Fund of the MPG (to H.A.B.).