The 59 Mexican landraces defined by CONABIO (2011) (http://www.biodiversidad.gob.mx/usos/maices/razas2012.html) were examined to select five (Zapalote Chico, Conejo, Zapalote Grande, Tehua, and Comiteco) that show contrasting phenotypic characteristics (Figure 1 and Table 1). These landraces are represented by 37 accessions from the Germplasm Bank of the International Maize and Wheat Improvement Center (CIMMYT) that were genotyped in this study (see below), plus 13 accessions from a previous study that used the same extraction protocol and SNP data generation (Arteaga et al., 2015). Our sample size for each landrace was 8–14 accessions, collected between 1946 and 2010 and covering most of the potential distribution of each landrace (see Supplementary Table 1 and Figure 2). Potential distributions were taken from Perales and Golicher (2011).

From CIMMYT's accessions, one seed per accession was germinated in greenhouse conditions. DNA was extracted from leaf tips collected from 3-week-old seedlings with Thermo Fisher Scientific Genomic DNA Purification Kit (K0512) (Thermo-Fisher Scientific, Waltham, MA, USA). DNA quality was determined by agarose gel electrophoresis and with a Thermo-Fisher Scientific Nanodrop Spectrophotometer 2000 (Thermo-Fisher Scientific, Waltham, MA, USA). Accessions with a 260/280 ratio below 1.6 were re-extracted.

Single nucleotide polymorphisms were genotyped at INMEGEN (Instituto Nacional de Medicina Genómica) with the Illumina MaizeSNP50 BeadChip on an Infinium HD assay (Illumina, San Diego, Ca, USA). Automated allele calling was done in GenomeStudio 2010.1 (Genotyping module 1.7.4; Illumina); loci with GenTrain score <0.3 and with more than 30% missing data were excluded. Data was exported to PLINK format (Purcell et al., 2007). A maximum of 10% missing data per individual was allowed. Data from the 13 accessions obtained from Arteaga et al. (2015) was added to this dataset. Only SNPs present in both groups of accessions were retained.

In order to guide the search for landrace-distinguishing SNPs, three sets of SNPs were defined: (1) all the SNPs available (35,909), (2) SNPs of domestication and improvement (974 SNPs), and (3) SNPs with high F_ST (435 SNPs). Set one refers to all the SNPs retained after filtering and joining our dataset with Arteaga's. The SNPs of domestication and improvement were determined using the maize domestication and improvement genes reported by Hufford et al. (2012) and Meyer and Purugganan (2013). To do this we used the Panzea Genotype Search Tool (http://cbsuss05.tc.cornell.edu/hdf5/hdf5.asp) to search which of our 35,909 SNPs were located in loci reported by Hufford et al. (2012) and Meyer and Purugganan (2013). SNPs with high F_ST were defined after calculating the value of F_ST among landraces for every SNP using the SNPStats package (Clayton, 2015) in R (R Core Team, 2008) and taking the 1% of SNPs with the highest F_ST.

We used the three previously defined set of SNPs to examine if they could be useful for distinguishing the landraces in a principal components analysis (PCA) and discriminant analysis of principal components (DPCA). The PCA was performed for each set using the package SNPRelate (Zheng et al., 2012). The first two components were used for plotting, coloring accessions by landrace. A K-means clustering analysis was performed to the three sets of SNPs to identify which set better recovered the landrace pattern using the packages adegenet (Jombart, 2008) and ape (Paradis et al., 2004). The DAPC (Jombart et al., 2010) was carried out to classify samples according to their landrace. This analysis was performed with the set of SNPs with high F_ST using the R package adegenet. The contribution of each SNP to the DAPC classification model was calculated. The SNP distribution showed a long-tailed distribution, the top contributing SNPs (~2%) were identified as landrace-informative SNPs.

To discard the possible confounding effect of altitude in the SNPs identified as landrace-informative, we repeated the previous analyses but grouping the accessions by altitude instead of landrace. For this, first the accessions were grouped according to their altitude in two groups: high altitude (>750 masl) and low altitude (0-750 masl). The limit between the altitudinal groups was defined following the altitudinal distribution of the landraces included, separating the interquartile range of the landraces typically growing at low altitude (Zapalote Chico and Zapalote Grande) from the interquantile range of the landraces growing at high altitude (Tehua and Comiteco). Still, all landraces except Zapalote Chico had accessions from both altitude categories. All Zapalote Chico accessions belonged to the low altitude group (see Table 1).

F_ST index of each SNP between these two altitude groups was calculated using the SNPStats R package, and a subset of the 1% of SNPs with highest F_ST was extracted to continue the analysis. Using this group of SNPs, a DAPC was performed as earlier, but so as to have a model that classifies the samples according to their altitude group. The contribution of each SNP to the classification model was obtained, and the SNPs with the highest contribution were identified as altitude-informative SNPs. The groups of landrace informative SNPs and altitude associated SNPs were compared to identify which landrace-informative SNPs are independent from altitude.

SNPs identified as landrace- or altitude-informative were reviewed in order to obtain their genomic position and biological annotation. This was accomplished by searching for the corresponding loci through the locus lookup tool at MaizeGDB (www.maizegdb.org/locus_lookup; Andorf et al., 2010). Further information about genes associated with the obtained SNPs was recovered through literature review.

All SNP data and R scripts used for this study is available in the Dryad online repository: http://dx.doi.org/10.5061/dryad.j2n8q.

After data filtering, the final SNPs data set consisted of 35,909 SNPs. The domestication and improvement SNP set consists of 974 SNPs located within the loci associated with domestication and improvement reported by Hufford et al. (2012) and Meyer and Purugganan (2013). A total of 84 SNPs were located in the loci reported by Meyer and Purugganan. In the loci reported by Hufford et al. 501 SNPs were located in loci associated with improvement and 389 SNPs in loci associated with domestication.

Values of F_ST among landraces for every SNP ranged from 0 to 0.597 (mean = 0.056, SD = 0.055). Most SNPs present F_ST values close to 0, indicating that they are poorly differentiated among landraces. The high F_ST SNP set consists of 435 SNPs with F_ST values higher than 0.2283, representing 1% of SNPs with highest F_ST (Figure 3). Most domestication and improvement SNPs have relatively small values of F_ST, and only 6 SNPs are common between the two sets.

PCA using each of the three sets of SNPs reflects that the genomic information each set carries is structured differently. Figure 4 shows the different distribution patterns of the first two principal component among the three SNP sets. The percent of variance explained by the first three PCs of the high F_ST SNP set is 20.91, 6.30, and 5.67, respectively, which is higher than the variance explained by the first three PCs of the set of all SNPs (4.39, 3.12, 2.89) and the domestication and improvement SNPs (5.04, 3.80, 3.60).

The utility of the three sets of SNPs for grouping the accessions according to their landrace was assessed by a k-means clustering analysis, testing from one to five clusters and determining the appropriate number of clusters through the Bayesian information criterion (BIC) (Figure 5). For the all SNPS and domestication and improvement SNPs sets, the best clustering groups all the accessions into a single cluster. The best clustering for the high F_ST SNP set identifies three groups: one group containing only accessions of Zapalote Chico, a second group with accessions of landraces Zapalote Grande and Conejo, and a third group with accessions of landraces Tehua and Comiteco. Since the only SNP set capable of identifying different clusters in the accessions was the high F_ST SNP set, we continued the analysis only with these SNPs.

The classification model from the DAPC analysis used the first four principal components of the high F_ST SNP set and was able to classify the accessions using two discriminant functions (Figure 6A). The first discriminant function (horizontal axis) recovers an altitude distribution pattern for the accessions, sorting them successively from the landrace growing in the lowest altitudes (Zapalote Chico) to the landrace growing in the highest places (Comiteco). The second discriminant function separates races growing in similar altitudes. The two discriminant functions are able to separate the clusters of accessions according to their landrace.

This model has a high degree of accuracy in assigning the correct landrace to each sample, 94% of the samples are correctly classified (Figure 6B). This accuracy is achieved despite the high degree of admixture present among some samples (Figure 6C). Tehua and Comiteco are the two landraces with more admixture between them; and accordingly, the ones in which the landrace reassignment by the model is less discriminant. Samples from Conejo and Zapalote Grande show certain degree of admixture, but the model is able to reassign their accessions accurately. Zapalote Chico is the landrace with less admixture, forming a clearly separated cluster from the rest. This classification model gains more validity considering that the samples we analyzed are representative of each landrace's geographic distribution and their time of collection spans from 1946 to 2010.

The contribution of each SNP to each discriminant function was calculated. For both functions, the SNP contributions showed an exponential distribution, with most SNPs having very little contribution and a few SNPs with very high contribution (Supplementary Figure 1). The top contributing SNPs of each function were identified, giving a group of 11 informative SNPs for the first function and 9 informative SNPs for the second function (Supplementary Table 2). Together, these 20 SNPs are identified as the landrace-informative SNPs. Genomic position and biological annotation of these SNPs was obtained: 9/20 are associated to seven genes, four of which have not yet been characterized. The three characterized genes are a polygalacturonase enzyme, heat shock factor 30 (HSF30), and transcription factor IIIA (Supplementary Table 2).

A set of SNPs differentiated between the two altitude groups was identified using F_ST index values. F_ST values between altitude groups for all SNPs ranged from 0 to 0.5251 (mean = 0.0345, SD = 0.043). Data for the 1% of SNPs with the highest F_ST (above 0.1728) was extracted. The best k-means grouping using this group of SNPs resulted in two groups, which almost perfectly recovered the high and low altitude groups previously defined, clustering only one accession incorrectly. The DAPC was done using the first principal component and one discriminant function (Supplementary Figure 2). The classification model had an almost perfect performance, misclassifying only that same accession that was incorrectly clustered before (Figure 6C). The contribution of each SNP to the discriminant function followed an exponential distribution. The top 14 contributing SNPs were identified as altitude-informative SNPs (Supplementary Figure 3), and their genomic position and biological annotation was obtained (Supplementary Table 3). Of the 14 altitude associated SNPs, eight fall within loci of genes of which only three genes have been annotated. These genes correspond to the auxin responsive protein SAUR40, Transformer-2, and the heat shock factor 30 (HSF30).

Comparing the altitude-informative SNPs with the landrace-informative SNPs, there are five SNPs shared between the two sets (Supplementary Table 2), two of which are associated to an uncharacterized gene (GRMZM2G100103). Additionally, both analyses identified a SNP, albeit a different one each, part of a gene associated with the heat shock factor protein 30. So this locus was also considered to be shared between the analyses. In total, for the landrace analysis there are 14 uniquely informative SNPs and six SNPs shared with the altitude analysis. For the altitude analysis, there are eight unique SNPs that are not shared with the landrace analysis.

Both the PCA and the k-means clustering analysis show that accession clustering depends on the set of SNPs used. In particular, the set of SNPs with high F_ST groups the samples according to their landrace while neither of the other two sets (“domestication and improvement SNPs” and “ all SNPs”) does.

The k-means clustering analysis grouped all the samples in a single cluster when using all the SNPs, meaning that the whole-genomic distribution of variation is not structured by landrace. In other words, from a whole-genome perspective, these five landraces are not distinguishable among them. This is consistent with previous studies showing that maize genomic structure is not explained by landrace, and is better correlated with environmental factors (van Heerwaarden et al., 2011; Breña Ochoa, 2013; Arteaga et al., 2015). Clustering analysis with the domestication and improvement SNPs also failed in finding structure among the samples, grouping them all in one cluster. Domestication SNPs are defined as loci with differentiated between teosinte and maize, thus involved in maize domestication from its wild relative. In this sense, it is expected that domestication SNPs are not differentiated among maize landraces. On the other hand, the improvement SNPs we used were found as loci differentiated between maize landraces and modern improved lines (Hufford et al., 2012). Modern improved lines derive from distinct landrace founders and have been subject to subsequent selection for agronomic traits, which are not the same that constitute the differences among these landraces. Thus, genomic variation in modern lines not necessarily maintains the structure of variation found on landraces, as shown by the low values of F_ST found in most domestication and improvement SNPs (Figure 3) indicates that these loci are not differentiated among landraces and so do not explain morphological differentiation.

Using the set of high F_ST SNPs (435 SNPs, Figure 2) recovers the landrace assemblages as recognized by morphology. This suggests that this set of SNPs are under different selective pressures depending on the landrace, therefore varying in a different way from the rest of the genome. These selective pressures are probably a combination of environmental differences among the areas the landraces are grown, and those guided by the smallholder farmers that continuously select seed looking for certain characteristics of the crop: desired morphology, taste, compatibility with their agricultural practices, among others (Pressoir and Berthaud, 2004; Boege, 2008; Bellon et al., 2011).

We were able to classify samples by landrace using the high F_ST SNPs. Sample classification was achieved with high accuracy (94%) despite broad landrace admixture. Most samples show some degree of mixing among landraces (Figure 6C), this pattern is explained by the known high degree of gene flow among maize population. There is a high degree of gene flow in maize populations. It is caused both by open pollination of maize fields flowering at the same time, and by seed transportation and exchange by farmers (Dyer and Taylor, 2008; Orozco-Ramírez et al., 2016). Gene flow among landraces is also strengthened by the fact that 85.5% of the total average surface cultivated with maize in Mexico (covering an average of 8.4 million ha between 1996 and 2006) is cultivated by smallholders using landraces (SIAP, 2008, 2016) therefore creating a large and genetically diverse landscape for open pollination to occur. Seed manipulation performed by producers maintains a phenotypic differentiation of maize landraces, in spite of high levels of gene flow, both by pollen and seed exchange (Pressoir and Berthaud, 2004; Dyer and Taylor, 2008; Mercer et al., 2008; Orozco-Ramírez et al., 2016).

Accurate models to classify accessions by landraces and altitude categories were attained through DAPC, allowing for the further inference of informative markers for both models. Informative SNPs were identified by their high contribution to the classification models, being the most informative SNPs when trying to classify by race or by altitude. SNPs association to landrace and altitude is based on the idea that if the model found these marker to be the most useful when separating the landrace/altitude groups, it is because there are underlying differences among groups in these markers, probably due to selection. We identified 14 landrace-informative SNPs that are independent from altitude. Even though our altitude analysis does not account for all the environmental factors, previous studies establish altitude as the principal factor shaping maize genomic variation (van Heerwaarden et al., 2011; Breña Ochoa, 2013; Arteaga et al., 2015; Romero Navarro et al., 2017), thus making it a valid approximation of environmentally driven natural selection.

A potential caveat of our analyses is that the genomic data used was obtained from the Illumina MaizeSNP50 BeadChip. This SNP chip was developed from the genome of maize B73 line to identify genomic sites of agronomic interest (Ganal et al., 2011). Therefore, it presents a bias toward identifying variation present in the B73 genome, and it is highly probable that there are other genomic regions differentiated among landraces that are not recognized by the SNP chip. In other words, our analysis may be missing variation at other alleles, which make up a substantial proportion of natural variation (Huang and Han, 2014).

The proposed landrace and altitude informative loci identified here represent an opportunity to further confirm if SNPs identified by genomic analyses can explain landraces phenotypic characteristics. This can be assessed by testing their direct influence in the plant phenotype, and analyzing their variation among other maize landraces. However, through searching for the functional annotation of associated genes, our approximation showed interesting results hinting functional importance of the identified loci.

One of the five annotated genes associated with the informative markers, the heat shock factor 30 HSF30, was shared between the landrace and the altitude analyses. HSF30 is part of the heat shock transcription factors that drive the cellular heat stress response in plants (Westerheide et al., 2012). HSF30 has not been functionally characterized in maize, but it is homologous to heat shock factor A2 present in other plant species where its function has been studied. For example, in Arabidopsis thaliana, HSF2A is the heat shock factor with highest expression under heat stress conditions, and acts as regulator of several stress-responsive genes (Schramm et al., 2006); and in Solanum lycopersicum, it has been shown that HSF2A transcription is strictly dependent on heat stress, and acts as a strong transcriptional activator for heat shock proteins (Nover et al., 2001). Considering the homology between HSF30 and HSF2A and the relationship between altitude and temperature, our results suggest that landraces hold a component of local adaptation to the altitudinal range where they are grown.

The other two annotated genes identified among the landrace-informative SNPs correspond to a polygalacturonase enzyme and transcription factor IIIA. An homologous gene of transcription factor IIIA has been characterized in rice, where its function was associated with plant development and tolerance to abiotic stress related to salinity, low temperatures, and drought (Huang et al., 2012). We consider transcription factor IIIA of special interest, because of the non-linear phenotypic effects a transcription factor can exert on a developing plant in the context of a transcription regulatory network (Davila-Velderrain and Alvarez-Buylla, 2014). Further maize research focusing on the effects of variation in this gene grown under different environmental conditions would be very interesting.

An interesting result in the classification model is that the samples of Zapalote Chico were the most distant cluster and presented a high homogeneity within the landrace (Figure 6). This has already been highlighted by a study focusing in the role of MuDR transposons, which are a characteristic of the Mutator maize line and are associated with a high rate of somatic mutations (Gutiérrez-Nava et al., 1998). Specifically, Gutierrez-Nava and collaborators reported that Zapalote Chico has an elevated rate of MuDR transposons, and that, when inter-crossed with germplasm from another landrace, the genomic elements present in Zapalote Chico cause somatic mutations, an effect known as hybrid dysgenesis. These mutations do not happen when two samples of Zapalote Chico are inter-crossed, showing that somehow the mutagenic activity of the transposons is repressed. This phenomenon explains Zapalote Chico's genetic isolation relative to other landraces. Interestingly, Zapotecs, the indigenous people growing this landrace, interpret this phenomenon in the form of a myth, which holds that when Zapalote Chico is crossed with another kind of maize it “kills" it (Gutiérrez-Nava et al., 1998).

In this study, we focused on genomic differences among maize landraces; we did not consider other factors that cause phenotypic differences and that could have a strong influence on crop phenotypes. Notably, the manipulation of the agroecological system by the people that grow maize is an important source of variation that could explain part of the differences seen among landraces (Jardón Barbolla and Benítez, 2012). Nonetheless, landraces, characterized in terms of both local adaptation and satisfaction of cultural preferences, have an associated genomic variation underlying their differences (Bellon et al., 2009). Several tens of millions of rural populations in Mexico depend solely on their maize production. Mexican cultural richness is reflected in the wide differences of uses people give to maize, resulting in a great variety in terms of flavor, shape, color, texture, and other organoleptic characteristics (Boege, 2008). Although natural and artificial selection exert their influence in particular loci, the concerted action of whole genomic variation and external conditions during a plant's development determine its phenotype, and consequently, overall success for its further sowing by the people using it.

Mexican maize is grown in complex agroecological systems where biological (interspecies relations, genetic flow, pollinators activities, soil microbiome, etc.), social, and cultural (farming techniques, seed interchange, land tenure patterns, religious, and symbolic uses, etc.) factors interact. This demands the study of maize and the agroecosystems in which it grows, including the groups of small-scale maize producers, as a whole for the understanding of maize diversity.

Mexican maize landraces contain a high genomic diversity within and among them, which results from the underlying multiplicity of environmental and biocultural factors mediating their existence (Arteaga et al., 2015). These factors entail selective pressures, both natural/environmental and artificial/biocultural, acting on different genomic regions. In this study we found that, even though from a whole-genome perspective maize landraces are indistinguishable, it is possible to identify genomic regions differentiated among landraces. Furthermore, combining the landrace analysis with an environmental association analysis on the genomic data, it is possible to recognize loci influenced by natural or artificial selection. Our results are evidence of the genomic variation that exists in maize landraces and their importance for maize diversity conservation. Nevertheless, it is important to keep in mind that landraces are evolving entities, sustained through a process based on a great genomic, environmental and bio-cultural diversity. Furthermore, focusing on preserving landraces by monitoring their phenotypic diversity through a small number of representative samples may lead to losing the genetic diversity underlying local adaptation. Understanding the genomic foundation of landrace diversity and the processes influencing it is imperative for its maintenance. Given the different factors shaping maize biodiversity, a program aimed at conserving the genomic and phenotypic richness of maize requires the maintenance of the biocultural processes in which it is grown and used.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer SR and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.