The Genomic Basis for Short-Term Evolution of Environmental Adaptation in Maize.

Understanding the evolutionary capacity of populations to adapt to novel environments is one of the major pursuits in genetics. Moreover, for plant breeding, maladaptation is the foremost barrier to capitalizing on intraspecific variation in order to develop new breeds for future climate scenarios in agriculture. Using a unique study design, we simultaneously dissected the population and quantitative genomic basis of short-term evolution in a tropical landrace of maize that was translocated to a temperate environment and phenotypically selected for adaptation in flowering time phenology. Underlying 10 generations of directional selection, which resulted in a 26-day mean decrease in female-flowering time, [Formula: see text] of the heritable variation mapped to [Formula: see text] of the genome, where, overall, alleles shifted in frequency beyond the boundaries of genetic drift in the expected direction given their flowering time effects. However, clustering these non-neutral alleles based on their profiles of frequency change revealed transient shifts underpinning a transition in genotype-phenotype relationships across generations. This was distinguished by initial reductions in the frequencies of few relatively large positive effect alleles and subsequent enrichment of many rare negative effect alleles, some of which appear to represent allelic series. With these genomic shifts, the population reached an adapted state while retaining [Formula: see text] of the standing molecular marker variation in the founding population. Robust selection and association mapping tests highlighted several key genes driving the phenotypic response to selection. Our results reveal the evolutionary dynamics of a finite polygenic architecture conditioning a capacity for rapid environmental adaptation in maize.

AFTER ∼150 years of progress toward understanding evolution—since The Origin of Species (Darwin 1859)—burgeoning experimental results fueled by advances in genomic technology are shedding light on still unresolved questions about the nature of phenotypic change, including: the impact of mutation (e.g., Levy ) and standing variation (e.g., Burke ; Jones ); the role of epistasis (e.g., Tenaillon ); and the relationship between natural and artificial selection (e.g., Chan ). A key question, especially in the face of biological invasions and climate change, is how genomes confer and constrain the capacity for organisms to adapt to new environments (Orr 2005). Genetic dissection of experimentally evolved populations is a tractable framework for elucidating adaptive evolution since the experimenter can control selection and mating in particular environmental settings (Barrick and Lenski 2013; Schlötterer ). As a new extension to this framework, we implemented an efficient study design for dual inference about the population and quantitative genomic basis of phenotypic evolution (Wisser ). This was used to investigate the response to a decade of directional phenotypic selection for tropical-to-temperate adaptation in maize—a model species for plant genetics and a crop of global importance.

Genomic Basis of Response to Phenotypic Selection

The rate and history of mutations, the numbers and positions of functional variants, the distribution of allele effects, and the modes of gene action are among the genetic factors that shape the response to selection and influence the maintenance of phenotypic and genetic variability (Barton and Keightley 2002).

Considering theoretical population and quantitative genetic expectations for the response to directional selection, alleles at one or few loci with large effects on a selected trait should rapidly change in frequency, resulting in a corresponding phenotypic response (Falconer and Mackay 1996). As these alleles approach fixation, genetic variance is reduced and the response diminishes. Sustained responses may be attributed to standing polygenic variation, new mutations, epistatic interactions, or heritable epigenetic effects. For polygenic traits controlled by numerous loci of small effects, modeled at the extreme of infinite loci (Barton ; Fisher 1918), responses to selection can arise from subtle changes in allele frequencies across many loci. Consequently, allelic variation is retained and the causal-genic variance is expected to undergo negligible change. However, directional selection also creates negative disequilibrium covariance between allele effects across loci, resulting in temporary reductions in genetic variance for the trait under selection, a phenomenon referred to as the Bulmer effect (Bulmer 1971; Walsh and Lynch 2018). Qualitatively similar expectations arise under a so-called finite polygenic architecture where tens or more loci with allele effects of varying magnitudes are at play (Chevalet 1994; Fernando ; Turelli and Barton 1994).

An empirical understanding about the genetics of adaptation has been advanced through experimental population and quantitative genetic approaches (Savolainen ). The relative importance of genes with major and minor effects varies among traits, populations and species. At one extreme, relatively rapid or dramatic phenotypic changes have resulted from a few alleles at loci with large effects on traits such as flowering behavior (Lowry and Willis 2010) and toxin resistance (Baxter ). In contrast, other dramatic shifts in adaptive phenotypes have been ascribed to a polygenic architecture (Burke ; Berg and Coop 2014). Drawing a clear line of distinction between the two is not straightforward and is partially confounded by differences in experimental systems and their statistical power, but adaptation from a mixture of genes with major and minor effects have been reported (Levy ). Moreover, “evolve-and-resequence” studies have exposed unforeseen outcomes in the genomic changes underlying phenotypic evolution, including unique patterns of allele frequency change and the maintenance of molecular genetic diversity (Burke and Long 2012), the direct causes of which are unresolved.

Phenological Adaptation in Maize

Pivotal to adaptation and productivity in crop species is synchrony between the growing season and flowering time (Jung and Müller 2009). Numerous studies have investigated the genetic architecture of natural variation in flowering time for maize using a variety of methods, including genic analysis, linkage and association mapping, ecogeographical genetics and historical genetic analysis. Emerging from this body of literature is a consensus that allele effects dispersed across a finite polygenic architecture capture the major proportion of genotypic variation in flowering time (e.g., Chardon ; Buckler ; Li ), and that certain flowering time genes—Vgt1 (Salvi ; Ducrocq ), ZmCCT10 (Hung ; Yang ), ZmCCT9 (Huang ), and ZCN8 (Guo )—appear to have been instrumental to the postdomestication spread of maize from its tropical origin to many different environments. However, with one exception (Durand ), to our knowledge, the genomic basis of this adaptive trait has not been investigated in experimentally evolved populations. This could fill gaps in knowledge about the evolution of adaptation and lay a foundation to innovate breeding methods for rapidly adapting populations to new environments.

In this study, we investigated the genomic basis of adaptation from a distinct vantage point, where the entire period of evolution to an adapted state was captured in a single, multi-generational population—“Hallauer’s Tusón” (Teixeira ). Selection was initiated within an admixed founder population formed by intermating separate seed bank populations of Tusón, a landrace historically adapted to lowland tropical environments (Goodman and Brown 1988). Remarkably, 10 generations of directional phenotypic selection for early female-flowering time, with secondary selection for other traits in a temperate U.S. environment (Ames, IA; 42.03° N latitude), recapitulated the temperate-adapted state of maize achieved by early farmers; this is thought to have occurred over the course of thousands of years (Swarts ), albeit by less deliberate breeding methodology. A study design, in which families derived from genotyped individuals sampled across generations were phenotypically evaluated, allowed us to employ methods for dissecting both the population and quantitative genomic basis of phenotypic evolution (Figure 1). Because tropical maize is a rich resource of potentially useful genetic diversity that has largely been under used by maize breeders in temperate environments (Goodman 1998), findings from this study can guide future maize breeding for climate change and address fundamental questions about the genomic basis of environmental adaptation in plants.

Figure 1

Study design adapted from Wisser . Filled circles represent populations, filled squares represent individuals, and filled rectangles represent families. Multiple populations of a historically adapted landrace of maize, sampled from different tropical regions, were randomly intermated to form . Artificial truncation selection for early female-flowering time was performed in a novel temperate environment to which the source populations had not previously been exposed. Individuals were sampled from among the generations of selection, and genotyped. Allele frequency mapping was applied along the generational axis to characterize the population genetic basis of the phenotypic response to selection. The genotyped individuals were self-pollinated to generate families that were evaluated in a common garden experiment. Using the parental genotype data and corresponding family mean phenotype data, association mapping was applied to characterize the quantitative genetic basis of variation in the selected trait.

Materials and Methods

Front matter

Unless otherwise noted, data analysis was performed using R (R Core Team 2016); R packages are cited accordingly. The following abbreviations are used: AFPC (allele frequency profile cluster); BLUEs (best linear unbiased estimates); FDR (false discovery rate); FITR (frequency increment test with reference loci); GWA (genome-wide association); LD (linkage disequilibrium); SIM (simulation test statistic).

Plant material

The subject of this study was Hallauer’s Tusón, a multigenerational population of maize derived from a landrace of tropical origins that was subjected to 10 generations of phenotypic truncation selection for early female-flowering time in a temperate environment (Ames, IA; 42.03° N latitude) (Teixeira ; Hallauer and Carena 2016). Figure 1 depicts the breeding scheme for Hallauer’s Tusón and our study design.

As described by Hallauer and Carena (2016), the base population used to initiate selection was produced in Iowa by isolated, open pollination (no intentional selection) among multiple seed-bank accessions of the Tusón landrace sampled from different countries; however, the origins of some of these accessions remain unclear. It should be noted that maize is a monoecious species with female and male organs on separate parts of the same plant, such that open pollination includes the possibility of self-fertilization. Using the base population, selection ensued where plants were grown in isolation and allowed to intermate at random, which again included the possibility of selfing. During flowering, of the earliest female-flowering individuals (based on their silk-emergence phenotype), secondarily selected for other traits, were tagged and later harvested. An equal number of seeds per ear were mixed to form the subsequent generation of individuals for selection. This recurrent selection scheme was applied for 10 generations (from 1995 to 2004) until the population was deemed phenologically adapted by comparison with other temperate adapted maize.

Phenotype data

Previously, Teixeira phenotypically evaluated 297 self-pollinated families derived from the even numbered generations of Hallauer’s Tusón ( and ) for two years at multiple locations in North America, including the Iowa location where Hallauer’s Tusón was originally selected. Under the design depicted in Figure 1, the present study combined the available phenotypic data for female-flowering time measured in the selection environment (Ames, IA) and a highly correlated environment (Newark, DE) with new genotype data.

Genotype data

DNA was isolated from the 297 parents (noninbred) of the families that were evaluated phenotypically, plus an additional 90 random individuals from that were not phenotyped (this was done to provide a larger sample size of for more reliable genomic inference). Lyophilized leaf tissue was pulverized using a Geno/Grinder 2000 and extracted with the DNeasy 96 Plant Kit (Qiagen). Genotyping with Illumina’s MaizeSNP50 Beadchip (Ganal ) was performed by DNA LandMarks (Québec, Canada), producing genotype data at SNP sites with an average of 2.8% missing data per sample (min: 0.2%; max: 20.9%). Additional genotype data were produced for variant sites upstream of ZmCCT10 (Zm00001d024909) (Supplemental Material, Supplemental Methods), including a presence–absence causal variant for photoperiod sensitivity (Hung ; Yang ). See Supplemental Methods for genotypic data quality control (Table S1) and projection of markers onto the consensus linkage map for a maize nested association mapping population (McMullen ).

Analysis of genetic diversity

Specific subsets of markers were used in examining different aspects of genetic diversity (File S1 and Table S2). Some subsets use marker names prefixed with “PZA” or “PZE,” which have the lowest ascertainment bias among the MaizeSNP50 SNPs (Frascaroli ).

Summary statistics of genetic diversity were computed with hierfstat v. 0.04–14 (Goudet 2005) to calculate (average observed heterozygosity within generations), (average expected heterozygosity within generations) and (average expected heterozygosity for the total population), as well as (average inbreeding coefficient within generations), (average differentiation between generations), and (average inbreeding coefficient for the total population) according to Nei (1986)—this corresponds to the sample-level scope of inference. HardyWeinberg v. 1.5.5 (Graffelman 2015) was used to perform exact tests for Hardy–Weinberg equilibrium.

The relationship of Hallauer’s Tusón with maize more broadly was assessed using 934 samples representative of global maize germplasm for which genotype data for the same markers was available (Ganal ). A two-dimensional projection of relationships among samples was computed using PHATE (Potential of Heat-diffusion for Affinity-based Trajectory Embedding) (Moon ) implemented in phateR v. 0.2.7, using default settings with a precomputed simple matching distance matrix as input for the knn.dist.method.

Results from PHATE using only Hallauer’s Tusón suggested some structure was present among samples from but not other generations. Therefore, STRUCTURE (Pritchard ) was used to examine subpopulations among samples assuming admixture and independent allele frequencies (Table S3; see Supplemental Methods for details). The method informed our selection of K (Evanno ). Although the exact accessions of Tusón used to form remain unclear, several were used (Hallauer and Carena 2016). Therefore, we ignored the stronger signal of suggesting K = 2 and chose the next highest peak at (Figure S1). To examine whether any one subpopulation was favored in the initial generations of selection, STRUCTURE was also used to ascertain the ancestry present in samples, assuming the aforementioned subpopulation classification for individuals in . A paired t-test was used to test the difference between the proportion of individuals per subpopulation in (Table S3) and the average per individual admixture proportion estimated by STRUCTURE for (Table S4).

LD (Hill and Robertson 1968), measured as , was computed with genetics v. 1.3.8.1 (Warnes ). The structure of LD between chromosomes was characterized per generation using low-ascertainment biased markers with no missing data and a minimum allele count of 12 in the given generation (Table S2). The structure of LD within chromosomes was similarly characterized, but with a larger number of markers (Table S2) and examined at different intervals of genetic distance: (0,1],(1,5],(5,10],(10,50],>50. In addition, was computed among generations (all samples) between sequential pairs of markers per chromosome. Using these latter estimates of LD, lokern 1.1.8 (Herrmann 2016) was used to perform kernel regression of the values as a function of the midpoint basepair coordinate between each pair of markers. A plot of pairwise between markers flanking the ZmCCT10 associated causal site for photoperiodism was made using LDheatmap v. 0.99.2 (Shin ).

Population genetic analysis

Allele frequency mapping:

Three statistical tests were used to detect markers with nonrandom patterns in allele frequency change across generations: (i) a customized whole genome SIM test for departures from genetic drift; (ii) the Bayenv test for robust correlations between allele frequencies and generations (Coop ); and (iii) the FITR for robust departures from genetic drift (Nishino 2013).

A detailed description of the SIM test can be found in Supplemental Methods. Briefly, by customizing simuPOP (Peng and Kimmel 2005), a simulator was constructed to generate in silico genomes (genome-wide genotypes) for individuals constituting a population that undergoes breeding according to the design for Hallauer’s Tusón (with random selection to model genetic drift) and sampling according to our study design (to account for sampling variance). The simulator used the fixed genotype matrix, fixed recombination rates estimated from projection onto the genetic map (Supplemental Methods) and the fixed STRUCTURE matrix (Table S3) to initially generate random in silico genomes (derived based on the structured sample data) from which simulated breeding ensued. At each marker across the genome, the probability of the observed sample allele frequency change was computed relative to the expected distribution created from replicates of simulation. Marker p-values were adjusted for multiple testing (Benjamini and Hochberg 1995), and those with a FDR are referred to as markers, while the remaining markers are referred to as .

Bayenv 2.0 was used to identify robust correlations between allele frequencies (response variable) and generations (explanatory variable). Bayenv was originally designed to test for correlations with an environmental variable; here, generation numbers (standardized) were used instead. The covariance matrix used to model the background expectation of allele frequency change was estimated from the low-ascertainment biased SNPs (Table S2). Markers with the top 1% Bayes factor values were considered robust outliers, which we refer to as markers.

The FITR statistic is conditioned on the variance in allele frequency change estimated from neutral markers in the sample (Nishino 2013). For this, we used markers declared at a 10% FDR. To minimize bias from using a single set of reference markers, we modified the test by bootstrap resampling markers and computing the proportion of times in which each marker was significant (at a 1% FDR) among bootstrap samples. In each bootstrap sample, the reference data comprised 1% of the markers (n = 444) sampled at random and without replacement. The test is not definable for fixed sites, yet the alleles at markers used for testing were not always observed in each generation. Therefore, alleles with an observed frequency of zero in a particular generation were set to , where is the number of individuals in the generation. Markers with a FDR in at least of the bootstrap sample tests are henceforth referred to as .

Localizing footprints of selection:

Chromosomal regions with a local footprint of selection (i.e., deviation from genetic drift across a segment of the genome) were delimited using chromosome-specific kernel regression functions of the values from the SIM test on the physical coordinates of markers. To obtain a definable input value for markers where the SIM test p-value equaled 0 (where the observed data fell outside the limits of the null distribution of simulated drift), was set to 0.00057, which was half the minimum p-value among all SIM tests. Regions along each chromosome where the kernel regression line surpassed a threshold of were considered local footprints of selection, and are henceforth referred to as .

Characterizing features of allele frequency change:

Divisive analysis of hierarchical clustering (Kaufman and Rousseeuw 1990) was performed with cluster v. 2.0.7 (Maechler ) in order to group markers with similar profiles of allele frequency change for the minor allele in . The number of clusters was determined using clusGap based on the Tibs2001SEmax criterion (Tibshirani ).

Additional data summaries, using the minor allele in as the reference allele, were used to compare features of allele frequency change between and markers, including: (a) the slope and intercept from regression of allele frequency change on generations; (b) the mean absolute change in allele frequency among the highest to lowest ranking changes in frequency per marker across sequential pairs of generations (the largest amount of change between a given pair of generations was assigned a ranking of 1; the least amount of change was assigned a ranking of 5); (c) the distribution of the longest run (i.e., number of generations) of positive and negative monotonic change in allele frequency across sequential generations; and (d) the rank distribution for the amount of allele frequency change across generation pairs.

Quantitative genetic analyses

Genetic differentiation:

, a measure of the proportion of genetic variance distributed among populations for quantitative traits (Spitze 1993), was used to estimate genetic differentiation in female-flowering time between and each subsequent generation of Hallauer’s Tusón, as well as to examine relative to the distributions of for and markers. Knowing that flowering time was under selection, the comparison was used to characterize the relationship between population genetic and quantitative trait divergence (Le Corre and Kremer 2012). Following from Spitze (1993):where and are estimates of the among-generation and average within-generation additive genetic variances, respectively. See the next section and Supplemental Methods for details on the estimation of variance components. For comparison, we used the Hudson estimator, (Bhatia ), which is compatible with restricted maximum likelihood estimation of genetic variances used for , since both estimates correspond to the population-scope of inference in the broad sense.

Partitioning of the genotypic variance:

Using ASReml v. 3 (Gilmour ), the following mixed linear model was used to partition the phenotypic variance and decompose the genotypic variance into additive, dominance, and residual genetic variance components:where y corresponds to the vector of observations (female-flowering time), is the fixed overall mean effect, and and correspond to the vectors of random environment effects, incomplete block nested in replication within environment effects, additive genetic family effects, dominance genetic family effects, residual genetic family effects, family × environment interaction effects, and residuals, respectively. The effect of replication nested in environments was excluded as it was not significant according to a likelihood ratio test. The respective design matrices [ (this has same structure for , and ] relate observations to their corresponding vectors of effects. The additive, dominance, and residual genetic family effects were assumed to be distributed independently of one another, where: . The G matrix was computed according to VanRaden (2008) and the D matrix according to Su , while I is an identity matrix.

Variance component estimates from Equation 2 were used to compute heritability in the broad and narrow sense on an entry mean-basis according to Holland . In addition, extensions of Equation 2 were used to examine the amount of genetic variance in female-flowering time explained by each chromosome and for vs. markers (Supplemental Methods).

Genome-wide association mapping:

Kinship-controlled GWA was performed using a standard two-step procedure (Yu ). First, BLUEs of female-flowering time for the families were estimated using Equation 2, but replacing additive, dominance, and independent random genetic effects with a single term for families fit as a fixed effect. Second, using BLUEs as the response variable, markers were tested for trait association using the mixed linear model in TASSEL (Bradbury ) standalone v. 5.2.12, while controlling for the random polygenic background with the aforementioned G matrix.

To reduce false-positive associations due to rare genotypes co-occurring with outlier phenotypes, if the sample size of phenotyped lines for a given genotypic class at a marker was less than five, individuals with the corresponding genotypic state (typically the homozygous minor allele class) were set to missing for that marker. The QQ plot of GWA p-values is shown in Figure S2. A 10% FDR was used to declare significant trait-marker associations; henceforth, these are referred to as markers.

When estimating additive allele effects, some markers had only two genotypic classes (heterozygous and one homozygous class). The effects of minor variants at these loci were estimated as the difference between the heterozygous class and the homozygous class. For markers with three genotypic classes, the additive effect was uniformly reported as half the difference between the homozygous variant class corresponding to the minor variant in and the homozygous variant class of the alternative variant.

Using the closest-features program of BEDOPS (Neph ) v. 2.4.15, the flowering time candidate gene [Dataset S8 in Hung ] nearest to each marker was determined.

Synthesis map

A graphical map of the maize genome integrating multiple results was created. The map included local linkage disequilibrium estimated by kernel regression, the difference in heterozygosity between generations 0 and 10 , SIM test results (-transformed q-values), Bayenv test results (-transformed Bayes factor values), FITR test results (bootstrap values >75%), GWA test results (-transformed q-values), previously mapped QTL associated with flowering time {flowering time per se [Table S3 in Buckler ] and photoperiod sensitivity [Dataset S3 in Hung ]}, and candidate genes for flowering time [Dataset S8 in Hung ].

Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article. Supplemental methods, tables, and figures are available at figshare: https://doi.org/10.25386/genetics.9936284. Supplemental Tables: Tables S1 and S2 details results from quality control filtering and lists subsets of the genotype data used for analysis; Tables S3 and S4 contain results from STRUCTURE analysis; Table S5 summarizes identified by kernel regression; Table S6 shows chromosome-specific genetic variance component estimates; Table S7 lists the markers detected by GWA; and Table S8 shows the candidate gene for flowering time nearest to each marker. Supplemental Figures: Figure S1 shows results from STRUCTURE used to select K; Figure S2 is a QQ plot of observed vs. expected p-values for GWA tests; Figure S3 is a Venn diagram of markers detected when the SIM test was applied to sequential pairs of generations; Figure S4 summarizes various features of allele frequency change for and markers; Figure S5 shows the allele frequency profile clusters and corresponding distributions of additive allele effects for markers; Figure S6 shows the structure of LD within and between chromosomes for pairwise combinations of and markers; Figure S7 shows the synthesis map of multiple analysis results. Supplemental Files: File_S1.txt contains a list of the quality control markers, their map locations on B73 AGPv2 and AGPv4 reference assemblies, and the analysis-specific, Table S2 subset to which they belong. File_S2.txt contains summary statistics and test results for each marker, including: allele frequencies per generation for the corresponding minor allele in ; observed heterozygosity per generation; between and ; p-values and q-values for the SIM test; Bayes factor values and correlation statistics for Bayenv, bootstrap values for FITR; p-values and q-values for GWA; and the estimated additive allele effect. Phenotype and genotype data are available at Dryad: phenotype - https://doi.org/10.5061/dryad.8f64f; genotype - https://doi.org/10.5061/dryad.q573n5tdt. Python code used for genome simulation is available via GitHub: https://github.com/maizeatlas/saegus.

Results

Artificial selection generated a tropical genome with a temperate-adapted phenome

Hallauer’s Tusón population was founded by intermating multiple seed bank accessions of a maize landrace historically adapted to tropical environments (Figure 1). Teixeira demonstrated the capacity of this population to become phenologically adapted to a temperate environment within 10 generations of artificial selection, based primarily on selection for early female-flowering time. Here, we found that the population was highly diverse; nearly the entire set (96%) of SNPs on the MaizeSNP50 chip (Ganal ) segregated within or among generations (Table 1).

Table 1

Summary statistics for molecular genetic diversity in Hallauer’s Tusón

Generation	Sample size	Proportion polymorphica	H_Ob	H^_Sc	H^_T^′d	F^_ISe	HWD^+f	HWD^+:H_O≤H^_Sg
0	105	0.897	0.285	0.320	NA	0.111	0.256	0.784
2	56	0.893	0.323	0.315	NA	−0.026	0.071	0.400
4	55	0.902	0.325	0.314	NA	−0.037	0.067	0.373
6	54	0.914	0.315	0.320	NA	0.016	0.072	0.475
8	56	0.923	0.328	0.325	NA	−0.006	0.072	0.440
10	55	0.929	0.333	0.327	NA	−0.017	0.069	0.411
0−10	381	0.963	0.318	0.320	0.305	0.021	0.262	0.704

Results in the column are based on 49,477 markers (Table S1), while all remaining results are based on markers with a minor variant count among all samples.

Average observed heterozygosity within generations. For all samples, this corresponds to the average among generations.

Average expected heterozygosity within generations. For all samples, this corresponds to the average among generations.

Average expected heterozygosity for the total population.

Average inbreeding coefficient within generations. For all samples, this is , the average inbreeding coefficient for the total population.

The proportion of markers that significantly deviated from Hardy–Weinberg equilibrium determined at a 5% FDR.

The proportion of markers where was less than .

When compared to a global sample of maize, all of the individuals across generations of Hallauer’s Tusón clustered with tropical germplasm (Figure 2a), notwithstanding the temperate-adapted phenome of individuals belonging to the later generations (Teixeira ). Thus, to tackle challenges associated with crop vulnerability through plant breeding, Hallauer’s Tusón highlights the adaptive potential of maize landrace populations, and provides a unique source of germplasm likely to contain novel alleles for temperate maize breeding programs.

Figure 2

Genetic diversity and population structure in Hallauer’s Tusón. (A) 2D PHATE plot showing relationships for Hallauer’s Tusón and a broad range of maize germplasm. Aside from distinguishing samples from Hallauer’s Tusón and teostine inbreds, population structure groups were assigned to samples based on Table S1 in Romay (several samples were unclassified but retained in the analysis). (B) Population structure in . Admixture profiles for six subpopulations are shown. Values at the top of the STRUCTURE plot correspond to days to female flowering time for individuals with phenotype data. Subpopulation frequencies of ZmCCT10-s, the deletion (“lack-of-insertion”) variant associated with photoperiod sensitivity is shown.

Retrospective analysis reveals admixture during selection on a structured founder population

Genomic inference combined with knowledge about the population development scheme was used to illuminate the breeding history of Hallauer’s Tusón. Although generations showed some SNP-based differentiation , as expected for progeny generated by random mating among selected individuals, overall inbreeding was minimal . However, on a per generation basis, was noticeably higher in , where of the markers significantly deviated from Hardy–Weinberg equilibrium, with of those deviations being due to an excess of homozygotes (Table 1). This suggested a Wahlund effect (Hartl and Clark 2007) from sampling geographically separated accessions that remained after the intermating step to form , but subsequent generations showed evidence of random mating during selection (Table 1).

STRUCTURE analysis also indicated that individuals in had formed largely by intermating within separate founder accessions and hybridization between specific pairs of accessions (Figure 2b and Table S3). These findings are congruent with the way in which was bred, whereby the original Tusón accessions were planted in adjoining blocks and allowed to open-pollinate, which would have favored mating within and between pairs of subpopulations. Although selection could potentially favor specific subpopulations under these conditions, the genomic ancestral composition for individuals showed admixture profiles that were proportional to that of the subpopulation sizes in (Table S4). Thus, randomized bulking and planting of seed between each generation of artificial selection minimized subsequent inbreeding and population structure during selection.

Differentiation across a fraction of the genome potentiated strong phenotypic change

A decade of directional phenotypic selection, resulting in an overall mean decrease of 26 days to female-flowering time, caused generations to become strongly differentiated phenotypically, as measured by (Figure 3). Simulation of neutral allele frequency changes that would occur under the breeding scheme used for phenotypic selection allowed us to identify 6115 of markers with non-neutral allele frequency changes (referred to as markers). These markers were widely dispersed across the genome (but with some clusters of linked markers, as described later), and were distinguished from markers by their increasing levels of across generations relative to (Figure 3). The identification of a sizeable fraction of genome-wide markers as suggested a finite polygenic architecture (i.e., possibly tens to hundreds of loci affecting flowering time) could underlie the phenotypic response to selection. Similarly, based on comparisons, the very large increases in could be explained by much smaller levels of amplified across a large number of loci.

Figure 3

comparisons. (vertical lines) for female-flowering time is compared to the distributions of (boxplots) for and markers. Color coding is according to the pair of generations that were compared as indicated on the y-axis. Black points show of the ZmCCT10_CACTA marker which reached its maximum value of 0.25 by , the generation by which ZmCCT10-s was purged from the population.

Population genetic analysis pinpoints shifts in the genetic architecture underlying response to selection

Although increased overall across generations at markers, changes in the frequencies of alleles at a locus varied mostly among generations. Using the simulator to test for non-neutral allele frequency changes between sequential pairs of generations showed that a majority of marker-specific departures were exclusive to one consecutive pair (Figure S3). This was coincident with the common observation of “bursts” in allele frequency change within a few generations, rather than monotonic changes across all generations (Figure S4, A–C).

Because transient changes in allele frequency were a prominent feature of the genomic response to selection, we used clustering to examine the temporal structure of allele frequencies among markers. A total of 15 allele frequency profile clusters (AFPCs) were resolved with a high degree of overall clustering structure (divisive coefficient = 0.98; Figure S5). Minor alleles in with negative frequency trajectories were captured in AFPCs comprising of the markers. The remaining minor alleles, however, were enriched from starting frequencies that spanned the minor-allele frequency spectrum. Across clusters, some notable transitions in allele frequency responses were observed: (i) alleles in AFPC1, which included the photoperiod sensitive allele (ZmCCT10-s), substantially reduced in frequency to become rare or removed within the first four generations; (ii) several AFPCs (6, 10, 11, 13, 14, and 15) showed clear increases in allele frequency within the first few generations that were limited thereafter; and (iii) for AFPC4, a dominant cluster comprising of all markers, initially rare alleles appreciably increased after .

There were few instances, like ZmCCT10-s, where alleles were purged from the population. The loss of SNP variants in Hallauer’s Tusón was actually rare— of all SNPs in were purged during selection (10 of these were ), and this was 3.4 times less than the average proportion of SNPs purged across replicate simulations of neutral frequency change . Taken together, population genetic analysis suggested that selected genotype–phenotype relationships temporally shifted across a finite polygenic architecture with a predominant enrichment of initially minor alleles, resulting in the maintenance of genetic variation.

Quantitative genetic analysis contextualizes genome-wide population genetic dynamics

Our experimental design included phenotypic data on self-pollinated families of genotyped individuals sampled across generations evaluated in common environments (Figure 1), permitting estimation of quantitative genetic parameters and interpretation in the context of population genetic results. Underlying a 50-day range for female-flowering time in Hallauer’s Tusón, high broad and narrow sense heritabilities indicated a large fraction of the phenotypic variance could be explained by genotypic effects, and the genotypic variance partitioned into additive and dominance variance with no residual genetic variance remaining. Including an additive-by-additive epistatic relationship matrix did not improve the model fit nor did it explain any of the genotypic variance.

The genetic variance explained by individual chromosomes varied widely (Table S6). For instance, chromosome 10, in which the ZmCCT10-associated causal variant is located, was an outlier that accounted for a large proportion of the additive variance, while chromosome 2 included no additive or dominance variance. This supports and extends the population genetic inference of a finite polygenic architecture, showing variability in the genetic effects across the genome available for selection.

Genome-wide additive allele effects (estimated among families across generations) correlated with average changes in allele frequency per generation , whereby, as expected, alleles with negative effects on female-flowering time (contributing to early flowering) tended to have positive slopes in allele frequency change and vice versa (Figure 4; the Bayenv and GWA hits highlighted in the Figure are discussed in the next section). This relationship was largely driven by markers, which accounted for of the additive variance (none of the dominance variance)—an excessive enrichment given these constituted of the SNPs.

Figure 4

Genome-wide relationship between slopes in allele frequency change and additive allele effects. Each point corresponds to the minor allele in and indicates the slope in allele frequency change estimated across all generations (y-axis) vs. the additive allele effect estimated among genotypes from all generations (scaled in days; x-axis). Marker effects are estimated separately, and do not account for covariances among loci, such that their combined effects tend to be less than the sum of individual effects. Marginal histograms show their distributions. The coefficient of determination and slope of the relationship are shown. Color coding depicts test statistic results ( FDR GWA hits; top Bayenv hits; and FDR q-values for the SIM test). The ZmCCT10-s allele is marked by an open circle.

Nonlinear changes in the phenotypic mean and additive variance for female-flowering time across generations reflected some of the observed dynamics in allele frequency change (Figure 5). Changes in the mean could be modeled as a cubic function with significant coefficients , where female-flowering time decreased across all generations but to a larger degree between generations and . These generations contained greater numbers of markers with larger magnitudes of allele frequency change (Figure S4d). The corresponding change in additive variance could be modeled as a quadratic function of generations with significant coefficients , where initially the variance was greatly reduced but later increased across generations ; by , the additive variance exceeded that of . These changes were coincident with initial reductions in the frequencies of few relatively large positive effect alleles (AFPC1) and subsequent enrichment of many rare negative effect alleles (AFPC4) (Figure 5 and Figure S5).

Figure 5

Quantitative and population genetic components of the response to selection. (A) Zero-centered (not standardized) “Z” values corresponding to the mean (BLUE) and additive variance (Va) for female-flowering time per generation. (B) Box plots of allele frequencies per generation for markers in AFPCs 1 and 4. (C) Histograms of additive allele effects for markers in AFPCs 1 and 4. (D) SIM hits (vertical lines colored by chromosome) within (x-axis) that belonged to AFPCs 1 and 4. The silhouette width (y-axis) is a measure of a markers fit to the cluster (values can be considered a good fit). Facet labels indicate the AFPC identifier and the proportion of markers per cluster [in (B) and (C) this corresponds to the proportion among all SIM hits; in (D) this corresponds to the proportion among regional SIM hits].

Robust selection and association mapping identify associations with key flowering time genes

The simulation test clearly enriched for markers that differentiated generations (Figure 3), but not all of these are necessarily linked to causal variants underlying female-flowering time selection. For instance, chromosome 2 captured none of the genetic variance in female-flowering time; however, it contained markers (Table S6). The model used to simulate breeding events does not include components of potentially important sources of variation, such that departures from the null distribution may also be due to factors such as individual differences in gametic fitness, or the secondary selection that was exerted for other traits (Hallauer and Carena 2016). On the other hand, covariance between causal and neutral allele frequency changes may generate false positives (Coop ). Although we could not control for deviations due the former, selection and association mapping tests that control for genomic background effects were used to identify markers exhibiting robust changes in allele frequency across generations and robust associations with trait variation among generations, respectively.

Prior to applying these tests, we assessed the transgenerational structure of LD across the genome using subsets of low-ascertainment biased markers with a standardized minimum allele count per generation (Table S2). Linkage disequilibrium was examined for each combination of and markers. Based on all pairwise estimates of LD, median showed little variation across generations, and did not exceed 0.04 within chromosomes and 0.01 between chromosomes (data not shown). At increasingly higher percentiles of the distribution, LD between pairs of markers within chromosomes (but not between chromosomes) was elevated relative to other combinations of and markers (Figure S6). We found that the heightened LD between markers was restricted to local linkage blocks, and tended to increase across generations (Figure S6)—a hallmark footprint of selection.

Given that the structure of LD gave rise to nonindependent sets of linked SIM hits, kernel regression was used to delimit 29 encompassing 1008 of the markers (Table S5). The Bayenv test, which controls for genome-wide covariance in sample allele frequencies, produced Bayes factor values that were correlated with the SIM test . The top 1% of Bayenv hits were located on all chromosomes and were present in most , with regions on chromosome 9 being heavily populated with markers (Figure 6 and Figure S7). Similarly, the FITR test, which is conditioned by variance in allele frequency change estimated from the sample, primarily implicated on chromosome 9 as robust outliers.

Figure 6

Synthesis map of chromosomes 8, 9, and 10. Multiple results are plotted on the physical map of each chromosome, with the y-axis corresponding to values for each of the following metrics: (i) kernel regression estimate of for LD between sequential pairs of markers (black line); (ii) kernel regression estimate of for the SIM test (orange line: delimited are enumerated and encompass the orange shaded areas); (iii) value for markers (orange vertical lines); (iv) value for complete-sweep markers (orange filled box); (v) difference in observed heterozygosity between and for markers in (black-filled triangles: pointing up if the change was positive and down if the change was negative); (vi) (Bayes factor) values for markers (cyan-filled points); (vii) bootstrap values for markers (blue-filled points); (viii) value for markers (red-outlined points); (ix) QTL previously identified for photoperiod sensitivity (gray shaded areas corresponding to QTL intervals) and flowering time per se (green vertical lines corresponding to QTL peaks); (x) candidate genes for flowering time (yellow vertical lines and labels); and (xi) centromeres (lilac-colored boxes).

GWA mapping performed on mean female-flowering time resulted in few genome-wide significant associations, notwithstanding the SNP-based polygenic model that explained essentially all of the genotypic variance. Between 2 and 12 GWA hits were detected across 1–10% FDR thresholds (Table S7). All but one of these showed the expected relationship between the sign for the additive allele effect and slope in frequency change (we note this one marker was on chromosome 2, which explained none of the genetic variance when the whole chromosome was modeled under a polygenic architecture). However, no GWA hits were detected on chromosome 9 nor within any of the . The top GWA hit was the presence–absence causal variant for ZmCCT10-regulated photoperiodism, which was also detected as hits by the SIM and Bayenv tests but not the FITR test. Otherwise, the strength of signals (slopes vs. effects) for top hits by selection and association mapping tests tended to differ (Figure 4).

Taken together, robust tests to dissect the genetic basis of the response to selection implicated a number of genes previously associated (causally or as a candidate gene) with variation in flowering time and photoperiodism in maize, several of which are highlighted in Figure 6.

Evidence for multiple local haplotypes underlying the phenotypic response to selection

The transition from selection on common to rare alleles occurred at some of the same regions of the genome. For instance, on different chromosomes included markers in both AFPC1 and AFPC4, in which allele frequencies showed strong shifts during different periods of selection (Figure 5). Similarly, at the ZmCCT10 locus, robust associations were detected for SNPs that responded to selection even after the elimination of ZmCCT10-s. Pairwise LD between significant markers at the ZmCCT10 locus indicated two separate haplotypes were responsive to selection (Figure 7). These results reinforce the conclusion of a finite polygenic architecture underlying the response to selection, and extend that to suggest selection on multiple local haplotypes was an important aspect of short-term evolution.

Figure 7

Allele frequency change and linkage disequilibrium at ZmCCT10. (A) Variant frequency change (y-axis) across generations (x-axis) for the nearest five markers flanking the ZmCCT10_CACTA causal site. Points are color coded according to the combination of test results (ZmCCT10_CACTA was positive for the SIM, Bayenv and GWA tests and is color coded red). (B) Pairwise LD between the markers in (A). Marker distances from the insertion site for ZmCCT10_CACTA (AGPv4 94,435,768) are indicated. Negative values are in the direction of the ZmCCT10 transcription start site.

Discussion

Genetic analysis of adaptation in crop species provides a lens into evolution and generates relevant information for plant breeding. Although flowering time phenology has been widely studied in plants (Jung ), we are aware of no study (in plants) that has dissected the transgenerational genomic basis of adaptive evolution (here, for flowering phenology) in a population translocated to a new environment. We help close this knowledge gap by investigating a tropical landrace of maize that was adapted to a temperate environment across a decade of artificial selection (Teixeira ; Hallauer and Carena 2016). Using an efficient study design (Figure 1; Wisser ), we simultaneously elucidated population and quantitative genetic components underlying the 10 generations of selection required for the population to reach a state of phenological adaptation similar to modern temperate maize lines.

The evolutionary capacity of the tropical Tusón landrace to become rapidly adapted to a temperate environment was attributed to a finite polygenic architecture, yet two genomic phases underlying the phenotypic response to selection could be discerned. The first phase, from generations was distinguished by an oligogenic-like architecture, where marked reductions of a relatively small number of moderate-frequency minor variants in (AFPCs 1 and 2), with relatively large positive effects on flowering time, contributed to an initial strong response to phenotypic selection and a large reduction in genetic variance (Figure 5). Afterward, the genomic basis of the response transitioned to become dominated by the enrichment of a large number of rare-minor variants in with smaller-sized negative effects on flowering time, leading to a genome-wide increase in heterozygosity (Table 1) and consequent increase in additive variance (Figure 5).

The observed changes in phenotypic mean and variance are similar to expected outcomes theorized for a finite polygenic architecture with additive allele effects (Chevalet 1994). Consistent with an additive genetic model, several AFPCs showed linear trends across all generations reflective of unconditionally (un)favorable alleles in Hallauer’s Tusón (Figure S5). However, AFPCs with transient shifts in allele frequency were also detected, such as mid-to-late generational responses and plateaus in allele frequency change, highlighting a context-dependent component of the genetic architecture underlying the response to phenotypic selection. The same pattern of plateauing allele frequencies after an initially strong shift was found by temporal analysis of natural populations of Drosophila melanogaster adapted to a novel laboratory environment, which Orozco-terWengel reasoned was due to overdominant or antagonistic pleiotropic effects. It has been demonstrated (mathematically) that the selection coefficient for an additive allele can vary across generations also as a result of changes in background polygenic variance (Chevin and Hospital 2008). The genotypic variance in Hallauer’s Tusón partitioned into additive (primarily) and dominance genetic variance with no apparent epistatic genetic variance, but epistatic genetic effects will contribute to the additive genetic variance component in many cases, such that inferences about gene action should not be drawn from variance components estimates (Hill ; Huang and Mackay 2016). Thus far, genetic studies on flowering time in maize have described an architecture with predominantly additive genetic variance (e.g., Buckler ; Coles ; note that these studies use inbred lines, which precludes estimation of dominance variance), but reports of dominant, overdominant (Coles ) and epsistatic (Blanc ; Durand ) allele effects on variation in flowering characteristics also exist. With a limited sample size for quantitative genetic dissection per generation, our study is unable to clarify the causes or relative contribution of context-dependent effects on the response to phenotypic selection.

Maize is highly diverse (Buckler ), and landraces of maize are locally adapted to a wide range of environments (Committee on The Preservation of Indigenous Strains of Maize 1952–1963). Still, it was surprising that Hallauer’s Tusón captured nearly all of the SNPs on the MaizeSNP50 chip (Ganal ). This high level of molecular genetic variation, as well as the detection of an increasing proportion of polymorphic platform SNPs across generations (Table 1), led us to question whether migrant pollen had entered the population, particularly since it was open-pollinated during selection; although the population was bred in spatial and temporal isolation of other maize populations. Separate lines of evidence indicate the population could have very high diversity while remaining a closed system with no migration or pollen flow. First, similar to our finding, another study has found that individual populations of maize landraces can capture of the SNPs on the same MaizeSNP50 platform (Arteaga ). Because the base population of Hallauer’s Tusón was admixed from multiple, geographically dispersed populations of the landrace Tusón, there is a greater likelihood for the level of diversity to be high. Second, the binomial sampling probabilities for our study limited detection of rare variants within generations despite their putative presence in the population. For instance, based on our sample sizes (which was larger for ), SNPs at a frequency of 1% have an 11% chance of being undetected in and a 33% chance of being undetected in the other generations, but subtle increases above 1% result in large increases in the probability of their detection. Therefore, variants that were not detected in one generation but detected in another may exist at low frequencies, and the transgenerational increase in polymorphic platform SNPs can be explained by selection of initially rare alleles. Finally, considering the most frequent migrant sources in Iowa where the population was selected would be of temperate origin, all of the genotyped individuals in Hallauer’s Tusón clustered with other maize samples of tropical rather than temperate origin (Figure 2).

A fundamental question in genetics is how populations acquire and maintain variation that conditions them with the capacity to adapt to a novel environment. At the locations where the source populations were already adapted and grown, we presume stabilizing selection occurred on flowering time, as flowering time affects fitness in an environmentally dependent manner (Hall and Willis 2006; Mercer and Perales 2019). Stabilizing selection is expected to deplete genetic variation (Barton and Keightley 2002), such that the extensive functional variation for flowering time in Hallauer’s Tusón suggests evolutionary forces beyond mutation and selection for a single optimal flowering time affected the founding populations. Teixeira showed that flowering time variation in Hallauer’s Tusón is under strong genetic and environmental control, with relatively little genotype-by-environment interaction (however, GxE effects were present across latitude and greater in the initial generations). Therefore, seasonal fluctuations that affect the relationship between flowering time and fitness (Giauffret ) and multivariate constraints to evolution (Walsh and Blows 2009) likely contributed to the maintenance of substantial standing variation for this trait, and therefore its capacity for adaptation to a novel environment.

As the population was subjected to directional selection in a temperate environment, alleles contributing to earlier flowering tended to be enriched (Figure 4). Although alleles with negative effects on flowering time spanned the full allele frequency spectrum in , including initially frequent variants such as the photoperiod insensitive allele ZmCCT10-i (present at 75% in ), most of these existed in the minor frequency domain. This was in contrast to minor alleles in at markers, for which had negative effects on flowering time. Hence, in the base generation of Hallauer’s Tusón, favorable alleles for temperate adaptation primarily exist in the minor frequency spectrum.

Due to the admixture of multiple Tusón populations to form , however, the allele frequencies reflect those among (not within) the founder populations, such that native population allele frequencies of temperate-adaptive variants are confounded. To address this, subpopulation assignments for samples from (assumed to correspond to the founding populations) were used to compute subpopulation-specific allele frequencies for the corresponding minor allele in the whole sample (data not shown). Across all markers, most alleles were shared among multiple subpopulations; only of these markers included private alleles. This suggests that geographically separated populations of the landrace Tusón have retained shared alleles that enable latitudinal adaptation, a finding that is congruent with geographical association results from a diverse sample of maize landraces (Romero Navarro ).

Our inference from allele frequency data, however, should be considered with caution. The inferred number of loci (finite polygenic architecture) and effect of selection on variants across the allele frequency spectrum is based on markers that are unlikely to be causal variants themselves, but are expected to tag local haplotype blocks containing causal variants (Nuzhdin and Turner 2013; Kelly and Hughes 2019). Moreover, any initially rare variants with late flowering time effects would likely have been purged or kept at low frequencies, leading to low power of detection for both selection and association mapping, despite the importance of such potential variants for environmental adaptation. This represents a bias to evolutionary inference in experimentally evolved populations and highlights the need for deeper sampling within generations and other approaches to elucidate the structure of natural variation for adaptation.

Our study design allows for the mapping and characterization of specific genomic loci underlying phenotypic evolution. Although the relatively small sample and low marker density (tens of thousands of SNPs) limited the power of our study, we nevertheless detected robust associations with known genes involved in flowering time adaptation, raising confidence in our discovery of other unique loci (Figure S7). For instance, selection against photoperiod sensitivity was a major component of adaptation in Hallauer’s Tusón (Teixeira ). Regions encompassing ZCN8, CONZ1, COL9, CRY2, and ZmCCT9 on chromosomes 8 and 9, and a causal regulatory variant of ZmCCT10 on chromosome 10, were detected by selection or association mapping (Figure 6). These genes are regulators of photoperiodism in plants (Guo ; Cheng and Wang 2005; Miller ; Meng ) that have contributed to latitudinal adaptation (Yang ; Guo ; Huang ), all of which were identified in a large-scale mapping study on photoperiod sensitivity in maize (cf. Figure 6 and Figure 2 in Hung ). Moreover, we found that different local haplotypes were responsive to selection at several of these same loci where Hung detected allelic series [e.g., and 25 (Figure 5) and the ZmCCT10 locus (Figure 7)]. On chromosome 8, no robust association was detected for the Vegetative to Generative Transition 1 (VGT1) gene (Salvi ) involved in earliness per se, but a maize homolog of Arabidopsis thaliana de-etiolated 1 (DET1) involved in photomorphogenesis (Pepper ), not previously highlighted for maize adaptation, mapped to a conspicuous across the centromere of the chromosome.

Development of next-generation crop varieties is crucial to ensuring ample production and quality of plant-based products for society. Reliance on limited pools of diversity for breeding and the creation of monocultures can lead to constrained and vulnerable production systems, but genetic diversity is ecologically structured whereby alleles that could be useful in a target production environment reside in germplasm that suffers from maladaptive syndromes (Teixeira ). Our study shows the potential for maize landraces to be adapted to temperate environments by simple recurrent selection for early flowering time while resulting in minimal loss of molecular genetic variation to reach the adapted state. The lack of a linked footprint of selection encompassing the ZmCCT10 causal site, which had the strongest association with flowering time and underwent a rapid complete-sweep, suggests that some critical adaptive mutations in maize are embedded in regions of low LD, which would permit the maintenance of linked variation during directional selection. Moreover, the multidimensional nature of the genetic architecture underlying response to phenotypic selection, involving multiple loci and alleles with context-dependent effects, appears to enable a rapid shift toward an adapted state with limited loss in diversity. We anticipate that further characterization of these additional layers of the genetic architecture and dynamics of the genomic response to selection will lead to new advances in genomic prediction across generations. By bridging population and quantitative genetic inference, this study advances our understanding of short-term evolution, providing unique insights that aid in developing approaches to adapt crops to climate change.