Trait-by-trait analysis of phenotypic variation within and among populations

In order to investigate phenotypic variation, we set up two common garden experiments located in Mexico to evaluate individuals from 11 teosinte populations (Fig 1). The two experimental fields were chosen because they were located at intermediate altitudes (S1 Fig). Although natural teosinte populations are not typically encountered around these locations [52], we verified that environmental conditions were compatible with both subspecies (S2 Fig). The 11 populations were sampled among 37 populations (S1 Table) distributed along two altitudinal gradients that range from 504 to 2176 m in altitude over ~460 kms for gradient a, and from 342 to 2581m in altitude over ~350 kms for gradient b (S1 Fig). Lowland populations of the subspecies parviglumis (n = 8) and highland populations of the subspecies mexicana (n = 3) were climatically contrasted as can be appreciated in the Principal Component Analysis (PCA) computed on 19 environmental variables (S2 Fig). The corresponding set of individuals grown from seeds sampled from the 11 populations formed the association panel.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Geographical location of sampled populations and experimental fields.

The entire set of 37 Mexican teosinte populations is shown with parviglumis (circles) and mexicana (triangles) populations sampled along gradient a (white) and gradient b (black). The 11 populations indicated with a purple outline constituted the association panel. This panel was evaluated in a four-block design over two years in two experimental fields located at mid-elevation, SENGUA and CEBAJ. Two major cities (Mexico City and Guadalajara) are also indicated. Topographic surfaces have been obtained from International Centre for Tropical Agriculture (Jarvis A., H.I. Reuter, A. Nelson, E. Guevara, 2008, Hole-filled seamless SRTM data V4, International Centre for Tropical Agriculture (CIAT), available from http://srtm.csi.cgiar.org).

https://doi.org/10.1371/journal.pgen.1008512.g001

We gathered phenotypic data during two consecutive years (2013 and 2014). We targeted 18 phenotypic traits that included six traits related to plant architecture, three traits related to leaves, three traits related to reproduction, five traits related to grains, and one trait related to stomata (S2 Table). Each of the four experimental assays (year-field combinations) encompassed four blocks. In each block, we evaluated one offspring (half-sibs) of ~15 mother plants from each of the 11 teosinte populations using a semi-randomized design. After filtering for missing data, the association panel included 1664 teosinte individuals. We found significant effects of Field, Year and/or their interaction for most traits, and a highly significant Population effect for all of them (model 1, S3 Table).

We investigated the influence of altitude on each trait independently. All traits, except for the number of nodes with ears (NoE), exhibited a significant effect of altitude (S3 Table, model 4). Note that after accounting for elevation, the population effect remained significant for all traits, suggesting that factors other than altitude contributed to shape phenotypic variation among populations. Traits related to flowering time and tillering displayed a continuous decrease with elevation, and traits related to grain size increased with elevation (Fig 2 & S3 Fig). Stomata density also diminished with altitude (Fig 2). In contrast, plant height, height of the highest ear, number of nodes with ear in the main tiller displayed maximum values at intermediate altitudes (highland parviglumis and lowland mexicana) (S3 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Population-level box-plots of adjusted means for four traits.

Traits are female flowering time (A), male flowering time (B), grain length (C) and stomata density (D). Populations are ranked by altitude. Parviglumis populations are shown in green and mexicana in red, lighter colors are used for gradient ‘a’ and darker colors for gradient ‘b’. In the case of male and female flowering time, we report data for 9 out of 11 populations because most individuals from the two lowland populations (P1a and P2b) did not flower in our common gardens. Covariation with elevation was significant for the four traits. Corrections for the experimental setting are detailed in the Material and Methods section (Model 2).

https://doi.org/10.1371/journal.pgen.1008512.g002

We estimated narrow-sense heritabilites (additive genotypic effect) per population for all traits using a mixed animal model. Average per-trait heritability ranged from 0.150 for tassel branching to 0.664 for female flowering time, albeit with large standard errors (S2 Table). We obtained higher heritability for grain related traits when mother plant measurements were included in the model with 0.631 (sd = 0.246), 0.511 (sd = 0.043) and 0.274 (sd = 0.160) for grain length, weight and width, respectively, suggesting that heritability was under-estimated for other traits where mother plant values were not available.

Multivariate analysis of phenotypic variation and correlation between traits

Principal component analysis including all phenotypic measurements highlighted that 21.26% of the phenotypic variation scaled along PC1 (Fig 3A), a PC axis that is strongly collinear with altitude (Fig 3B). Although populations partly overlapped along PC1, we observed a consistent tendency for population phenotypic differentiation along altitude irrespective of the gradient (Fig 3C). Traits that correlated the most to PC1 were related to grain characteristics, tillering, flowering and to a lesser extent to stomata density (Fig 3B). PC2 correlated with traits exhibiting a trend toward increase-with-elevation within parviglumis, but decrease-with-elevation within mexicana (Fig 3D). Those traits were mainly related to vegetative growth (Fig 3B). Together, both axes explained 37.4% of the phenotypic variation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Principal component analysis on phenotypic values corrected for the experimental setting.

Individuals factor map (A) and corresponding correlation circle (B) on the first two principal components with altitude (Alt) added as a supplementary variable (in blue). Individual phenotypic values on PC1 (C) and PC2 (D) are plotted against population ranked by altitude and color-coded following A. For populations from the two subspecies, parviglumis (circles) and mexicana (triangles), color intensity indicates ascending elevation in green for parviglumis and red for mexicana. Corrections for experimental setting are detailed in the Material and Methods (Model 3).

https://doi.org/10.1371/journal.pgen.1008512.g003

We assessed more formally pairwise-correlations between traits after correcting for the experimental design and population structure (K = 5). We found 82 (54%) significant correlations among 153 tested pairs of traits. The following pairs of traits had the strongest positive correlations: male and female flowering time, plant height and height of the highest ear, height of the highest and lowest ear, grain length with grain weight and width (S4 Fig). The correlation between flowering time (female or male) with grain weight and length were among the strongest negative correlations (S4 Fig).

Neutral structuring of the association panel

We characterized the genetic structure of the association panel using SSRs. The highest likelihood from Bayesian classification was obtained at K = 2 and K = 5 clusters (S5 Fig). At K = 2, the clustering separated the lowland of gradient a from the rest of the populations. From K = 3 to K = 5, a clear separation between the eight parviglumis and the three mexicana populations emerged. Increasing K values finally split the association panel into the 11 populations it encompassed (S6 Fig). The K = 5 structure reflected both altitude (lowland parviglumis versus highland mexicana) and gradients a and b (Fig 4A & 4B). TreeMix analysis for a subset of 10 of these populations confirmed those results with an early split separating the lowlands from gradient a (cf. K = 2, S6 Fig) followed by the separation of the three mexicana from the remaining populations (Fig 4C). TreeMix results further supported three migration edges, a model that explained 98.75% of the variance and represented a significant improvement over a model without admixture (95.7%, S7 Fig). This admixture model was consistent with gene flow between distant lowland parviglumis populations from gradient a and b, as well as between parviglumis and mexicana populations (Fig 4C). Likewise, Structure analysis also suggested admixture among some of the lowland populations, and to a lesser extent between the two subspecies (Fig 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Genetic clustering, historical splits and admixture among populations of the association panel.

Genetic clustering visualization based on 38 SSRs is shown for K = 5 (A). Colors represent the K clusters. Individuals (vertical lines) are partitioned into colored segments whose length represents the membership proportions to the K clusters. Populations (named after the subspecies M: mexicana, P: parviglumis and gradient ‘a’ or ‘b’) are ranked by altitude indicated in meters above sea level. The corresponding geographic distribution of populations along with their average membership probabilities are plotted (B). Historical splits and admixture between populations were inferred from neutral SNP data for a subset of 10 populations of the association panel (C). Admixtures are colored according to their weight.

https://doi.org/10.1371/journal.pgen.1008512.g004

Identification of traits evolving under spatially-varying selection

We estimated the posterior mean (and 95% credibility interval) of genetic differentiation (F_ST) among the 11 populations of the association panel using driftsel. Considering 1125 plants for which we had both individual phenotypes and individual genotypes for 38 SSRs (S4 Table), we estimated the mean F_ST to 0.22 (0.21–0.23). Note that we found a similar estimate on a subset of 10 of these populations using 1000 neutral SNPs (F_ST (CI) = 0.26 (0.25–0.27)). To identify traits whose variation among populations was driven primarily by local selection, we employed the Bayesian method implemented in driftsel, that infers additive genetic values of traits from a model of population divergence under drift [59]. Selection was inferred when observed phenotypic differentiation exceeded neutral expectations for phenotypic differentiation under random genetic drift. Single-trait analyses revealed evidence for spatially-varying selection at 12 traits, with high consistency between SSRs and neutral SNPs (Table 1). Another method that contrasted genetic and phenotypic differentiation (Q_ST-F_ST) uncovered a large overlap with nine out of the 12 traits significantly deviating from the neutral model (Table 1) and one of the remaining ones displaying borderline significance (Plant height = PL, S8 Fig). Together, these two methods indicated that phenotypic divergence among populations was driven by local selective forces.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Signals of selection (posterior probability S) for each trait considering SSR markers (11 populations) or SNPs (10 populations).

https://doi.org/10.1371/journal.pgen.1008512.t001

Altogether, evidence of spatially varying selection at 10 traits (Table 1) as well as continuous variation of a subset of traits across populations in both elevation gradients (Fig 2, S3 Fig) was consistent with a syndrome where populations produced less tillers, flowered earlier, displayed lower stomata density and carried larger, longer and heavier grains with increasing elevation.

Outlier detection and correlation with environmental variables

We successfully genotyped 218 (~81%) out of 270 outlier SNPs on a broad set of 28 populations, of which 141 were previously detected in candidate regions for local adaptation [58]. Candidate regions were originally identified from re-sequencing data of only six teosinte populations (S1 Table) following an approach that included high differentiation between highlands and lowlands, environmental correlation, and in some cases their intersection with genomic regions involved in quantitative trait variation in maize. The remaining outlier SNPs (77) were discovered in the present study by performing F_ST-scans on the same re-sequencing data (S5 Table). We selected outlier SNPs that were both highly differentiated between highland and lowland populations within gradients (high/low in gradient a or b or both), and between highland and lowland populations within subspecies in gradient b (high/low within parviglumis, mexicana or both). F_ST-scans pinpointed three genomic regions of particularly high differentiation (S9 Fig) that corresponded to previously described inversions [55, 56]: one inversion on chromosome 1 (Inv1n), one on chromosome 4 (Inv4m) and one on the far end of chromosome 9 (Inv9e).

A substantial proportion of outlier SNPs was chosen based on their significant correlation among six populations between variation of allele frequency and their coordinate on the first environmental principal component [58]. We extended environmental analyses to 171outlier SNPs (MAF>5%) on a broader sample of 28 populations (S1 Table) and used the two first components (PCenv1 and PCenv2) to summarize environmental information. When considering all 37 populations, the first component that explained 56% of the variation, correlated with altitude but displayed no correlation to either latitude or longitude. PCenv1 was defined both by temperature- and precipitation- related variables (S2B Fig) including Minimum Temperature of Coldest Month (T6), Mean Temperature of Driest and Coldest Quarter (T9 and T11) and Precipitation of Driest Month and Quarter (P14 and P17). The second PC explained 20.5% of the variation and was mainly defined (S2B Fig) by Isothermality (T3), Temperature Seasonality (T4) and Temperature Annual Range (T7).

We first employed multiple regression to test for each SNP, whether the pairwise F_ST matrix across 28 populations correlated to the environmental (distance along PCenv1) and/or the geographical distance. As expected, we found a significantly greater proportion of environmentally-correlated SNPs among outliers compared with neutral SNPs (χ² = 264.07, P-value = 2.2 10⁻¹⁶), a pattern not seen with geographically-correlated SNPs. That outlier SNPs displayed a greater isolation-by-environment than isolation-by-distance, indicated that patterns of allele frequency differentiation among populations were primarily driven by adaptive processes. We further tested correlations between allele frequencies and environmental variation. Roughly 60.8% (104) of the 171 outlier SNPs associated with at least one of the two first PCenvs, with 87 and 33 associated with PCenv1 and PCenv2, respectively, and little overlap (S5 Table). As expected, the principal component driven by altitude (PCenv1) correlated to allele frequency for a greater fraction of SNPs than the second orthogonal component. Interestingly, we found enrichment of environmentally-associated SNPs within inversions both for PCenv1 (χ² = 14.63, P-value = 1.30 10⁻⁴) and PCenv2 (χ² = 33.77, P-value = 6.22 10⁻⁹).

Associating genotypic variation to phenotypic variation

We tested the association between phenotypes and 171 of the outlier SNPs (MAF>5%) using the association panel. For each SNP-trait combination, the sample size ranged from 264 to 1068, with a median of 1004 individuals (S6 Table). We used SSRs to correct for both structure (at K = 5) and kinship among individual genotypes. This model (model 6) resulted in a uniform distribution of P-values when testing the association between genotypic variation at SSRs and phenotypic trait variation (S10 Fig). Under this model, we found that 126 outlier SNPs (73.7%) associated to at least one trait (Fig 5 and S11 Fig) at an FDR of 10%. The number of associated SNPs per trait varied from 0 for leaf and grain coloration, to 55 SNPs for grain length, with an average of 22.6 SNPs per trait (S5 Table). Ninety-three (73.8%) out of the 126 associated SNPs were common to at least two traits, and the remaining 33 SNPs were associated to a single trait (S5 Table). The ten traits displaying evidence of spatially varying selection in the Q_ST-F_ST analyses displayed more associated SNPs per trait (30.5 on average), than the non-spatially varying traits (12.75 on average).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Manhattan plots of associations between 171 outlier SNPs and 6 phenotypic traits.

X-axis indicates the positions of outlier SNPs on chromosomes 1 to 10, black and gray colors alternating per chromosome. Plotted on the Y-axis are the negative Log₁₀-transformed P values obtained for the K = 5 model. Significant associations (10% FDR) are indicated considering either a structure matrix at K = 5 (pink dots), 11 populations (blue dots) or both K = 5 and 11 populations (purple dots) models.

https://doi.org/10.1371/journal.pgen.1008512.g005

A growing body of literature stresses that incomplete control of population stratification may lead to spurious associations [61]. Hence, highly differentiated traits along environmental gradients are expected to co-vary with any variant whose allele frequency is differentiated along the same gradients, without underlying causal link. We therefore expected false positives in our setting where both phenotypic traits and outlier SNPs varied with altitude. We indeed found a slightly significant correlation (r = 0.5, P-value = 0.03) between the strength of the population effect for each trait—a measure of trait differentiation (S3 Table)–and its number of associated SNPs (S5 Table).

To verify that additional layers of structuring among populations did not cause an excess of associations, we repeated the association analyzes considering a structuring with 11 populations (instead of K = 5) as covariate (model 7), a proxy of the structuring revealed at K = 11 (S6 Fig). With this level of structuring, we retrieved much less associated SNPs (S5 Table). Among the 126 SNPs associating with at least one trait at K = 5, only 22 were recovered considering 11 populations. An additional SNP was detected with structuring at 11 populations that was absent at K = 5. Eight traits displayed no association, and the remaining traits varied from a single associated SNP (Leaf length—LeL and the number of tillers—Til) to 8 associated SNPs for grain weight (S5 Table). For instance, traits such as female or male flowering time that displayed 45 and 43 associated SNPs at K = 5, now displayed only 4 and 3 associated SNPs, respectively (Fig 5). Note that one trait (Leaf color) associated with 4 SNPs considering 11 populations while displaying no association at K = 5. Significant genetic associations were therefore highly contingent on the population structure. Noteworthy, traits under spatially varying selection still associated with more SNPs (2.00 on average) than those with no spatially varying selection (1.25 SNPs on average).

Altogether the 23 SNPs recovered considering a neutral genetic structure with 11 populations corresponded to 30 associations, 7 of the SNPs being associated to more than one trait (S5 Table). For all these 30 associations except in two cases (FFT with SNP_7, and MFT with SNP_28), the SNP effect did not vary among populations (non-significant SNP-by-population interaction in model 7 when we included the SNP interactions with year*field and population). For a subset of two SNPs, we illustrated the regression between the trait value and the shift of allele frequencies with altitude (Fig 6A & 6B). We estimated corresponding additive and dominance effects (S7 Table). In some cases, the intra-population effect corroborated the inter-population variation with relatively large additive effects of the same sign (Fig 6). Note that in both examples shown in Fig 6, one or the other allele was dominant. In other cases, the results were more difficult to interpret with negligible additive effect but extremely strong dominance (S7 Table, SNP_210 for instance).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Regression of phenotypic average value on SNP allele frequency across populations, and within-populationaverage phenotypic value for each SNP genotype.

Per-population phenotypic average values of traits are regressed on allele frequencies at SNP_149 (A) and SNP_179 (B) with corresponding within-population average phenotypic value per genotype (C & D). In A and B, the 11 populations of the association panel are shown with parviglumis (green circles) and mexicana (red triangles) populations sampled along gradient a and gradient b. Phenotypic average values were corrected for the experimental design (calculated as the residues of model 3). Pval refers to the P-value of the linear regression represented in blue. In C and D, genotypic effects from model 7 are expressed as the average phenotypic value of heterozygotes (1) and homozygotes for the alternative allele (2) as compared to the homozygous for the reference allele (0). FDR values were obtained from the association analysis on 171 SNPs with correction for genetic structure using 11 populations.

https://doi.org/10.1371/journal.pgen.1008512.g006

Independence of SNPs associated to phenotypes

We computed the pairwise linkage disequilibrium (LD) as measured by r² between the 171 outlier SNPs using the R package LDcorSV [62]. Because we were specifically interested by LD pattern between phenotypically-associated SNPs, as for the association analyses we accounted for structure (K = 5) and kinship computed from SSRs while estimating LD [63]. The 171 outlier SNPs were distributed along the 10 chromosomes of maize, and exhibited low level of linkage disequilibrium (LD), except for SNPs located on chromosomes eight, nine, and a cluster of SNPs located on chromosome 4 (S12 Fig).

Among the 171, the subset of 23 phenotypically-associated SNPs (detected when considering the 11-populationstructure) displayed an excess of elevated LD values—out of 47 pairs of SNPs phenotypically-associated to a same trait, 16 pairs were contained in the 5% higher values of the LD distribution of all outlier SNP pairs. Twelve out of the 16 pairs associated to grain weight, the remaining four to leaf coloration, and one pair of SNPs associated to both traits. Noteworthy was that inversions on chromosomes 1, 4, and 9, taken together, were enriched for phenotypically-associated SNPs (χ² = 8.95, P-value = 0.0028). We recovered a borderline significant enrichment with the correction K = 5 (χ² = 3.82, P-value = 0.051).

Finally, we asked whether multiple SNPs contributed independently to the phenotypic variation of a single trait. We tested a multiple SNP model where SNPs were added incrementally when significantly associated (FDR < 0.10). We found 2, 3 and 2 SNPs for female, male flowering time and height of the highest ear, respectively (S5 Table). For the two former traits, the SNPs were located on different chromosomes. For the latter trait, the SNPs were both located on chromosome 5 but displayed no LD (SNP_25 and SNP_30, S12 Fig).

The syndrome of altitudinal adaptation results from selection at multiple co-adapted traits

Common garden studies along elevation gradients have been conducted in European and North American plants species [66]. Together with other studies, they have revealed that adaptive responses to altitude are multifarious [67]. They include physiological responses such as high photosynthetic rates [68], tolerance to frost [69], biosynthesis of UV-induced phenolic components [70]; morphological responses with reduced stature [71, 72], modification of leaf surface [73], increase in leaf non-glandular trichomes [74], modification of stomata density; and phenological responses with variation in flowering time [75], and reduced growth period [76].

Our multivariate analysis of teosinte phenotypic variation revealed a marked differentiation between teosinte subspecies along an axis of variation (21.26% of the total variation) that also discriminated populations by altitude (Fig 2A & 2B). The combined effects of assortative mating and environmental elevation variation may generate, in certain conditions, trait differentiation along gradients without underlying divergent selection [77]. While we did not measure flowering time differences among populations in situ, we did find evidence for long distance gene flow between gradients and subspecies (Fig 4A & 4C). In addition, several lines of arguments suggest that the observed clinal patterns result from selection at independent traits and is not solely driven by differences in flowering time among populations. First, two distinct methods accounting for shared population history concur with signals of spatially-varying selection at ten out of the 18 traits (Table 1). Nine of them exhibited a clinal trend of increase/decrease of population phenotypic values with elevation (S3 Fig) within at least one of the two subspecies. This number is actually conservative, because these approaches disregard the impact of selective constraints which in fact tend to decrease inter-population differences in phenotypes. Second, while male and female flowering times were positively correlated, they displayed only subtle correlations (|r|<0.16) with other spatially-varying traits except for grain weight and length (|r| <0.33). Third, we observed convergence at multiple phenotypes between the lowland populations from the two gradients that occurred despite their geographic and genetic distance (Fig 4) again arguing that local adaptation drives the underlying patterns.

Spatially-varying traits that displayed altitudinal trends, collectively defined a teosinte altitudinal syndrome of adaptation characterized by early-flowering, production of few tillers albeit numerous lateral branches, production of heavy, long and large grains, and decrease in stomata density. We also observed increased leaf pigmentation with elevation, although with a less significant signal (S3 Table), consistent with the pronounced difference in sheath color reported between parviglumis and mexicana [78, 79]. Because seeds were collected from wild populations, a potential limitation of our experimental setting is the confusion between genetic and environmental maternal effects. Environmental maternal effects could bias upward our heritability estimates. However, our results corroborate previous findings of reduced number of tillers and increased grain weight in mexicana compared with parviglumis [80]. Thus, although maternal effects could not be fully discarded, we believe they were likely to be weak.

The trend towards depleted stomata density at high altitudes (Fig 2) could arguably represent a physiological adaptation as stomata influence components of plant fitness through their control of transpiration and photosynthetic rate [81]. Indeed, in natural accessions of A. thaliana, stomatal traits showed signatures of local adaptation and were associated with both climatic conditions and water-use efficiency [82]. Furthermore, previous work has shown that in arid and hot highland environments, densely-packed stomata may promote increased leaf cooling in response to desiccation [83] and may also counteract limited photosynthetic rate with decreasing pCO₂ [84]. Accordingly, increased stomata density at high elevation sites has been reported in alpine species such as the European beech [85] as well as in populations of Mimulus guttatus subjected to higher precipitations in the Sierra Nevada [86]. In our case, higher elevations display both arid environment and cooler temperatures during the growing season, features perhaps more comparable to other tropical mountains for which a diversity of patterns in stomatal density variation with altitude has been reported [87]. Further work will be needed to decipher the mechanisms driving the pattern of declining stomata density with altitude in teosintes. Altogether, the altitudinal syndrome was consistent with natural selection for rapid life-cycle shift, with early-flowering in the shorter growing season of the highlands and production of larger propagules than in the lowlands. This altitudinal syndrome evolved in spite of detectable gene flow.

Although we did not formally measure biomass production, the lower number of tillers and higher amount and size of grains in the highlands when compared with the lowlands may reflect trade-offs between allocation to grain production and vegetative growth [88]. Because grains fell at maturity and a single teosinte individual produces hundreds of ears, we were unable to provide a proxy for total grain production. The existence of fitness-related trade-offs therefore still needs to be formally addressed.

Beyond trade-offs, our results more generally question the extent of correlations between traits. In maize, for instance, we know that female and male flowering time are positively correlated and that their genetic control is in part determined by a common set of genes [89]. They themselves further increase with yield-related traits [90]. Response to selection for late-flowering also led to a correlated increase in leaf number in cultivated maize [91], and common genetic loci have been shown to determine these traits as well [92]. Here we found strong positive correlations between traits: male and female flowering time, grain length and width, plant height and height of the lowest or highest ear. Strong negative correlations were observed instead between grain weight and both male and female flowering time. Trait correlations were therefore partly consistent with previous observations in maize, suggesting that they were inherited from wild ancestors [93].

Footprints of past adaptation are relevant to detect variants involved in present phenotypic variation

The overall level of differentiation in our outcrossing system (F_ST ≈22%) fell close to the range of previous estimates (23% [94] and 33% [55] for samples encompassing both teosinte subspecies). It is relatively low compared to other systems such as the selfer Arabidopsis thaliana, where association panels typically display maximum values of F_ST around 60% within 10kb-windows genome-wide [95]. Nevertheless, correction for sample structure is key for statistical associations between genotypes and phenotypes along environmental gradients. This is because outliers that display lowland/highland differentiation co-vary with environmental factors, which themselves may affect traits [96]. Consistently, we found that 73.7% SNPs associated with phenotypic variation at K = 5, but only 13.5% of them did so when considering a genetic structure with 11 populations. Except for one, the latter set of SNPs represented a subset of the former. Because teosinte subspecies differentiation was fully accounted for at K = 5 (as shown by the clear distinction between mexicana populations and the rest of the samples, Fig 4A), the inflation of significant associations at K = 5 is not due to subspecies differentiation, but rather to residual stratification among populations within genetic groups. Likewise, recent studies in humans, where global differentiation is comparatively low [97] have shown that incomplete control for population structure within European samples strongly impacts association results [61, 98]. Controlling for such structure may be even more critical in domesticated plants, where genetic structure is inferred a posteriori from genetic data (rather than a priori from population information) and pedigrees are often not well described. Below, we show that considering more than one correction using minor peaks delivered by the Evanno statistic (S5 Fig) can be informative.

Considering a structure with 5 genetic groups, the number of SNPs associated per trait varied from 1 to 55, with no association for leaf and grain coloration (S5 Table). False positives likely represent a greater proportion of associations at K = 5 as illustrated by a slight excess of small P-values when compared with a correction with 11 populations for most traits (S11 Fig). Nevertheless, our analysis recovered credible candidate adaptive loci that were no longer associated when a finer-grained population structure was included in the model. For instance, at K = 5 we detected Sugary1 (Su1), a gene encoding a starch debranching enzyme that was selected during maize domestication and subsequent breeding [99, 100]. We found that Su1 was associated with variation at six traits (male and female flowering time, tassel branching, height of the highest ear, grain weight and stomata density) pointing to high pleiotropy. A previous study reported association of this gene to oil content in teosintes [101]. In maize, this gene has a demonstrated role in kernel phenotypic differences between maize genetic groups [102]. Su1 is therefore most probably a true-positive. That this gene was no longer recovered with the 11-population structure correction indicated that divergent selection acted among populations. Indeed, allelic frequency was highly contrasted among populations, with most populations fixed for one or the other allele, and a single population with intermediate allelic frequency. With the 11-population correction, very low power is thus left to detect the effect of Su1 on phenotypes.

Although the confounding population structure likely influenced the genetic associations, experimental evidence indicates that an appreciable proportion of the variants recovered with both K = 5 and 11 populations are true-positives (S5 Table). One SNP associated with female and male flowering time, as well as with plant height and grain length (at K = 5 only for the two latter traits) maps within the phytochrome B2 (SNP_210; phyB2) gene. Phytochromes are involved in perceiving light signals and are essential for growth and development in plants. The maize gene phyB2 regulates the photoperiod-dependent floral transition, with mutants producing early flowering phenotypes and reduced plant height [103]. Genes from the phosphatidylethanolamine-binding proteins (PEBPs) family–Zea mays CENTRORADIALIS (ZCN) family in maize—are also well-known to act as promotor and repressor of the floral transition in plants [104]. ZCN8 is the main floral activator of maize [105], and both ZCN8 and ZCN5 strongly associate with flowering time variation [102, 106]. Consistently, we found associations of male and female flowering time with PEBP18 (SNP_15). It is interesting to note that SNPs at two flowering time genes, phyB2 and PEBP18, influenced independently as well as in combination both female and male flowering time variation (S5 Table).

The proportion of genic SNPs associated to phenotypic variation was not significantly higher than that of non-genic SNPs (i.e, SNPs >1kb from a gene) (χ²_{(df = 1)} = 0.043, P-value = 0.84 at K = 5 and χ²_{(df = 1)} = 1.623, P-value = 0.020 with 11 populations) stressing the importance of considering both types of variants [107]. For instance, we discovered a non-genic SNP (SNP_149) that displayed a strong association with leaf width variation as well as a pattern of allele frequency shift with altitude among populations (Fig 6B).

Physically-linked and independent SNPs both contribute to the establishment of adaptive genetic correlations

We found limited LD among our outlier SNPs (S12 Fig) corroborating previous reports (LD decay within <100bp, [58, 94]). However, the subset of phenotypically-associated SNPs displayed greater LD, a pattern likely exacerbated by three Mb-scale inversions located on chromosomes 1 (Inv1n), 4 (Inv4m) and 9 (Inv9e) that, taken together, were enriched for SNPs associated with environmental variables related to altitude and/or SNPs associated with phenotypic variation. Previous work [55, 56] has shown that Inv1n and Inv4m segregate within both parviglumis and mexicana, while two inversions on chromosome 9, Inv9d and Inv9e, are present only in some of the highest mexicana populations; such that all four inversions also follow an altitudinal pattern. Our findings confirmed that three of these inversions possessed an excess of SNPs with high F_ST between subspecies and between low- and high-mexicana populations for Inv9e [57]. Noteworthy Inv9d contains a large ear leaf width quantitative trait locus in maize [107]. Corroborating these results, we found consistent association between the only SNP located within this inversion and leaf width variation in teosinte populations (S5 Table). Overall, our results further strengthen the role of chromosomal inversions in teosinte altitudinal adaptation.

Because inversions suppress recombination between inverted and non-inverted genotypes, their spread has likely contributed to the emergence and maintenance of locally adaptive allelic combinations in the face of gene flow, as reported in a growing number of other models (reviewed in [108]) including insects [109], fish [110], birds [111] and plants [26, 112]. But we also found three cases of multi-SNP determinism of traits (male and female flowering time and height of the highest ear, S5 Table) supporting selection on genetically independent loci. Consistently with Weber et al. [101], we found that individual SNPs account for small proportions of the phenotypic variance (S7 Table). Altogether, these observations are consistent with joint selection of complex traits determined by several alleles of small effects, some of which being maintained in linkage through selection of chromosomal rearrangements.

Conclusion

Elevation gradients provide an exceptional opportunity for investigating variation of functional traits in response to continuous environmental factors at short geographical scales. Here we documented patterns indicating that local adaptation, likely facilitated by the existence of chromosomal inversions, allows teosintes to cope with specific environmental conditions in spite of gene flow. We detected an altitudinal syndrome in teosintes composed of sets of independent traits evolving under spatially-varying selection. Because traits co-varied with environmental differences along gradients, however, statistical associations between genotypes and phenotypes largely depended on control of population stratification. Yet, several of the variants we uncovered seem to underlie adaptive trait variation in teosintes. Adaptive teosinte trait variation may be relevant for maize evolution and breeding. Whether the underlying SNPs detected in teosintes bear similar effects in maize or whether their effects differ in domesticated backgrounds will have to be further investigated.

Ethics statement

All the field work has been done in Mexico in collaboration with Instituto Nacional de InvestigacionesForestales, Agrícolas y Pecuarias, Celaya in Celaya.

Description of teosinte populations and sampling

We used 37 teosinte populations of mexicana (16) and parviglumis (21) subspecies from two previous collections [57, 58, 113] to design our sampling. These populations (S1 Table) are distributed along two altitudinal gradients (Fig 1). We plotted their altitudinal profiles using R ‘raster’ package [114] (S1 Fig). We further obtained 19 environmental variable layers from http://idrisi.uaemex.mx/distribucion/superficies-climaticas-para-mexico. These high-resolution layers comprised monthly values from 1910 to 2009 estimated via interpolation methods [115]. We extracted values of the 19 climatic variables for each population (S1 Table). Note that high throughput sequencing (HTS) data were obtained in a previous study for six populations out of the 37 (M6a, P1a, M7b, P2b, M1b and P8b; Fig 1, S1 Table) to detect candidate genomic regions for local adaptation [58]. The four highest and lowest of these populations were included in the association panel described below.

We defined an association panel of 11 populations on which to perform a genotype-phenotype association study (S1 Table). Our choice was guided by grain availability as well as the coverage of the whole climatic and altitudinal ranges. Hence, we computed Principal Component Analyses (PCA) from environmental variables using the FactoMineR package in R [116] and added altitude to the PCA graphs as a supplementary variable. Our association panel comprised five populations from a first gradient (a)–two mexicana and three parviglumis, and six populations from a second gradient (b)–one mexicana and five parviglumis (Fig 1).

Finally, we extracted available SNP genotypes generated with the MaizeSNP50 Genotyping BeadChipfor 28 populations out of our 37 populations [57] (S1 Table). From this available SNP dataset, we randomly sampled 1000 SNPs found to display no selection footprint [57], hereafter neutral SNPs. Data for neutral SNPs (S1 Data) are available at: https://doi.org/10.6084/m9.figshare.9901472. We used this panel of 28 populations to investigate correlation with environmental variation. Note that 10 out of the 28 populations were common to our association panel, and genotypes were available for 24 to 34 individuals per population, albeit different from the ones of our association mapping panel.

Common garden experiments

We used two common gardens for phenotypic evaluation of the association panel (11 populations). Common gardens were located at INIFAP (Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias) experimental field stations in the state of Guanajuato in Mexico, one in Celaya municipality at the Campo Experimental Bajío (CEBAJ) (20°31’20” N, 100°48’44”W) at 1750 meters of elevation, and one in San Luis de la Paz municipality at the Sitio Experimental Norte de Guanajuato (SENGUA) (21°17’55”N, 100°30’59”W) at 2017 meters of elevation. These locations were selected because they present intermediate altitudes (S1 Fig). The two common gardens were replicated in 2013 and 2014.

The original sampling contained 15 to 22 mother plants per population. Eight to 12 grains per mother plant were sown each year in individual pots. After one month, seedlings were transplanted to the field. Each of the four fields (2 locations, 2 years) was divided into four blocks encompassing 10 rows and 20 columns. We evaluated one offspring of ~15 mother plants from each of the 11 teosinte populations in each block, using a semi-randomized design, i.e. each row containing one or two individuals from each population, and individuals being randomized within row, leading to a total of 2,640 individual teosinte plants evaluated.

SSR genotyping and genetic structuring analyses on the association panel

In order to quantify the population structure and individual kinship in our association panel, we genotyped 46 SSRs (microsatellites, S4 Table). Primers sequences are available from the maize database project [117] and genotyping protocol were previously published [118]. Genotyping was done at the GENTYANE platform (UMR INRA 1095, Clermont-Ferrand, France). Allele calling was performed on electropherograms with GeneMapper Software Applied Biosystems. Allele binning was carried out using Autobin software [119], and further checked manually.

We employed STRUCTURE Bayesian classification software to compute a genetic structure matrix on individual genotypes. Individuals with over 40% missing data were excluded from analysis. We applied the same criterion on SSRs success rate and restricted all analyses to a subset of 38 SSRs (S4 Table). For each number of clusters (K from 2 to 13), we performed 10 independent runs of 500,000 iterations after a burn-in period of 50,000 iterations, and combined these 10 replicates using the LargeKGreedy algorithm from the CLUMPP program [120]. We plotted the resulting clusters using DISTRUCT software. We then used the Evanno method [121] to choose the optimal K value.

We inferred a kinship matrix K from the same SSRs using SPAGeDI [122]. Kinship coefficients were calculated for each pair of individuals as the correlation between allelic states [123]. Since teosintes are outcrossers and expected to exhibit an elevated level of heterozygosity, we estimated intra-individual kinship to fill in the diagonal. We calculated ten kinship matrices, each excluding the SSRs from one out of the 10 chromosomes. Microsatellite data (S2 Data) are available at: https://doi.org/10.6084/m9.figshare.9901472.

In order to gain insights into population history of divergence and admixture, we used 1000 neutral SNPs (i.e. SNPs genotyped by Aguirre-Liguori and collaborators [57] and that displayed patterns consistent with neutrality among 49 teosinte populations) genotyped on 10 out of the 11 populations of the association panel to run a TreeMix analysis (TreeMix version 1.13 [124]). TreeMix models genetic drift to infer populations’ splits from an outgroup as well as migration edges along a bifurcating tree. We oriented the SNPs using the previously published MaizeSNP50 Genotyping BeadChip data from the outgroup species Tripsacumdactyloides [55]. We tested from 0 to 10 migration edges. We fitted both a simple exponential and a non-linear least square model (threshold of 1%) to select the optimal number of migration edges as implemented in the OptM R package [125]. We further verified that the proportion of variance did not substantially increase beyond the optimal selected value.

Phenotypic trait measurements

We evaluated a total of 18 phenotypic traits on the association panel (S2 Table). We measured six traits related to plant architecture (PL: Plant Height, HLE: Height of the Lowest Ear, HHE: Height of the Highest Ear, Til: number of Tillers, LBr: number of Lateral Branches, NoE: number of Nodes with Ears), three traits related to leave morphologies (LeL: Leaf Length, LeW: Leaf Width, LeC: Leaf Color), three traits related to reproduction (MFT: Male Flowering Time, FFT: Female Flowering Time, TBr: Tassel Branching), five traits related to grains (Gr: number of Grains per ear, GrL: Grain Length, GrWi: Grain Width, GrWe: Grain Weight, GrC: Grain Color), and one trait related to Stomata (StD: Stomata Density). These traits were chosen because we suspected they could contribute to differences among teosinte populations based on a previous report of morphological characterization of 112 teosintes grown in five localities [126].

We measured the traits related to plant architecture and leaves after silk emergence. Grain traits were measured at maturity. Leaf and grain coloration were evaluated on a qualitative scale. For stomata density, we sampled three leaves per plant and conserved them in humid paper in plastic bags. Analyses were undertaken at the Institute for Evolution and Biodiversity (University of Münster) as follows: 5mm blade discs were cut out from the mid length of one of the leaves and microscopic images were taken after excitation with a 488nm laser. Nine locations (0.15mm²) per disc were captured with 10 images per location along the z-axis (vertically along the tissue). We automatically filtered images based on quality and estimated leaf stomata density using custom image analysis algorithms implemented in Matlab. For each sample, we calculated the median stomata density over the (up to) nine locations. To verify detection accuracy, manual counts were undertaken for 54 random samples. Automatic and manual counts were highly correlated (R² = 0.82), indicating reliable detection (see S1 Annex Stomata Detection, Dittberner and de Meaux, for a detailed description). The filtered data set of phenotypic measurements (S3 Data) is available at: https://doi.org/10.6084/m9.figshare.9901472.

Statistical analyses of phenotypic variation

In order to test for genetic effects on teosinte phenotypic variation, we decomposed phenotypic values of each trait considering a fixed population effect plus a random mother-plant effect (model 1): (1) where the response variable Y is the observed phenotypic value, μ is the total mean, α_i is the fixed year effect (i = 2013, 2014), β_j the fixed field effect (j = field station: SENGUA, CEBAJ), θ_ij is the year-by-field interaction, γ_k/ij is the fixed block effect (k = 1, 2, 3, 4) nested within the year-by-field combination, δ_l is the fixed effect of the population of origin (l = 1 to 11),χ_il is the year-by-population interaction, ψ_jl is the field-by-population interaction, P_m/l is the random effect of mother plant (m = 1 to 15) nested within population, and ε_ijklm is the individual residue. Identical notations were used in all following models. For the distribution of the effects, the same variance was estimated within all populations. Mixed models were run using ASRemlv.3.0 [127] and MM4LMM v2.0.1 [https://rdrr.io/cran/MM4LMM/man/MM4LMM-package.html, update by F. Laporte] R packages, which both gave very similar results, and fixed effects were tested through Wald tests.

For each trait, we represented variation among populations using box-plots on mean values per mother plant adjusted for the experimental design following model 2: (2) where mother plant within population is considered as a fixed effect. We used the function predict to obtain least-square means (ls-means) of each mother plant, and looked at the tendencies between population’s values. All fixed models were computed using lm package in R, and we visually checked the assumptions of residues independence and normal distribution.

We performed a principal component analysis (PCA) on phenotypic values corrected for the experimental design, using FactoMineR package in R [116] from the residues of model 3 computed using the lm package in R: (3)

Finally, we tested for altitudinal effects on traits by considering the altitude of the sampled population (l) as a covariate (ALT) and its interaction with year and field in model 4: (4) where all terms are equal to those in model 1 except that the fixed effect of the population of origin was replaced by a regression on the population altitude (ALT_l).

Detection of selection acting on phenotypic traits

We aimed at detecting traits evolving under spatially varying selection by comparing phenotypic to neutral genotypic differentiation. Q_ST is a statistic analogous to F_ST but for quantitative traits, which can be described as the proportion of phenotypic variation explained by differences among populations [19, 128]. Significant differences between Q_ST and F_ST can be interpreted as evidence for spatially-varying selection when Q_ST>F_ST [128]. We used the R package Q_STF_STComp[129] that is adequate for experimental designs with randomized half-sibs in outcrossing species. We used individuals that were both genotyped and phenotyped on the association panel to establish the distribution of the difference between statistics (Q_ST-F_ST) under the neutral hypothesis of evolution by drift—using the half-sib dam breeding design and 1000 resamples. We next compared it to the observed difference with 95% threshold cutoff value in order to detect traits under spatially-varying selection.

In addition to Q_ST-F_STanalyses, we employed the driftsel R package [130] to test for signal of selection of traits while accounting for drift-driven population divergence and genetic relatedness among individuals (half-sib design). driftsel is a Bayesian method that compares the probability distribution of predicted and observed mean additive genetic values. It provides the S statistic as output, which measures the posterior probability that the observed population divergence arose under divergent selection (S~1), stabilizing selection (S~0) or genetic drift (intermediate S values) [59]. It is particularly powerful for small datasets, and can distinguish between drift and selection even when Q_ST-F_ST are equal [59]. We first applied rafm to estimate the F_ST value across populations, and the population-by-population coancestry coefficient matrix. We next fitted both the rafm and driftsel models with 15,000 MCMC iterations, discarded the first 5,000 iterations as transient, and thinned the remaining by 10 to provide 1000 samples from the posterior distribution. Note that driftsel was slightly modified because we had information only about the dams, but not the sires, of the phenotyped individuals. We thus modified driftsel with the conservative assumption of all sires being unrelated. Because driftsel does not require that the same individuals were both genotyped and phenotyped, we used SSRs and phenotype data of the association panel as well as the set of neutral SNPs and phenotype data on 10 out of the 11 populations. For the SNP analyses, we selected out of the 1000 neutral SNPs the 465 most informative SNPs based on the following criteria: frequency of the less common variant at least 10%, and proportion of missing data at most 1%. Finally, we estimated from driftsel the posterior probability of the ancestral population mean for each trait as well as deviations of each population from these values.

Both Q_ST-F_ST and driftsel rely on the assumption that the observed phenotypic variation was determined by additive genotypic variation. We thus estimated narrow-sense heritability for each trait in each population to estimate the proportion of additive variance in performance. We calculated per population narrow-sense heritabilites as the ratio of the estimated additive genetic variance over the total phenotypic variance on our common garden measurements using the MCMCglmm R package [131] where half sib family is the single random factor, and the design (block nested within year and field) is corrected as fixed factor. For three grain-related traits, we also ran the same model but including mother plants phenotypic values calculated from the remaining grains not sown. We ran 100,000 iterations with 10,000 burn-in, inverse gamma (0.001; 0.001) as priors. We then calculated the mean and standard deviation of the 11 per population h² estimates.

Pairwise correlations between traits

We evaluated pairwise-correlations between traits by correlating the residues obtained from model 5, that corrects the experiment design (year, field and blocks) as well as the underlying genetic structure estimated from SSRs: (5) where b_n is the slope of the regression of Y on the n^th structure covariate Cⁿ. Structure covariate values (Cⁿ covariates, from STRUCTURE output) were calculated at the individual level, i.e. for each offspring of mother plant m from population l, grown in the year i field j and block k. Cⁿ are thus declared with ijklm indices, although they are purely genetic covariates.

Genotyping of outlier SNPs on 28 populations

We extracted total DNA from each individual plant of the association panel as well as 20 individuals from each of the 18 remaining populations that were not included in the association panel (Table 1). Extractions were performed from 30 mg of lyophilized adult leaf material following recommendations of DNeasy 96 Plant Kit manufacturer (QIAGEN, Valencia, CA, USA). We genotyped outlier SNPs using Kompetitive Allele Specific PCR technology (KASPar, LGC Group) [132]. Data for outlier SNPs (S4 and S5 Data) are available at: https://doi.org/10.6084/m9.figshare.9901472.

Among SNPs identified as potentially involved in local adaptation, 270 were designed for KASPar assays, among which 218 delivered accurate quality data. Of the 218 SNPs, 141 were detected as outliers in two previous studies using a combination of statistical methods—including F_ST-scans [133], Bayescan [32] and Bayenv2 [35, 134], Bayescenv [135]–applied to either six of our teosinte populations [58] or to a broader set of 49 populations genotyped by the Illumina MaizeSNP50 BeadChip [57]. The remaining outlier SNPs (77) were detected by F_ST-scans from six populations (S9 Fig, S5 Table), following a simplified version of the rationale in [58] by considering only differentiation statistics: SNPs were selected if they displayed both a high differentiation (5% highest F_ST values) between highland and lowland populations in at least one of the two gradients, and a high differentiation (5% highest F_ST values) between highland and lowland populations either within parviglumis (P2b and P8b) or within mexicana (M7b and M1b) or both in gradient b (Fig 1). We thereby avoided SNPs fixed between the two subspecies.

Association mapping

We tested the association of phenotypic measurements with outlier SNPs on a subset of individuals for which (1) phenotypic measurements were available, (2) at least 60% of outlier SNPs were adequately genotyped, and (3) kinship and cluster membership values were available from SSR genotyping. For association, we removed SNPs with minor allele frequency lower than 5%.

In order to detect statistical associations between outlier SNPs and phenotypic variation, we used the following mixed model derived from [136]: (6) where ζ is the fixed bi-allelic SNP factor with one level for each of the three genotypes (o = 0, 1, 2; with o = 1 for heterozygous individuals), and u_ijklm is the random genetic effect of the individual. We assumed that the vector of u_ijklm effects followed a N(0,K σ²u) distribution, where K is the kinship matrix computed as described above.

A variant of model 6 was employed to test for SNP association to traits, while correcting for structure as the effect of population membership (δ_l), δ being a factor with 11 levels (populations): (7)

In order to avoid overcorrection of neutral genetic structure and improve power, we ran the two models independently for each chromosome using a kinship matrix K estimated from all SSRs except those contained in the chromosome of the tested SNP [137]. We tested SNP effects through the Wald statistics, and applied a 10% False Discovery Rate (FDR) threshold for each phenotype separately. In order to validate the correction for genetic structure, the 38 multiallelic SSR genotypes were transformed into biallelic genotypes, filtered for MAF > 5%, and used to run associations with the complete 6 and 7 models, as well as 6-type models excluding either kinship or both structure and kinship. For each trait, we generated QQplots of P-values for each of these models.

Multiple SNP models were built by successively adding at each step the most significant SNP, as long as its FDR was lower than 0.10. We controlled for population structure considering 11 populations and used the kinship matrix that excluded the SSRs on the same chromosome as the last tested SNP.

Environmental correlation of outlier SNPs

We tested associations between allelic frequency at 171 outlier SNPs and environmental variables across 28 populations, using Bayenv 2.0 [40]. Because environmental variables are highly correlated, we used the first two principal component axes from the environmental PCA analysis (PCenv1 and PCenv2) to run Bayenv 2.0. This software requires a neutral covariance matrix, that we computed from the available dataset of 1000 neutral SNPs (S1 Table). We performed 100,000 iterations, saving the matrix every 500 iterations. We then tested the correlation of these to the last matrix obtained, as well as to an F_ST matrix calculated with BEDASSLE [138], as described in [57].

For each outlier SNP, we compared the posterior probability of a model that included an environmental factor (PCenv1 or PCenv2) to a null model. We determined a 5% threshold for significance of environmental association by running 100,000 iterations on neutral SNPs. We carried out five independent runs for each outlier SNP and evaluated their consistency from the coefficient of variation of the Bayes factors calculated among runs.

In order to test whether environmental distance was a better predictor of allele frequencies at candidate SNPs than geography, we used multiple regression on distance matrices (MRM, [139]) implemented in the ecodist R package [140] for each outlier SNP. We used pairwise F_ST values as the response distance matrix and the geographic and environmental distance matrices as explanatory matrices. We evaluated the significance of regression coefficients by 1000 permutations and iterations of the MRM. We determined the total number of environmentally and geographically associated SNPs (P-value<0.05) among outliers. We employed the same methodology for our set of 1000 neutral SNPs.