The interplay of demography and selection during maize domestication and expansion

Abstract

Background

The history of maize has been characterized by major demographic events, including population size changes associated with domestication and range expansion, and gene flow with wild relatives. The interplay between demographic history and selection has shaped diversity across maize populations and genomes.

Results

We investigate these processes using high-depth resequencing data from 31 maize landraces spanning the pre-Columbian distribution of maize, and four wild teosinte individuals (Zea mays ssp. parviglumis). Genome-wide demographic analyses reveal that maize experienced pronounced declines in effective population size due to both a protracted domestication bottleneck and serial founder effects during post-domestication spread, while parviglumis in the Balsas River Valley experienced population growth. The domestication bottleneck and subsequent spread led to an increase in deleterious alleles in the domesticate compared to the wild progenitor. This cost is particularly pronounced in Andean maize, which has experienced a more dramatic founder event compared to other maize populations. Additionally, we detect introgression from the wild teosinte Zea mays ssp. mexicana into maize in the highlands of Mexico, Guatemala, and the southwestern USA, which reduces the prevalence of deleterious alleles likely due to the higher long-term effective population size of teosinte.

Conclusions

These findings underscore the strong interaction between historical demography and the efficiency of selection and illustrate how domesticated species are particularly useful for understanding these processes. The landscape of deleterious alleles and therefore evolutionary potential is clearly influenced by recent demography, a factor that could bear importantly on many species that have experienced recent demographic shifts.

Background

Genomes are shaped over the course of their evolutionary history through a complex interaction of demography and selection. Neutral processes that comprise a species’ demographic history, such as stochastic changes in population size and migration events, influence both the pool of diversity upon which selection can act and its efficiency. Selection and genetic drift then jointly determine the fate of this diversity.

After the development of agriculture, both crops and humans have experienced profound demographic shifts that left clear signatures in genome-wide patterns of diversity [1, 2]. Early agriculturalists sampled a subset of the diversity present in crop wild relatives, resulting in an initial demographic bottleneck for many domesticates [3]. Subsequent to domestication, humans and their crops experienced a process of global expansion facilitated by the rise of agriculture [4]. In many cases expansion was accompanied by gene flow with close relatives, a demographic process that further altered patterns of diversity [5, 6].

Recent interest in the effects of demography on functional variation has led to a growing body of theory that is increasingly supported by empirical examples. To date, the relationship between demography and selection has been most thoroughly explored in the context of deleterious alleles. While theory suggests mutation load (i.e., the reduction in mean fitness caused by the presence of deleterious alleles) may be insensitive to demography over long periods [7, 8], empirical results are consistent with load being shaped by demography over shorter timescales [9–13]. For example, evidence in both plant and animal species has revealed increased mutation load in populations that have undergone recent, sudden declines in effective population size (N _e ) [10–12, 14]. Similarly, in geographically expanding populations, repeated sub-sampling of diversity (i.e., serial founder effects) can occur during migration away from a center of origin [15, 16], a phenomenon shown to have decreased genetic diversity and increased counts of deleterious alleles in human populations more distant from Africa [17, 18]. Finally, gene flow may also affect genome-wide patterns of deleterious variants, particularly when occurring between populations with starkly contrasting N _e . For instance, during the Out-of-Africa migration, modern humans inter-mated with the Neanderthal species, a close relative with substantially lower N _e and higher mutation load [9]. The higher mutation load in Neanderthals presented a cost of gene flow, and subsequent purifying selection appears to have limited the amount of Neanderthal introgression near genes in the modern human genome [9, 19].

The domesticated plant maize (Zea mays ssp. mays) has a history of profound demographic shifts accompanied by selection for agronomic performance and adaptation to novel environments, making it an ideal system in which to study the interaction between demography and selection. Maize was domesticated in a narrow region of southwest Mexico from the wild plant teosinte (Zea mays ssp. parviglumis; hereafter, parviglumis [20–22]) and experienced an associated genetic bottleneck that removed a substantial proportion of the diversity found in its progenitor [23, 24]. Archaeological evidence suggests that after initial domestication, maize spread across the Americas, reaching the southwestern USA by approximately 4500 years before the present (BP) [25] and coastal South America as early as 6700 years BP [26]. Gene flow into maize from multiple teosinte species has been documented in geographical regions outside of its center of origin [5, 27]. To date, genetic studies of demography and selection in maize have primarily focused on initial domestication [28], only broadly considering the effects of subsequent change in population size on diversity [2] and largely disregarding the spatial effects of geographic expansion and gene flow (but see [29]). Furthermore, the effect of maize demography on the prevalence of deleterious alleles has yet to receive in-depth attention.

Here, we investigate the genome-wide effects of demographic change in maize during domestication and subsequent expansion using high-depth resequencing data from a panel of maize landraces. We present evidence for a protracted domestication bottleneck, further loss of diversity during crop expansion, and gene flow between maize and its wild relatives outside of its center of origin. We then explore how this demographic history has shaped genome-wide patterns of deleterious alleles.

Results

Maize population size change during domestication and expansion

We resequenced 31 maize individuals, each from one open-pollinated landrace, representing six geographical regions that span the pre-Columbian range of maize cultivation (southwestern US highlands, 6 individuals; Central Mexican Plateau, 6 individuals; Mexican lowlands, 5 individuals; Guatemalan highlands, 3 individuals; South American lowlands, 6 individuals; Andes, 5 individuals). In addition, we resequenced four wild parviglumis individuals from a single population located in the Balsas River Valley in Mexico (Fig. 1a). The median sequencing depth was 29X, with a range of 24–53X, resulting in a data set consisting of 49,508,640 single nucleotide polymorphisms (SNPs). Landrace accessions were selected to broadly reflect the diversity of maize in the Americas and to be representative of defined ecogeographic regions based on consultation with experts on landrace germplasm (Major Goodman, personal communication) and on descriptions in the Races of Maize handbooks [30].

Fig. 1

Maize domestication and expansion. a Sampling locations. b Estimates of effective population size over time (mutation rate =3∗10⁻⁸, generation time = 1 year). Dashed lines represent bootstrapping results. The x axis is l o g10 scaled when time is less than 10,000 generations BP and linear when greater than 10,000 generations BP as indicated by the gray background. c The percentage of polymorphic sites versus distance from the maize domestication center. Abbreviations for populations: GuaHigh Guatemalan highlands, MexHigh Mexican highlands, MexLow Mexican lowlands, SA_Low South American lowlands, SW_US southwestern US highlands

We first estimated historical changes in effective population size (N _e ) of maize and parviglumis using the multiple sequentially Markovian coalescent (MSMC) [31]. Consistent with archaeological evidence [21], we find that the demographic histories of the various maize populations begin to diverge from one another approximately 10,000 years BP (Fig. 1b). Surprisingly, our single population of parviglumis diverges from maize much earlier, around 75,000 years BP. All maize populations show a gradual decline in diversity concomitant with divergence from parviglumis, but the slope becomes more pronounced around the time of domestication. This period of declining N _e continues until the recent past (≈ 1100−2400 years BP) and is followed by extremely rapid population growth, suggesting recovery from domestication post-dated expansion of maize across the Americas. In contrast to our results in maize, parviglumis shows an increase in N _e which also lasts until the recent past (≈ 1200−1800 years BP). To determine if linked selection associated with domestication could bias estimates of N _e in maize (see [32]), we masked previously identified domestication candidates [24] and observed nearly identical results (Additional file 1: Figure S1A).

One explanation for the prolonged population size reduction in maize following the onset of domestication would be repeated colonization bottlenecks during the spread of maize across the Americas. Genome-wide levels of heterozygosity across our maize samples are consistent with this idea, showing a strong negative correlation (R ²=0.3636,p=0.0004; Fig. 1c) with distance from the center of maize domestication in the Balsas River Basin. To confirm this trend, we performed a similar analysis with a much larger sample of published genotyping data (n=3520; Additional file 1: Figure S1B) [33] and observed similar results.

While the gradual decrease in genetic diversity seen with distance from the Balsas indicates serial founder effects, our analyses also point to a more extreme founder event in the Andean highlands of South America. Andean landraces show a deeper bottleneck in our MSMC analysis (Fig. 1b), have the lowest overall diversity (Additional file 1: Figure S2), and show both a distinct reduction of low frequency alleles and a greater proportion of derived homozygous alleles compared to other populations (Additional file 1: Figure S2). To shed light on the timing of this extreme founder event, we assessed evidence for recent inbreeding. Inbreeding coefficients in Andean samples were quite low and not statistically different from other populations (all F<0.002 and p>0.05 based on a Wilcoxon test). Likewise, no significant difference could be found across populations in the number of runs of homozygosity (ROHs) longer than 1 cM (p>0.05 in all cases, Wilcoxon test). Using simple conversions between generations and the genetic length of an inherited region in the genome [34], these results provide further evidence for limited recent (< 50 generations) inbreeding in the Andes. However, when ROHs were limited to those shorter than 0.05c M and longer than 0.005c M (inbreeding from approximately 1000–10,000 generations in the past), Andean samples demonstrated significantly greater cumulative ROHs compared to all (p<0.05, Wilcoxon test) but the South American lowland population (p=0.165, Wilcoxon test; Additional file 1: Figure S3). Together, these lines of evidence are consistent with an unusually strong founder event during colonization of the Andes.

Introgression from wild maize in highland populations

Adaptive introgression from the wild teosinte taxon Zea mays ssp. mexicana (hereafter, mexicana) has previously been observed in maize in the highlands of Mexico [5]. Our broad sampling allowed us to investigate whether introgressed mexicana haplotypes have spread to highland maize populations outside of Mexico, potentially playing a role in adaptation in other regions. In order to test this hypothesis, we calculated Patterson’s D statistic [35] across all maize populations. All individuals from both the Mexican and Guatemalan highlands exhibited highly significant evidence for shared ancestry with mexicana (Additional file 1: Figure S4). Maize from the southwestern USA also showed more limited evidence of introgression, consistent with findings from ancient DNA suggesting this region was originally colonized by admixed maize from the highlands of Mexico [36]. In contrast, the distribution of z-scores for South American populations overlapped zero, providing no evidence for substantial spread of mexicana haplotypes to this region.

We localized introgression to chromosomal regions through genome-wide calculation of the $\hat {f_{d}}$ statistic [37]. Megabase-scale regions of introgression were identified in both Mexican and Guatemalan highland populations that correspond to those reported by [5] on chromosomes 4 and 6 (Fig. 2; Additional file 1: Figure S5). On chromosome 3 (at around 75−90 Mb), a large, previously unidentified region of introgression can be found in the Mexican and southwestern US highlands (Fig. 2; Additional file 1: Figure S5). This region overlaps a putative chromosomal inversion associated with flowering time in maize landraces [38] and in the maize nested association mapping population [39] and may be an example of mexicana contribution to modern maize lines.

Fig. 2

Introgression from mexicana into maize landraces. Loess regression of $\hat {f_{d}}$ is plotted for all five populations on a chromosome 3 and b chromosome 4. Each plot highlights a single population, with other populations shown in gray. The Mexican lowlands population is used as a reference and thus not plotted. No significant introgression was detected in the South American lowlands or the Andes, and loess regressions for these populations are only shown as gray lines. The statistic $\hat {f_{d}}$ was calculated based on the tree in which P2 is varied across populations. mex mexicana, Trip Tripsacum

The influence of demography on accumulation of deleterious alleles

Population-specific changes in historical N _e should influence the efficiency of purifying selection and alter genome-wide patterns of deleterious variants [10]. Introgression from a species with substantially different N _e may also influence the abundance and distribution of deleterious alleles in the genome [9, 19]. Below we evaluate the effects of major demographic events during the pre-Columbian history of maize on patterns of deleterious alleles.

Domestication and deleterious alleles

We first compared counts of deleterious alleles in Mexican lowland maize individuals to four parviglumis individuals from a single population in the Balsas River Valley. Maize from the Mexican lowlands has not experienced substantial introgression from wild relatives and is near the center of maize origin [22], and thus best reflects the effects of domestication alone. After identifying putatively deleterious mutations using Genomic Evolutionary Rate Profiling (GERP) [40], we calculated the number of derived deleterious alleles per genome under both an additive and a recessive model across four levels of mutation severity (see Methods for details). Maize showed significantly more deleterious alleles than teosinte under both additive (<1 0% more; p=0.0079, Wilcoxon test; Additional file 1: Figure S6) and recessive (< 20−30% more; p=0.0079; Fig. 3) models across all categories (Additional file 1: Figure S7). Additionally, maize contained more than twice as many fixed deleterious alleles than teosinte (57,881 versus 26,947) and 10% fewer segregating deleterious alleles (429,837 versus 478,594), effects expected under a domestication bottleneck (Fig. 3c; [7]). GERP load (GERP score × frequency of deleterious alleles), a more direct proxy of mutation load quantified at the population level, revealed a similar trend (additive model: maize median =23.635, teosinte median =22.791, p=0.008, Wilcoxon test; recessive model: maize median =14.922, teosinte median =12.231, p=0.008). Maize, like other domesticates [12, 14, 41, 42], thus appears to have a higher mutation load compared to its wild progenitor parviglumis.

Fig. 3

Burden of deleterious mutations during maize domestication and expansion. Comparison of counts of deleterious alleles at the individual level a between parviglumis and maize (mean value in parviglumis population was used as the standard to calculate the relative burden) and b among maize populations (mean value in Mexican lowland population was utilized as the standard to calculate the relative burden) under a recessive model. Comparison of fixed versus segregating (seg) deleterious alleles at the population level c between parviglumis and maize and d among maize populations. A jackknife sub-sampling approach (n = 4) was utilized for maize in c and for individual maize populations (n = 3) in d

While the elevated mutation load we observe in maize relative to parviglumis may be driven primarily by the domestication bottleneck, positive selection on causal variants underlying domestication phenotypes may also fix nearby deleterious variants through genetic hitchhiking, which would result in a higher number of deleterious alleles in regions linked to domestication loci [41, 43]. To test this hypothesis, we first confirmed that 420 previously identified domestication candidates [24] showed evidence of selection in our data (Additional file 1: Figure S8), and then assessed the distribution of deleterious alleles in and near (5 kb upstream and downstream) these genes by calculating the number of deleterious alleles per base pair under both recessive and additive models. No significant difference was found in the prevalence of deleterious alleles near domestication and random sets of genes (Additional file 1: Figure S9), suggesting the increased mutation load we observe in maize has been driven primarily by the genome-wide effects of the domestication bottleneck rather than linkage associated with selection on specific genes.

The effect of the Andean founder event on deleterious alleles

The extreme founder event observed in the Andes could potentially alter genome-wide patterns of deleterious variants beyond the effects of domestication alone. Under a recessive model, maize from the Andes contains significantly more deleterious alleles than any other population (Fig. 3b; Additional file 1: Figure S7; all p values <0.02, Wilcoxon test), and this difference becomes more extreme when considering more severe (i.e., higher GERP score) mutations (Additional file 1: Figure S7). In contrast, we observe no significant difference under an additive model (Additional file 1: Figure S6; Additional file 1: Figure S7). The Andean founder event therefore appears to have resulted in higher mutation load than seen in other maize populations. This result is further supported by a higher proportion of fixed deleterious alleles within the Andes and fewer segregating deleterious alleles (Additional file 1: Figure S10; Fig. 3d), a result comparable to the differences observed between maize and parviglumis.

Introgression decreases the prevalence of deleterious alleles

Highly variable rates of mexicana introgression were detected across our landrace populations (Fig. 2; Additional file 1: Figure S4; Additional file 1: Figure S5). To explore the potential effects of introgression on the genomic distribution of deleterious alleles, we fit a linear model in which the number of deleterious sites is predicted by introgression (represented by $\hat {f_{d}}$) and gene density (exonic base pairs per centimorgan) in 10-kb non-overlapping windows in the Mexican highland population where we found the strongest evidence for mexicana introgression. Gene density was included as a predictor in the regression to control for the positive correlation observed between gene density and both introgression (p=3.48e−08) and deleterious alleles (p≈0). When identifying deleterious alleles under both additive and recessive models, we found a strong negative correlation with introgression (i.e., fewer deleterious alleles in introgressed regions; p≈0 under both models). These findings likely reflect the larger ancestral N _e and more efficient purifying selection in mexicana.

Discussion

Demographic studies in domesticated species have focused largely on identifying progenitor population(s) and quantifying the effect of the domestication bottleneck on genetic diversity [24, 44, 45]. It is likely, however, that the demographic history of domesticates is generally more complex than a simple bottleneck followed by recovery [46, 47]. Many crops and domesticated animals have expanded from defined centers of origin to global distributions, experiencing population size changes and gene flow from closely related taxa throughout their histories [48]. With this in mind, we have characterized maize demography from domestication through initial expansion in order to provide a more complete assessment of the influence of demography on deleterious variants.

Historical changes in maize population size

Early models of maize demography suggested the ratio of the domestication bottleneck size and duration was between ≈ 2.5:1 and ≈ 5:1, but little statistical support was found for specific estimates of these individual parameters [23, 28, 49]. Most recently, Beissinger et al. [2] fit a model assuming a bottleneck followed by instantaneous exponential recovery. While our results concur with the most recent model in finding a similar bottleneck size (≈ 10% compared to ≈ 5% in Beissinger et al.) and that the modern N _e of maize is quite large, the flexibility of MSMC also allowed us to estimate the duration of the bottleneck. We find that the domestication bottleneck may have lasted much longer than previously believed, spanning ≈ 9000 generations and only beginning to recover in the recent past (Fig. 1b). Analysis of bottlenecks during African rice and grape domestication have also suggested a duration of several thousand generations [46, 47], indicating that demographic bottlenecks during crop evolution may have generally occurred over substantial periods of time. Previous work has suggested population structure can generate spurious signals of population size change in methods like MSMC [50, 51], such that individuals sampled from a single deme of a highly structured population can falsely demonstrate signatures of a population bottleneck similar to what we observe in maize [51]. Given that our maize landraces are sampled from broad ecogeographic regions, however, this effect should be minimal. Moreover, a similar analysis in an Americas-wide sample of maize landraces demonstrated qualitatively similar results [2].

In addition to a strong bottleneck during domestication, our finding that levels of diversity decline in populations increasingly distant from the center of maize domestication are suggestive of serial founder effects during the spread of maize across the Americas (Fig. 1c; Additional file 1: Figure S1). Serial founder effects are the result of multiple sampling events during which small founder populations are repeatedly drawn from ancestral pools, leading to a stepwise increase in genetic drift and a concomitant decrease in genetic diversity. During maize range expansion, serial founder effects would have occurred if seed carried to each successive colonized region was limited to a sample of whole ears that contained a fraction of the diversity found in the source population [29]. Movement of entire ears involves a collective transfer of seeds that are either full or maternal half-siblings and could lead to more substantial founder effects than would be seen if dispersal were truly random. Such "kin-structured" migration, which is common in nature, has theoretically been demonstrated to increase inbreeding due to a reduction in the number of effective colonists [52]. Consistent with serial founder effects, other researchers have found a correlation between geographic and genetic distance in maize landraces [22, 53], though this was previously attributed to limited migration across the species range leading to isolation by distance (IBD). Neutral expectations of allele frequencies across populations under serial founder effects differ substantially from those predicted under equilibrium conditions [15]. For example, Slatkin and Excoffier [15] have demonstrated that allele frequency clines previously attributed to adaptation could be generated entirely by neutral processes under expansion. Many of the world’s crops have experienced such histories of expansion, and studies attempting to identify loci underlying crop adaptation during post-domestication spread to new environments may most accurately compare empirical data to neutral expectations under a serial founder effects demography [15].

While a history of serial founder effects partially explains the variation in diversity across maize landraces, there are deviations from this model. For example, our combined results showing increased ROHs (Additional file 1: Figure S3), lower nucleotide diversity (Additional file 1: Figure S2), and smaller effective population size (Fig. 1) in the Andes all suggest a pronounced, ancient founder event and are in agreement with previous work modeling demography in this region [54]. The founder event in the Andes may reflect initially limited cultivation due to the poor performance of maize in this region relative to established root and tuber staples [55]; maize cultivation may have only become widespread after an initial lag period necessary for adaptation. Additionally, we observe somewhat higher than expected nucleotide diversity in maize landraces from the highlands of Mexico and Guatemala (Fig. 1c), which may be linked to the introgression we have detected from mexicana.

In striking contrast to the bottleneck we observe in maize, the effective population size in parviglumis increases steadily from the time of initial maize domestication until the recent past. Multiple population genetic studies have reported negative genome-wide values of Tajima’s D in parviglumis from the Balsas River Valley [2, 23, 56], findings characteristic of an expanding population. Likewise, analyses of pollen content in sediment cores from Mexico suggest herbaceous vegetation and grasslands have expanded over the last 10,000 years due to changing environmental conditions during the Holocene and human management of vegetation with fire [57, 58]. While our parviglumis samples are drawn from a single population in the Balsas, these data collectively suggest parviglumis from this region has experienced expansion over the last several millennia.

Consistent with archaeological evidence of the timing of initial maize domestication [21], we find that maize demographies begin to diverge ≈ 10,000 generations BP, a time that appears to coincide with a steeper decline in maize N _e as well. In contrast, we estimate the timing of the split between maize and our single population of parviglumis to be ≈ 75,000 generations BP, potentially reflecting population structure in parviglumis. Beissinger et al. [2], using samples from additional populations, also find an estimate of maize-parviglumis divergence older than the probable onset of domestication, suggesting that currently available sequences of parviglumis may not sample well from the populations directly ancestral to domesticated maize.

The prevalence of gene flow during maize diffusion

Increasingly, range-wide analyses of crops and their wild relatives have identified evidence of gene flow during post-domestication expansion from newly sympatric populations of their progenitor taxa and closely related species [59–61]. Consistent with previous results from genotyping data [5, 22, 62], we find strong support for introgression from mexicana to maize in the highlands of Mexico. While mexicana is not currently found in the highlands of Guatemala, we also find strong evidence for mexicana introgression in maize from this region, suggesting either mexicana was at one time more broadly distributed, or, perhaps more likely, that highland maize from Mexico was introduced to the Guatemalan highlands. Support is also found for mexicana introgression in the southwestern USA at specific chromosomal regions such as a putative inversion polymorphism on chromosome 3 (Fig. 2). These results confirm previous findings suggesting maize from the highlands of Mexico originally colonized the southwestern USA [36]. The more limited signal of mexicana introgression here may be due to subsequent gene flow from lowland maize as suggested by [36]. Very little evidence is found for mexicana haplotypes extending into South America, as highland-adapted haplotypes would likely have been maladaptive and removed by selection as maize traversed the lowland regions of Central America [54].

Impacts of demography on accumulation of deleterious variants

Previous work in maize has characterized genome-wide trends in deleterious alleles across modern inbred maize lines, revealing that inbreeding during the formation of modern lines has likely purged many recessive deleterious variants [63] and that complementation of deleterious alleles likely underlies the heterosis observed in hybrid breeding programs [63, 64]. Additionally, [2] revealed that purifying selection has removed a greater extent of pairwise diversity (θ _π ) near genes in parviglumis than in maize due to the higher historical N _e in parviglumis, but that this trend is reversed when considering younger alleles due to the recent dramatic expansion in maize population size. To date, however, few links have been made between the historical demography of maize domestication and expansion and the prevalence of deleterious alleles. Our analysis reveals, for the first time, that demography has played a pivotal role in determining both the geographic and genomic landscapes of deleterious alleles in maize.

Population size and deleterious variants

Previous studies have suggested a “cost of domestication” in which a higher burden of deleterious alleles is found in domesticates compared to their wild progenitors [12, 41, 43, 65, 66]. Consistent with these results, we detect an excess of deleterious alleles in maize relative to parviglumis (Fig. 3; Additional file 1: Figure S6; Additional file 1: Figure S7), which could be caused by two potential factors. First, reduced population size during the domestication bottleneck could result in deleterious alleles drifting to higher allele frequency. Second, hitchhiking caused by strong positive selection on domestication genes could cause linked deleterious alleles to rise in frequency [12, 65]. While we find support for the former in maize, we see little evidence of the latter. Recent studies have reported contrasting results regarding the effect of selective sweeps in patterning the distribution of deleterious alleles. For example, putative selective sweeps in cassava showed a paucity of deleterious alleles, a result that was attributed to purifying selection [67]. Sweep regions in grape exhibited an overall decrease in the number of deleterious alleles but an increase in the ratio of deleterious mutations to synonymous variants, a pattern suggesting deleterious alleles may have hitchhiked along with the targets of positive directional selection [46]. Finally, selective sweeps in Asian rice contained a roughly equivalent ratio of deleterious mutations to synonymous mutations when compared to neutral regions [68]. Clearly, further exploration is warranted to clarify the effect of selection on the distribution of deleterious mutations. In addition to the cost of domestication, we find a cost of geographic expansion that is likely driven by serial founder effects. The increase in deleterious alleles during expansion is most pronounced in the Andes and may be symptomatic of the extreme founder event we propose above.

Differences in the number of deleterious alleles between maize and parviglumis and non-Andean and Andean maize are more dramatic under a recessive model than an additive model. This trend may indicate that the bulk of deleterious alleles in maize are at least partially recessive, such that heterozygous sites contribute less to a reduction in individual fitness. Previous work in human populations has shown that the majority of deleterious mutations are recessive or partially recessive [69], and data from knock-out mutations in yeast have revealed that large-effect mutations tend to be more recessive [70]. Likewise, both theory and empirical evaluation across a number of organisms suggest that mildly deleterious mutations are likely to be partially recessive [71]. In maize, Yang et al. [63] have found that most deleterious alleles are at least partially recessive and note a negative correlation between the severity of a deleterious variant and its dominance. Our results thus match nicely both with previous empirical data and theoretical expectations showing that recent population bottlenecks should only show strong differences in load under a recessive model [7].

Introgression and deleterious variants

Very few studies have investigated the effects of introgression from a taxon with substantially different N _e on the genomic landscape of deleterious variants. The best example is found in the human literature, where confirmation has been found that introgression from Neanderthals with low ancestral N _e increased the overall mutation load in modern humans [9, 19]. We report here the opposite pattern in maize, as introgression appears to have reduced the proportion of deleterious variants. Nonetheless, the underlying interpretation is similar: the taxon donating alleles mexicana has had a larger ancestral N _e than maize [27], and introgressed haplotypes have thus experienced more efficient long-term purging of deleterious alleles.

Conclusions

We have demonstrated that demography during the domestication and expansion of maize across the Americas has profoundly influenced putative functional variation across populations and within individual genomes. More generally, we have learned that population size changes and gene flow from close relatives with contrasting effective population size will influence the distribution of deleterious alleles in species undergoing rapid shifts in demography. The significance of our results extends far beyond maize. For example, invasive species that have recently experienced founder events followed by expansion, endangered species subject to precipitous declines in N _e , species with a history of post-glacial expansion, and new species expanding their range will all likely show clear genetic signals of the interplay between demography and selection. This interaction bears importantly on the adaptive potential of both individual populations and species. By fully characterizing this relationship, we can better understand how the current evolutionary trajectory of a species has been influenced by its history.

Methods

Samples, whole genome resequencing, and read mapping

A total of 31 maize landrace accessions were obtained from the US Department of Agriculture (USDA)’s National Plant Germplasm System and through collaborators (Additional file 2: Table S1). Samples were chosen from four highland populations (Andes, Mexican highlands, Guatemalan highlands, and southwestern US highlands) and two lowland populations (Mexican and South American lowlands) (Fig. 1a). In addition, four open-pollinated parviglumis samples were selected from a single population in the Balsas River Valley in Mexico. DNA was extracted from leaves using a standard cetyltrimethyl ammonium bromide (CTAB) protocol [72]. Library preparation and Illumina HiSeq 2000 sequencing (100-bp paired-end) were conducted by BGI (Shenzhen, China) following their established protocols. the Burrows-Wheeler Aligner (BWA) v.0.7.5.a [73] was used to map reads to the maize B73 reference genome v3 (GenBank BioProject PRJNA72137) [74] with default settings. The duplicate molecules in the realigned bam files were removed with MarkDuplicates in Picardtools v.1.106 (http://broadinstitute.github.io/picard), and indels were realigned with the Genome Analysis Toolkit (GATK) v.3.3-0 [75]. Sites with mapping quality less than 30 and base quality less than 20 were removed, and only uniquely mapped reads were included in downstream analyses.

Demography of maize domestication and diffusion

The MSMC method [31], which models ancestral relationships under recombination and mutation and has been used in several plant species [46, 47], was utilized to infer effective population size changes in both parviglumis and maize. SNPs were called via HaplotypeCaller and filtered via VariantFiltration in GATK [75] across all samples. SNPs with the following metrics were excluded from the analysis: QD <2.0; FS >60.0; MQ <40.0; MQRankSum <−12.5; ReadPosRankSum <−8.0. Vcftools v.0.1.12 [76] was used to further filter SNPs to include only bi-allelic sites. Following these data filtering steps, our data set consisted of 49 million SNPs. SNPs were phased using BEAGLE v.4.0 [77] with SNP data from the maize HapMap2 panel [78] used as a reference. Only sites with depth between half and twice of the mean depth were included in analyses. In addition, the software SNPable (http://lh3lh3.users.sourceforge.net/snpable.shtml) was used to mask genomic regions in which reads were not uniquely mapped. The mappability mask file for MSMC was generated by stepping in 1-bp increments across the maize genome to generate 100-bp single-end reads, which were then mapped back to the maize B73 reference genome [74]. Sites with the majority of overlapping 100-mers mapped uniquely without mismatch were determined to be “SNPable” sites and used for the MSMC analyses. For effective population size inference in MSMC, we used 5×4+25×2+5×4 as the pattern parameter, and the value m was set as half of the heterozygosity in parviglumis and maize populations, respectively.

In order to explore the trend of genetic diversity away from the domestication center, the correlation between the percentage of polymorphic sites and the geographic distance to the Balsas River Valley (latitude 18.099138, longitude –100.243303) was examined via linear regression. Geographical distance in kilometers was calculated based on great circle distance using the haversine transformation [17]. The correlation between percentage of heterozygous sites and distance away from domestication center was also explored in the SeeDs data set. SNPs with more than 50% missing samples and samples with more than 50% missing genotypes were removed from the SeeDs data set.

Population structure, genetic diversity, and inbreeding coefficients

We first evaluated population structure using principal component analysis (PCA) with ngsCovar [79] in ngsTools [80] based on the matrix of posterior probabilities of SNP genotypes produced in Analysis of Next Generation Sequencing Data (ANGSD) v.0.614 [81], and then utilized NGSadmix v.32 [82] to investigate the admixture proportion of each accession. The NGSadmix analysis was conducted based on genotype likelihoods for all individuals, which were generated with ANGSD (options -GL 2 -doGlf 2 -SNP_pval 1e−6), and K from 2 to 10 was set to run the analysis for sites present in a minimum of 77% of all individuals (24 in 31). A clear outlier in the Mexican highland population was detected, RIMMA0677, a sample from relatively low altitude, which was suspected to contain a divergent haplotype. A neighbor-joining tree of SNPs within an inversion polymorphism on chromosome 4 that includes a diagnostic highland haplotype [5] was constructed with the R package phangorn [83]. The sample RIMMA0677 was not clustered with other highland samples, but embedded within lowland haplotypes (Additional file 1: Figure S11), so it was removed from further analyses.

The genetic diversity measures Watterson’s θ and θ _π were calculated in ANGSD [81] for each population. The neutrality test statistic Tajima’s D was calculated with an empirical Bayes approach [84] implemented in ANGSD by first estimating a global site frequency spectrum (SFS) then calculating posterior sample allele frequencies using the global SFS as a prior. The three statistics were summarized across the genome using 10-kb non-overlapping sliding windows.

Inbreeding coefficients for each individual were estimated with ngsF [85] with initial values of F _IS set to be uniform at 0.01 with an epsilon value of 1e−5.

In addition, SNPs were polarized using the Tripsacum dactyloides genome to assess the frequency of derived homozygous sites in each maize landrace population. T. dactyloides short reads were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (SRR447804–SRR447807), mapped to the B73 reference genome v3 [74] with BWA [73], and incorporated into SNP calling as described above.

Runs of homozygosity

SNPs were down-sampled to contain one SNP in a 2-kb window to identify segments representing homozygosity by descent (i.e., autozygosity) rather than by chance. PLINK v.1.07 [86] was applied to identify segments of ROHs in a window containing 20 SNPs, among which the number of the maximum missing SNPs was set to 2 and the number of the maximum heterozygous sites was set to 1. The shortest length of final ROHs was set to be 300 kb. Physical distances were converted into genetic distances based on a recent genetic map [87].

Detection of introgression

To assess per-genome evidence of population admixture between maize landraces and teosinte, we calculated the D statistic using ANGSD [81]. The statistic was calculated using trees of the form (((X, low),mexicana),T. dactyloides). One accession from the Mexican lowland population was randomly sampled as the “low” taxon, and each sample from all other populations except the Mexican lowland was set as "X". The mexicana accession TIL25 from the maize HapMap2 project [78] was treated as the third column species. The D statistic was calculated in a 1-kb block, and then jackknife bootstrapping was conducted to estimate significance.

In addition, the $\hat {f_{d}}$ statistic [37] was calculated based on a similar tree form (((P ₁,P ₂),P ₃),O), but using allele frequencies across multiple individuals for each position on the tree. We fixed P ₁ as the Mexican lowland population, P ₃ as two lines of mexicana (TIL08 and TIL25), and T. dactyloides as the outgroup. P ₂ was set to each of the four highland populations and the South American lowland population.

The $\hat {f_{d}}$ statistic was calculated in 10-kb non-overlapping windows across the genome with the python script egglib_sliding_windows.py (https://github.com/johnomics/Martin_Davey_Jiggins_evaluating_introgression_statistics), which makes use of the EggLib library [88]. The input file was generated by first identifying genotypes using ANGSD (-doMajorMinor 1 -doMaf 1 -GL 2 -doGeno 4 -doPost 1 -postCutoff 0.95 -SNP_pval 1e−6) followed by format adjustments with a custom script (see “Availability of data and materials”). Outliers were detected by setting the 95% quantile of the $\hat {f_{d}}$ distribution in the South American lowland population as the cutoff.

Estimating burden of deleterious mutations

We estimated the individual burden of deleterious alleles based on GERP scores [89] for each site in the maize genome, which reflects the strength of purifying selection based on constraint in a whole genome alignment of 13 plant species [90]. The alignment and estimated GERP scores are available at iplant (https://doi.org/10.7946/P2WS60). Scores above 0 may be interpreted as historically subject to purifying selection, and mutations at such sites are likely deleterious. We identified Sorghum bicolor alleles in the multiple species alignment as ancestral and defined the non-Sorghum allele as the deleterious allele. Only biallelic sites were included for our evaluation. Inclusion of the maize B73 reference genome when calculating GERP scores [90] introduces a bias toward underestimation of the burden of deleterious alleles in maize versus teosinte populations. Therefore, we corrected the GERP scores of sites where the B73 allele is derived following [7]. Briefly, we divided SNPs where the B73 allele is ancestral into bins of 1% derived allele frequency based on maize HapMap3 [91] and used this frequency distribution to estimate the posterior probability of GERP scores for SNPs where the B73 allele is derived.

The sum of GERP scores multiplied by deleterious allele frequency for each SNP site was used as a proxy of individual burden of deleterious alleles under an additive model (H E T∗0.5+H O M∗1). This burden was calculated under a recessive model as the sum of GERP scores multiplied by one for each deleterious homozygous site (H O M∗1). For a better understanding of the variation of individual burden among sites under varied selection strength, we partitioned the deleterious SNPs into four categories (−2< GERP ≤0, nearly neutral; 0< GERP ≤2, slightly deleterious; 2< GERP ≤4, moderately deleterious; GERP >4, strongly deleterious) and recapitulated the preceding statistics.