Margaux-Alison Fustier

Natalia Elena

Martinez-Ainsworth Jonas

Andres Aguirre-Liguori

Juliette De Meaux

E Martínez-Ainsworth

J A Aguirre-Liguori

Anthony Venon

Hélène Corti

N E Martínez-Ainsworth

A Rousselet

F Dumas

H Dittberner

M G Camarena

D Grimanelli

O Ovaskainen

M Falque

L Moreau

S Montes

L E Eguiarte

Y Vigouroux

D Manicacci
email: domenica.manicacci@inra.fr

M I Tenaillon
email: maud.tenaillon@inra.fr

Common gardens in teosintes reveal the establishment of a syndrome of adaptation to altitude

Keywords: spatially-varying selection, F ST -scan, association mapping, altitudinal syndrome, pleiotropy, chromosomal inversions

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Local adaptation is key for the preservation of ecologically useful genetic variation [START_REF] Whitlock | Modern approaches to local adaptation[END_REF]. The conditions for its emergence and maintenance have been the focus of a long-standing debate nourished by ample theoretical work [START_REF] Bradshaw | Ecological significance of genetic variation between populations[END_REF][START_REF] Bulmer | Multiple niche polymorphism[END_REF][START_REF] Endler | Natural selection in the wild[END_REF][START_REF] Gay | Comparing clines on molecular and phenotypic traits in hybrid zones: A window on tension zone models[END_REF][START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF][START_REF] Lenormand | Gene flow and the limits to natural selection[END_REF][START_REF] Whitlock | Probability of fixation in a heterogeneous environment[END_REF][START_REF] Yeaman | Establishment and maintenance of adaptive genetic divergence under migration, selection and drift[END_REF]. On the one hand, spatially-varying selection promotes the evolution of local adaptation, provided that there is genetic diversity underlying the variance of fitness-related traits [START_REF] Rundle | Ecological speciation[END_REF]. On the other hand, opposing forces such as neutral genetic drift, temporal fluctuations of natural selection, recurrent introduction of maladaptive alleles via migration and homogenizing gene flow may hamper local adaptation (reviewed in [START_REF] Kawecki | Conceptual issues in local adaptation[END_REF]) . Metaanalyzes indicate that local adaptation is pervasive in plants, with evidence of nativesite fitness advantage in reciprocal transplants detected in 45% to 71% of populations [START_REF] Hereford | A quantitative survey of local adaptation and fitness trade-offs[END_REF][START_REF] Leimu | A meta-analysis of local adaptation in plants[END_REF].

While local adaptation is widespread, much has yet to be discovered about the traits affected by spatially-varying selection, their molecular determinants and the underlying ecological drivers [START_REF] Tiffin | Advances and limits of using population genetics to understand local adaptation[END_REF]. Local adaptation is predicted to increase with phenotypic, genotypic and environmental divergence among populations [START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF][START_REF] Garcia-Ramos | Genetic models of adaptation and gene flow in peripherial populations[END_REF][START_REF] Slatkin | Rare alleles as indicators of gene flow[END_REF].

Comparisons of the quantitative genetic divergence of a trait (Q ST ) with the neutral genetic differentiation (F ST ) can provide hints on whether trait divergence is driven by spatially-divergent selection [START_REF] Lande | Neutral theory of quantitative genetic variance in an island model with local extinction and colonization[END_REF][START_REF] Whitlock | Neutral additive genetic variance in a metapopulation[END_REF][START_REF] Spitze | Population structure in Daphina obtusa: quantitative genetic and allozymic variation[END_REF][START_REF] Wright | The genetical structure of populations[END_REF]. Striking examples of divergent selection include developmental rate in the common toad [START_REF] Luquet | Phenotypic divergence of the common toad (Bufo bufo) along an altitudinal gradient: evidence for local adaptation[END_REF], drought and frost tolerance in alpine populations of the European silver fir [START_REF] Roschanski | Evidence of divergent selection for drought and cold tolerance at landscape and local scales in Abies alba Mill. in the French Mediterranean Alps[END_REF], and traits related to plant phenology, size and floral display among populations of Helianthus species [START_REF] Kawakami | Natural selection drives clinal life history patterns in the perennial sunflower species, Helianthus maximiliani[END_REF][START_REF] Moyers | Remarkable life history polymorphism may be evolving under divergent selection in the silverleaf sunflower[END_REF]. These studies have reported covariation of physiological, morphological and/or life-history traits across environmental gradients which collectively define adaptive syndromes. Such syndromes may result from several non-exclusive mechanisms: plastic responses, pleiotropy, non-adaptive genetic correlations among traits (constraints), and joint selection of traits encoded by different sets of genes resulting in adaptive correlations.

In some cases, the latter mechanism may involve selection and rapid spread of chromosomal inversions that happen to capture multiple locally favored alleles [START_REF] Kirkpatrick | Chromosome inversions, local adaptation and speciation[END_REF] as exemplified in Mimulus guttatus [START_REF] Lowry | A widespread chromosomal inversion polymorphism contributes to a major life-history transition, local adaptation, and reproductive isolation[END_REF]. While distinction between these mechanisms is key to decipher the evolvability of traits, empirical data on the genetic bases of correlated traits are currently lacking [START_REF] Legrand | Evolution of a butterfly dispersal syndrome[END_REF].

The genes mediating local adaptation are usually revealed by genomic regions harboring population-specific signatures of selection. These signatures include alleles displaying greater-than-expected differentiation among populations [START_REF] Bierne | The coupling hypothesis: Why genome scans may fail to map local adaptation genes[END_REF] and can be identified through F ST -scans [START_REF] Lewontin | Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms[END_REF][START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF][START_REF] Vitalis | Interpretation of variation across marker loci as evidence of selection[END_REF][START_REF] Foll | A genome scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective[END_REF][START_REF] Lah | Detecting loci under selection in a hierarchically structured population[END_REF][START_REF] Bonhomme | Detecting selection in population trees: The Lewontin and Krakauer test extended[END_REF][START_REF] Günther | Robust identification of local adaptation from allele frequencies[END_REF]. However, F ST -scans and its derivative methods [START_REF] Bierne | The coupling hypothesis: Why genome scans may fail to map local adaptation genes[END_REF] suffer from a number of limitations, among them a high number of false positives (reviewed in [START_REF] Lotterhos | Evaluation of demographic history and neutral parameterization on the performance of FST outlier tests[END_REF][START_REF] Haasl | Fifteen years of genomewide scans for selection: trends, lessons and unaddressed genetic sources of complication[END_REF]) and the lack of power to detect true positives [START_REF] Corre | The genetic differentiation at quantitative trait loci under local adaptation[END_REF]. Despite these caveats, F ST -outlier approaches have helped in the discovery of emblematic adaptive alleles such as those segregating at the EPAS1 locus in Tibetan human populations adapted to high altitude [START_REF] Yi | Sequencing of fifty human exomes reveals adaptations to high altitude[END_REF]. An alternative to detect locally adaptive loci is to test for genotype-environment correlations [START_REF] Günther | Robust identification of local adaptation from allele frequencies[END_REF][START_REF] Coop | Using environmental correlations to identify loci underlying local adaptation[END_REF][START_REF] Guillot | A unifying model for the analysis of phenotypic, genetic, and geographic data[END_REF][START_REF] Frichot | Testing for associations between loci and environmental gradients using latent factor mixed models[END_REF][START_REF] Gautier | Genome-wide scan for adaptive divergence and association with population-specific covariates[END_REF][START_REF] Joost | A spatial analysis method (SAM) to detect candidate loci for selection: Towards a landscape genomics approach to adaptation[END_REF][START_REF] Poncet | Tracking genes of ecological relevance using a genome scan in two independent regional population samples of Arabis alpina[END_REF]. Correlation-based methods can be more powerful than differentiation-based methods [START_REF] De Mita | Detecting selection along environmental gradients: Analysis of eight methods and their effectiveness for outbreeding and selfing populations[END_REF], but spatial autocorrelation of population structure and environmental variables can lead to spurious signatures of selection [START_REF] Hoban | Finding the genomic basis of local adaptation: pitfalls, practical solutions, and future directions[END_REF].

Ultimately, to identify the outlier loci that have truly contributed to improve local fitness, a link between outliers and phenotypic variation needs to be established.

The most common approach is to undertake association mapping. However, recent literature in humans has questioned our ability to control for sample stratification in such approach [START_REF] Barton | Why structure matters[END_REF]. Detecting polymorphisms responsible for trait variation is particularly challenging when trait variation and demographic history follow parallel environmental (geographic) clines. Plants however benefit from the possibility of conducting replicated phenotypic measurements in common gardens, where environmental variation is controlled. Hence association mapping has been successfully employed in the model plant species Arabidopsis thaliana, where broadly distributed ecotypes evaluated in replicated common gardens have shown that fitness-associated alleles display geographic and climatic patterns indicative of selection [START_REF] Fournier-Level | A map of local adaptation in Arabidopsis thaliana[END_REF]. Furthermore, the relative fitness of A. thaliana ecotypes in a given environment could be predicted from climate-associated SNPs [START_REF] Hancock | Adaptation to climate across the Arabidopsis thaliana genome[END_REF]. While climatic selection over broad latitudinal scales produces genomic and phenotypic patterns of local adaptation in the selfer plant A. thaliana, whether similar patterns exist at shorter spatial scale in outcrossing species remains to be elucidated.

We focused here on a well-established outcrossing plant system, the teosintes, to investigate the relationship of molecular, environmental, and phenotypic variation in populations sampled across two elevation gradients in Mexico. The gradients covered a relatively short yet climatically diverse, spatial scale. They encompassed populations of two teosinte subspecies that are the closest wild relatives of maize, Zea mays ssp. parviglumis (hereafter parviglumis) and Z. mays ssp. mexicana (hereafter mexicana). The two subspecies display large effective population sizes [START_REF] Ross-Ibarra | Historical divergence and gene flow in the genus Zea[END_REF], and span a diversity of climatic conditions, from warm and mesic conditions below 1800 m for parviglumis, to drier and cooler conditions up to 3000 m for mexicana [START_REF] Hufford | Inferences from the historical distribution of wild and domesticated maize provide ecological and evolutionary insight[END_REF]. Previous studies have discovered potential determinants of local adaptation in these systems.

At a genome-wide scale, decrease in genome size correlates with increasing altitude, which likely results from the action of natural selection on life cycle duration [START_REF] Bilinski | Parallel altitudinal clines reveal trends in adaptive evolution of genome size in Zea mays[END_REF][START_REF] Diez | Genome size variation in wild and cultivated maize along altitudinal gradients[END_REF].

More modest structural changes include megabase-scale inversions that harbor clusters of SNPs whose frequencies are associated with environmental variation [START_REF] Pyhäjärvi | Complex patterns of local adaptation in teosinte[END_REF][START_REF] Fang | Megabase-scale inversion polymorphism in the wild ancestor of maize[END_REF]. Also, differentiation-and correlation-based genome scans in teosinte populations succeeded in finding outlier SNPs potentially involved in local adaptation [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF][START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF].

But a link with phenotypic variation has yet to be established.

In this paper, we genotyped a subset of these outlier SNPs on a broad sample of 28 teosinte populations, for which a set of neutral SNPs was also available; as well as on an association panel encompassing 11 populations. We set up common gardens in two locations to evaluate the association panel for 18 phenotypic traits over two consecutive years. Individuals from this association panel were also genotyped at 38 microsatellite markers to enable associating genotypic to phenotypic variation while controlling for sample structure and kinship among individuals. We addressed three main questions: What is the extent of phenotypic variation within and among populations? Can we define a set of locally-selected traits that constitute a syndrome of adaptation to altitude? What are the genetic bases of such syndrome? We further discuss the challenges of detecting phenotypically-associated SNPs when trait and genetic differentiation parallel environmental clines.

Results

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

In order to investigate phenotypic variation, we set up two common garden 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 M1, 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, M3 model). 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. 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. 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 M2).

Multivariate analysis of phenotypic variation and correlation between traits.

Principal component analysis including all phenotypic measurements

We assessed more formally pairwise-correlations between traits after correcting for experimental design and population structure (K=5). We found 82

Neutral structuring of the association panel.

We characterized the genetic structure of the association panel using SSRs.

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 with 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 [START_REF] Ovaskainen | A new method to uncover signatures of divergent and stabilizing selection in quantitative traits[END_REF]. 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. and marked by an asterisk when significant in Q ST -F ST Comp analysis. We considered the underlined traits as spatially varying. For a detailed description of traits see S2 Table : Values reported correspond to S from DRIFTSEL. S is the posterior probability that divergence among populations was not driven by drift only. Following [START_REF] Mckinney | Genetic and morphological divergence in three strains of brook trout Salvelinus fontinalis commonly stocked in Lake Superior[END_REF], we used here a conservative credibility value of S > 0.95 to declare divergent selection.

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

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 [START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF]. 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 [START_REF] Soularue | Evolutionary responses of tree phenology to the combined effects of assortative mating, gene flow and divergent selection[END_REF] 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 gradient (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 [START_REF] Pyhäjärvi | Complex patterns of local adaptation in teosinte[END_REF][START_REF] Fang | Megabase-scale inversion polymorphism in the wild ancestor of maize[END_REF]: 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 [START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF]. 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. and PCenv2 (χ² = 33.77, P-value=6. 22 10 -9 ).

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 ). 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).

A growing body of literature stresses that incomplete control of population stratification may lead to spurious associations [START_REF] Sohail | Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies[END_REF]. 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 ). 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 M5' 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 6 A&B). We estimated corresponding additive and dominance effects (S7 Table ). In some cases, the intra-population effect corroborated the interpopulation variation with relatively large additive effects of the same sign (

Independence of SNPs associated to phenotypes.

We computed the pairwise linkage disequilibrium (LD) as measured by r 2 between the 171 outlier SNPs using the R package LDcorSV [START_REF] Desrousseaux | Package ' LDcorSV[END_REF]. Because we were specifically interested by LD pattern between phenotypically-associated SNPs, as for the association analyses we accounted for structure and kinship computed from SSRs while estimating LD [START_REF] Mangin | Novel measures of linkage disequilibrium that correct the bias due to population structure and relatedness[END_REF]. 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-population structure) 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 related to grain weight, the remaining four to leaf coloration, and one pair of SNPs was 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 ). Except for the latter trait, the SNPs were located on different chromosomes.

Discussion

Plants are excellent systems to study local adaptation. First, owing to their sessile nature, local adaptation of plant populations is pervasive [START_REF] Leimu | A meta-analysis of local adaptation in plants[END_REF]. Second, environmental effects can be efficiently controlled in common garden experiments, facilitating the identification of the physiological, morphological and phenological traits influenced by spatially-variable selection [START_REF] Savolainen | Ecological genomics of local adaptation[END_REF]. Identification of the determinants of complex trait variation and their covariation in natural populations is however challenging [START_REF] Anderson | Evolutionary genetics of plant adaptation[END_REF]. While population genomics has brought a flurry of tools to detect footprints of local adaptation, their reliability remains questioned [START_REF] Sohail | Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies[END_REF]. In addition, local adaptation and demographic history frequently follow the same geographic route, making the disentangling of trait, molecular, and environmental variation, particularly arduous. Here we investigated those links on a well-established outcrossing system, the closest wild relatives of maize, along altitudinal gradients that display considerable environmental shifts over short geographical scales.

The syndrome of altitudinal adaptation results from selection at multiple coadapted traits.

Common garden studies along elevation gradients have been conducted in European and North American plants species [START_REF] Halbritter | Trait differentiation and adaptation of plants along elevation gradients[END_REF]. Together with other studies, they have revealed that adaptive responses to altitude are multifarious [START_REF] Körner | The use of 'altitude' in ecological research[END_REF]. They include physiological responses such as high photosynthetic rates [START_REF] Friend | Field measurements of photosynthesis, stomatal conductance, leaf nitrogen and δ 13 C along altitudinal gradients in Scotland[END_REF], tolerance to frost [START_REF] Neuner | Frost resistance in alpine woody plants[END_REF],

biosynthesis of UV-induced phenolic components [70]; morphological responses with reduced stature [START_REF] Byars | Local adaptation and cogradient selection in the alpine plant, Poa hiemata, along a narrow altitudinal gradient[END_REF][START_REF] Luo | The roles of genetic drift and natural selection in quantitative trait divergence along an altitudinal gradient in Arabidopsis thaliana[END_REF], modification of leaf surface [START_REF] Guerin | Leaf morphology shift linked to climate change[END_REF], increase in leaf nonglandular trichomes [START_REF] Kofidis | Contemporary seasonal and altitudinal variations of leaf structural features in oregano (Origanum vulgare L.)[END_REF], modification of stomata density; and phenological responses with variation in flowering time [START_REF] Mendez-Vigo | Altitudinal and climatic adaptation is mediated by flowering traits and FRI, FLC, and PHYC genes in Arabidopsis[END_REF], and reduced growth period [START_REF] Oleksyn | Growth and physiology of Picea abies populations from elevational transects: common garden evidence for altitudinal ecotypes and cold adaptation[END_REF].

Our multivariate analysis of teosinte phenotypic variation revealed a marked differentiation between teosinte subspecies along an axis of variation ( 21 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 [START_REF] Doebley | Maize introgression into teosinte --a reappraisal[END_REF][START_REF] Lauter | The inheritance and evolution of leaf pigmentation and pubescence in teosinte[END_REF]. 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 [START_REF] Smith | Variation within teosinte. I. Numerical analysis of morphological data[END_REF]. 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 (S3 Fig) could

arguably represent a physiological adaptation as stomata influence components of plant fitness through their control of transpiration and photosynthetic rate [START_REF] Raven | Selection pressures on stomatal evolution[END_REF]. 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 [START_REF] Dittberner | Natural variation in stomata size contributes to the local adaptation of water-use efficiency in Arabidopsis thaliana[END_REF]. Furthermore, previous work has shown that in arid and hot highland environments, densely-packed stomata may promote increased leaf cooling in response to desiccation [START_REF] Carlson | Intraspecific variation in stomatal traits, leaf traits and physiology reflects adaptation along aridity gradients in a South African shrub[END_REF] and may also counteract limited photosynthetic rate with decreasing pCO 2 [START_REF] Körner | Stomatal behaviour in alpine plant communities between 600 and 2600 metres above sea level[END_REF]. Accordingly, increased stomata density at high elevation sites has been reported in alpine species such as the European beech [START_REF] Bresson | To what extent is altitudinal variation of functional traits driven by genetic adaptation in European oak and beech?[END_REF] as well as in populations of Mimulus guttatus subjected to higher precipitations in the Sierra Nevada [START_REF] Kooyers | Replicate altitudinal clines reveal that evolutionary flexibility underlies adaptation to drought stress in annual Mimulus guttatus[END_REF]. 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 [START_REF] Körner | Functional morphology of mountain plants[END_REF]. 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 [START_REF] Jakobsson | A comparative study of seed number, seed size, seedling size and recruitment in grassland plants[END_REF]. 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 [START_REF] Buckler | The genetic architecture of maize flowering time[END_REF]. They themselves further increase with yield-related traits [START_REF] Moreau | Use of trial clustering to study QTL x environment effects for grain yield and related traits in maize[END_REF]. Response to selection for late-flowering also led to a correlated increase in leaf number in cultivated maize [START_REF] Durand | Flowering time in maize: Linkage and epistasis at a major effect locus[END_REF], and common genetic loci have been shown to determine these traits as well [START_REF] Li | The genetic architecture of leaf number and its genetic relationship to flowering time in maize[END_REF]. 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.

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 within the range of previous estimates (23% [START_REF] Aguirre-Liguori | Divergence with gene flow is driven by local adaptation to temperature and soil phosphorus concentration in teosinte subspecies ( Zea mays parviglumis and Zea mays mexicana)[END_REF]and 33% [START_REF] Pyhäjärvi | Complex patterns of local adaptation in teosinte[END_REF] 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 [START_REF] Consortium | 1,135 Genomes reveal the global pattern of polymorphism in Arabidopsis thaliana[END_REF].

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 [START_REF] Novembre | Tread lightly interpreting polygenic tests of selection[END_REF]. 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 [START_REF] Guo | Global genetic differentiation of complex traits shaped by natural selection in humans[END_REF] have shown that incomplete control for population structure within European samples strongly impacts association results [START_REF] Sohail | Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies[END_REF][START_REF] Berg | Reduced signal for polygenic adaptation of height in UK Biobank[END_REF]. 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 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 ).

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) (χ² Physically-linked and independent SNPs both contribute to the establishment of adaptive genetic correlations.

We found limited LD among our outlier SNPs (S10 Fig) corroborating previous reports (LD decay within <100bp, [START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF][START_REF] Aguirre-Liguori | Divergence with gene flow is driven by local adaptation to temperature and soil phosphorus concentration in teosinte subspecies ( Zea mays parviglumis and Zea mays mexicana)[END_REF]). 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 [START_REF] Pyhäjärvi | Complex patterns of local adaptation in teosinte[END_REF][START_REF] Fang | Megabase-scale inversion polymorphism in the wild ancestor of maize[END_REF] 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 highmexicana populations for Inv9e [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF]. Noteworthy Inv9d contains a large ear leaf width quantitative trait locus in maize [106]. 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 [107]) including insects [108], fish

[109], birds [110] and plants [START_REF] Lowry | A widespread chromosomal inversion polymorphism contributes to a major life-history transition, local adaptation, and reproductive isolation[END_REF]111]. But we also found three cases of multi-SNP determinism of traits (male and female flowering time and height of the highest ear, Table S5) supporting selection of genetically independent loci. Consistently with Weber et al.

[100], 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 is likely 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.

Finally, we extracted available SNP genotypes generated with the MaizeSNP50 Genotyping BeadChip for 28 populations out of our 37 populations [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF] (S1 Table ). From this available SNP dataset, we randomly sampled 1000 SNPs found to display no selection footprint [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF], hereafter neutral SNPs. Data for neutral SNPs (Data S1) are available at: 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

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 (S4 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. 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 [118]. We plotted the resulting clusters using DISTRUCT software. We then used the Evanno method [119] to choose the optimal K value. We followed the same methodology to compute a structure matrix from the outlier SNPs.

We inferred a kinship matrix K from the same SSRs using SPAGeDI 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 [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF] 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 [122]. 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 Tripsacum dactyloides [START_REF] Pyhäjärvi | Complex patterns of local adaptation in teosinte[END_REF]. 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 [123]. 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 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 followed: 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 2 ) 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 StomataDetection, Dittberner and de Meaux, for a detailed description). The filtered data set of phenotypic measurements (Data S3) is available at: 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 M1):

ܻ ൌ ߤ ߙ ߚ ߠ ߛ / ߝ ሺM2ሻ
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 M3:

ܻ ൌ ߤ ߙ ߚ ߠ ߛ / ܿ. ܶܮܣ ܽ . ܶܮܣ ܾ . ܶܮܣ ܲ / ߝ ሺM3ሻ
where all terms are equal to those in model M1 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 [START_REF] Spitze | Population structure in Daphina obtusa: quantitative genetic and allozymic variation[END_REF]107].

Significant differences between Q ST and F ST can be interpreted as evidence for spatially-varying (Q ST >F ST ) selection [126]. We used the R package QstFstComp

[127] 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 halfsib 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 spatiallyvarying selection.

In addition to Q ST -F ST analyses, we employed the DRIFTSEL R package [128] 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) [START_REF] Ovaskainen | A new method to uncover signatures of divergent and stabilizing selection in quantitative traits[END_REF]. It is particularly powerful for small datasets, and can distinguish between drift and selection even when Q ST -F ST are equal [START_REF] Ovaskainen | A new method to uncover signatures of divergent and stabilizing selection in quantitative traits[END_REF]. We first applied RAFM to estimate the 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 narrowsense 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 [129] 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 M4, that corrects the experiment design (year, field and blocks)

as well as the underlying genetic structure estimated from SSRs:

ܻ ൌ ߤ ߙ ߚ ߠ ߛ / ܾ . ܥ ସ ୀଵ ߝ ሺM4ሻ

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 associations with the complete M5 and M5' models, as well as the M5 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 SSR 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 [START_REF] Coop | Using environmental correlations to identify loci underlying local adaptation[END_REF]111]. 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 [135], as described in [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF].

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.

female flowering time, we report values for all 11 populations although very few individuals from the two most lowland populations (P1a and P2b) flowered.

Covariation with altitude was significant for all traits except for the number of nodes with ears on the main tiller (S3 Table ).

  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.

Figure 1 .

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

  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).

Figure 2 :

 2 Figure 2: Population-level box-plots of adjusted means for four traits. Traits are

highlighted that 21 .

 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 decreasewith-elevation within mexicana (Fig 3D). Those traits were mainly related to vegetative growth (Fig 3B). Together, both axes explained 37% of the phenotypic variation.

Figure 3 :

 3 Figure 3: Principal Component Analysis on phenotypic values corrected for the

(

  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).

  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 associated to both altitude (lowland parviglumis versus highland mexicana) and gradients a and b (Fig 4A & B). 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 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%, Figure S7). 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).

Figure 4 :

 4 Figure 4: Genetic clustering, historical splits and admixture among populations

  .

  b

gradients (Fig 2 ,

 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.

  We extended environmental analyses to 171 outlier SNPs (MAF>5%) on a broader sample of 28 populations (S2 Fig) and used the two first components (PCenv1 and PCenv2) to summarize environmental information. 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 (S2 B 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 (S2 B Fig) by Isothermality (T3), Temperature Seasonality (T4) and Temperature Annual Range (T7).

  Fig).Under this model, we found that 126 outlier SNPs (73.7%) associated to at least

Figure 5 :

 5 Figure 5: Manhattan plots of associations between 171 outlier SNPs and 6

  Fig 6 C&D). 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 (S7Table, SNP_210 for instance).

Figure 6 :

 6 Figure 6: Regression of phenotypic average value on SNP allele frequency across

  Fig) within at least one of the two subspecies. This number is actually conservative,

  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

  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 [102]. 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 [103]. ZCN8 is the main floral activator of maize [104], and both ZCN8 and ZCN5 strongly associate with flowering time variation [101, 105]. 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

  (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 [106]. 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).

  populations). Common gardens were located at INIFAP (Instituto Nacional de Investigaciones Forestales, Agricolas y Pecuaria) 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 in the field. Each of the four fields (2 locations, 2 years) was separated 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.

  [120]. Kinship coefficients were calculated for each pair of individuals as correlation between allelic states [121]. 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 (Data S2) are available at: 10.6084/m9.figshare.9901472

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 a 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.

Figure S4 :

 S4 Figure S4: Pairwise correlations between phenotypic traits. Pearson coefficient

Figure S5 .

 S5 Figure S5. Evanno method calculations for population number ∆K in the association panel genotyped for 38 SSRs.

Figure S6 .

 S6 Figure S6. Genetic clustering of ancestry proportions in the association panel

Figure S7 .

 S7 Figure S7. Determination of the migration edge number in the TreeMix model.

Figure S8 :

 S8 Figure S8: Significance of Q ST -F ST difference for each trait. The dotted blue line

Figure S9 : 9 (

 S99 Figure S9: Genomic F ST -scans on 6 teosinte populations. We computed 4 pairwise-

Figure S10 :

 S10 Figure S10: QQ-plots of observed P-values and expected P-values generated

Figure S11 :

 S11 Figure S11: Manhattan plots of associations between 171 outlier SNPs and 12

Figure S12 :

 S12 Figure S12: Pairwise Linkage Disequilibrium (LD) between outlier SNPs.Pairwise LD between 171 SNPs was estimated using r 2 , and corrected for structure at

Table 1 .

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

	Traits a	SSR b	SNP b
	Plant height	0.995	0.972
	Height of the lowest ear*	0.950	0.959
	Height of the highest ear	0.982	0.966
	Number of tillers*	1.000	1.000
	Number of lateral branches*	1.000	0.990
	Number of nodes with ears	0.682	0.699
	Leaf length	0.888	0.875
	Leaf width	0.999	0.996
	Leaf color	0.633	0.583
	Female flowering time*	1.000	1.000
	Male flowering time*	1.000	1.000
	Tassel branching*	0.925	0.908
	Number of grains per ear	0.832	0.622
	Grain length*	1.000	1.000
	Grain width*	0.995	0.984
	Grain weight*	1.000	0.999
	Grain color	0.717	0.689
	Stomata density*	0.999	0.999

a

: Traits displaying signal of selection (spatially-varying traits, S > 0.95) are indicated in bold,

  , 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 -4 )

2 10 -16

), a pattern not seen with geographically-correlated SNPs. That outlier SNPs displayed a greater isolation-byenvironment 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

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
	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).

structuring with 11 populations (instead of K=5) as covariate (M5'), 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)

Table ) .

 ) False positives likely represent a greater proportion of associations at K=5 as 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 [100]. In maize, this gene has a demonstrated role in kernel phenotypic differences between maize genetic groups [101]. 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.

	illustrated by a slight excess of small P-values when compared with a correction with
	11 populations for most traits (S10 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

[START_REF] Whitt | Genetic diversity and selection in the maize starch pathway[END_REF][START_REF] Jaenicke-Despres | Early allelic selection in maize as revealed by ancient DNA[END_REF]

. We found that Su1 was associated with variation at six traits (male and female flowering time, tassel branching, height With the 11-population correction, very low power is thus left to detect the effect of Su1 on phenotypes.

  Table). Primers sequences are available from the maize database project [115] and genotyping protocol were previously published [116]. Genotyping was done at the GENTYANE platform (UMR INRA 1095, Clermont-Ferrand, France). Allele calling was performed on electropherograms with the GeneMapper® Software Applied Biosystems®. Allele binning was carried out using Autobin software [117], and further checked manually.

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,

	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 on
	112 teosinte collections grown in five localities [124].
	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:

  where b n is the slope of the regression of Y on the n th structure covariate C n . Structure covariate values (C n 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 n are thus declared with ijklm indices, although they are purely genetic covariates.

Table. Description of 37 teosinte populations and sets of populations used in the present study by data types. S2 Table. List of the 18 phenotypic traits measured and estimates of narrow- sense heritabilities (h 2 ). S3 Table. Significance of main effects for each trait as determined by models M1 and M3.

Acknowledgments

We thank Jessica Melique for her help in gathering stomata data. We are very grateful to Angelica Cibrian at Langebio (Cinvestav, Irapuato, Mexico) who let us use her lab to lyophilize our samples, and store them. Valeria Souza greatly helped us with logistic support in Mexico. Insights from Delphine Legrand, Pierre de Villemereuil and Elodie Marchadier were important for analyzing correlations between traits and to estimate heritabilities. We would like to thank warmly five anonymous reviewers for their insightful comments on our work.

Material and Methods

Description of teosinte populations and sampling.

We used 37 teosinte populations of mexicana [START_REF] Slatkin | Rare alleles as indicators of gene flow[END_REF] and parviglumis [START_REF] Luquet | Phenotypic divergence of the common toad (Bufo bufo) along an altitudinal gradient: evidence for local adaptation[END_REF] subspecies from two previous collections [START_REF] Aguirre-Liguori | Connecting genomic patterns of local adaptation and niche suitability in teosintes[END_REF][START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF]112] 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 [113] [START_REF] Fustier | Signatures of local adaptation in lowland and highland teosintes from whole-genome sequencing of pooled samples[END_REF]. 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) for each gradient from environmental variables using the FactoMineR package in R [114] 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).

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 ASReml v.3.0 [125] and MM4LMM v2.0.1 [https://rdrr.io/cran/MM4LMM/man/MM4LMMpackage.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 M'1:

where mother plant within population is considered as fixed. 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 [114] from the residues of model M2 computed using the lm package in R:

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 [128]:

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 Ν(0,K σ 2 u) distribution, where K is the kinship matrix computed as described above.

A variant of model M5 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):

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 [134]. 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

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, [136]) implemented in the ecodist R package [137] 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.

Ethics Statement

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

Supporting information captions