Natural variation in teosinte at the domestication locus teosinte branched1 (tb1)

Sampling and genotyping

We sampled all individuals and populations that were available to us, consisting of 1,110 individuals from 350 populations (247 maize landraces, 17 mexicana populations, and 86 parviglumis populations) and assessed the presence or absence of the Hopscotch insertion (Table S1 and Table S2). Numbers of individuals sampled per population ranged from 1–43 for parviglumis, 1–35 for mexicana, and 1–18 for the maize landrace populations. Available samples did not allow us to sample evenly from populations, but did allow us to calculate Hopscotch frequency in a subset of populations, as well as elucidate the geographic distribution of the Hopscotch across multiple independent sampling sites. DNA was extracted from leaf tissue using a modified CTAB approach (Doyle & Doyle, 1990; Maloof et al., 1984). We designed primers using PRIMER3 (Rozen & Skaletsky, 2000) implemented in Geneious (Kearse et al., 2012) to amplify the entire Hopscotch element, as well as an internal primer allowing us to simultaneously check for possible PCR bias between presence and absence of the Hopscotch insertion due to its large size (∼5 kb). Two PCRs were performed for each individual, one with primers flanking the Hopscotch (HopF/HopR) and one with a flanking primer (HopF) and an internal primer (HopIntR). Primer sequences are HopF, 5′-TCGTTGATGCTTTGATGGATGG-3′; HopR, 5′-AACAGTATGATTTCATGGGACCG-3′; and HopIntR, 5′-CCTCCACCTCTCATGAGATCC-3′ (Fig. 1 and Fig. S1). Homozygotes for the no-Hopscotch allele show a single band for absence of the element (∼300 bp) produced by the HopF/HopR primer set, and 0 bands for the HopF/HopIntR primer set since they lack the LTR where the internal primer sequence is located. Homozygotes for the Hopscotch allele also show one band at 5 kb for the HopF/HopR PCR product as well as one band at 1.1 kb for the HopF/HopIntR PCR. Heterozygotes for the Hopscotch allele show three bands total; both a 300 bp band and a 5 kb band for the HopF/HopR PCR and a 1.1 Kb band for the HopF/HopIntR PCR (Table S2). Since we developed a PCR protocol for each allele, if only one PCR resolved well, we scored one allele for that individual rather than infer the diploid genotype. We used Phusion High Fidelity Enzyme (Thermo Fisher Scientific Inc., Waltham, Massachusetts, USA) and the following conditions for amplifications: 98 °C for 3 min, 30 cycles of 98 °C for 15 s, 65 °C for 30 s, and 72 °C for 3 min 30 s, with a final extension of 72 °C for 10 min. PCR products were visualized on a 1% agarose gel and scored for presence/absence of the Hopscotch based on band size.

Genotyping analysis

To calculate differentiation between populations (F_ST) and subspecies (F_CT) we used HierFstat (Goudet, 2005). These analyses only included populations (n = 32) in which eight or more chromosomes were sampled. To test the hypothesis that the Hopscotch insertion may be adaptive under certain environmental conditions, we looked for significant associations between Hopscotch frequency and environmental variables using the software BayEnv (Coop et al., 2010). BayEnv creates a covariance matrix of relatedness between populations and then tests a null model that allele frequencies in populations are determined by the covariance matrix of relatedness alone against the alternative model that allele frequencies are determined by a combination of the covariance matrix and an environmental variable, producing a posterior probability (i.e., Bayes Factor; Coop et al., 2010). We used teosinte (ssp. parviglumis and ssp. mexicana) genotyping and covariance data from Pyhäjärvi et al. (2013) for BayEnv, with the Hopscotch insertion coded as an additional biallelic marker. SNP data from Pyhäjärvi et al. (2013) were obtained using the MaizeSNP50 BeadChip and Infinium HD Assay (Illumina, San Diego, CA, USA) and phased using the program fastPHASE (Scheet & Stephens, 2006). Environmental data were previously obtained from www.worldclim.org and soil data were downloaded from the Harmonized World Soil Database (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012) at www.harvestchoice.org. Environmental data represent average values for the last several decades (climatic data) or are likely stable over time (soil data) and therefore represent conditions important for local adaptation of our samples. Information from these data sets was summarized by principle component analysis following Pyhäjärvi et al. (2013).

Sequencing

In addition to genotyping, we chose a subset of parviglumis individuals for sequencing. We chose twelve individuals from each of four populations from Jalisco state, Mexico (San Lorenzo, La Mesa, Ejutla A, and Ejutla B). For amplification and sequencing, we selected two regions approximately 600 bp in size from within the 5′ UTR of tb1 (Region 1) and from 1,235 bp upstream of the start of the Hopscotch (66,169 bp upstream from the start of the tb1 ORF; Region 2). We designed the following primers using PRIMER3 (Rozen & Skaletsky, 2000): for the 5′ UTR, 5′-GGATAATGTGCACCAGGTGT-3′ and 5′-GCGTGCTAGAGACACYTGTTGCT-3′; for the 66 kb upstream region, 5′-TGTCCTCGCCGCAACTC-3′ and 5′-TGTACGCCCGCCCCTCATCA-3′ (Table S1, See Supplemental Materials with the online version of this article). We used Taq polymerase (New England Biolabs Inc., Ipswich, Massachusetts, USA) and the following thermal cycler conditions to amplify fragments: 94 °C for 3 min, 30 cycles of 92 °C for 40 s, annealing for 1 min, 72 °C for 40 s, and a final 10 min extension at 72 °C. Annealing temperatures for Region 1 and Region 2 were 59.7 °C and 58.8 °C, respectively. To clean excess primer and dNTPs we added two units of Exonuclease1 and 2.5 units of Antarctic Phosphatase to 8.0 µL of amplification product. This mix was placed on a thermal cycler with the following program: 37 °C for 30 min, 80 °C for 15 min, and a final cool-down step to 4 °C.

We cloned cleaned fragments into a TOPO-TA vector (Life Technologies, Grand Island, New York, USA) using OneShot TOP10 chemically competent E. coli cells, with an extended ligation time of 30 min for a complex target fragment. We plated cells on LB agar plates containing kanamycin, and screened colonies using vector primers M13 Forward and M13 Reverse under the following conditions: 96 °C for 5 min; then 35 cycles at 96 °C for 30 s, 53 °C for 30 s, 72 °C for 2 min; and a final extension at 72 °C for 4 min. We visualized amplification products for incorporation of our insert on a 1% agarose TAE gel.

Amplification products with successful incorporation of our insert were cleaned using Exonuclease 1 and Antarctic Phosphatase following the procedures detailed above, and sequenced with vector primers M13 Forward and M13 Reverse using Sanger sequencing at the College of Agriculture and Environmental Sciences (CAES) sequencing center at UC Davis. We aligned and trimmed primer sequences from resulting sequences using the software Geneious (Kearse et al., 2012). Following alignment, we verified singleton SNPs by sequencing an additional one to four colonies from each clone. If the singleton was not present in these additional sequences it was considered an amplification or cloning error, and we replaced the base with the base of the additional sequences. If the singleton appeared in at least one of the additional sequences we considered it a real variant and kept it for further analyses.

Sequence analysis

For population genetic analyses of sequenced Region 1 and sequenced Region 2 we used the Analysis package from the Libsequence library (Thornton, 2003) to calculate pairwise F_ST between populations and to calculate standard diversity statistics (number of haplotypes, haplotype diversity, Watterson’s estimator ${\hat{\theta }}_W$, pairwise nucleotide diversity ${\hat{\theta }}_{\pi }$, and Tajima’s D). Significance of Tajima’s D results was gauged by comparing empirical data to 10,000 coalescent simulations conducted using the program ms (Hudson, 2002) under a standard neutral model based on observed estimates of the population mutation rate theta and assuming an identical value for the population recombination rate rho. Empirical results falling outside the 95% confidence interval of our simulated data were deemed significant. To produce a visual representation of differentiation between sequences and examine patterns in sequence clustering by Hopscotch genotype, we used Phylip (http://evolution.genetics.washington.edu/phylip.html) to create neighbor-joining trees with bootstrap-supported nodes (10,000 repetitions). For creation of trees we also included homologous sequence data from Maize HapMapV2 (Chia et al., 2012) for teosinte inbred lines (TILs), some of which are known to be homozygous for the Hopscotch insertion (TIL03, TIL17, TIL09), as well as 59 lines of domesticated maize.

Introgression analysis

In order to assess patterns of linkage disequilibrium (LD) around the Hopscotch element in the context of chromosomal patterns of LD we used Tassel (Bradbury et al., 2007) and calculated LD between SNPs across chromosome 1 using previously published data from twelve plants each of the Ejutla A (EjuA), Ejutla B (EjuB), San Lorenzo (SLO), and La Mesa (MSA) populations (Pyhäjärvi et al., 2013). We chose these populations because we had both genotyping data for the Hopscotch as well as chromosome-wide SNP data for chromosome 1. For each population, we filtered the initial set of 5,897 SNPs on chromosome 1 to accept only SNPs with a minor allele frequency of at least 0.1, resulting in 1,671, 3,023, 3,122, and 2,167 SNPs for SLO, EjuB, EjuA, and MSA, respectively. We then used Tassel (Bradbury et al., 2007) to calculate linkage disequilibrium (r²) across chromosome 1 for each population.

We examined evidence of introgression on chromosome 1 in these same four populations (EjuA, EjuB, MSA, SLO) using STRUCTURE (Falush, Stephens & Pritchard, 2003) and phased data from Pyhäjärvi et al. (2013), combined with the corresponding SNP data from a diverse panel of 282 maize lines (Cook et al., 2012). SNPs were anchored in a modified version of the IBM genetic map (Gerke et al., 2013). Since STRUCTURE does not account for LD due to physical linkage we created haplotype blocks using a custom Perl script from Hufford et al. (2013), code available at http://dx.doi.org/10.6084/m9.figshare.1165577. In maize, LD decays over an average distance of 5,500 bp (Chia et al., 2012); because LD decay is even more rapid in teosinte (Chia et al., 2012) we used a conservative haplotype block size of 5 kb. We ran STRUCTURE at K = 2 under the linkage model, with the assumption being that individuals fall into either a maize or teosinte cluster, performing three replicates with an MCMC burn-in of 10,000 steps and 50,000 steps post burn-in.

Phenotyping of parviglumis

To investigate the phenotypic effects of the Hopscotch insertion in teosinte we conducted a phenotyping trial in which we germinated 250 seeds of parviglumis collected in Jalisco state, Mexico (population San Lorenzo; Hufford, 2010) where the Hopscotch insertion is segregating at highest frequency (0.44) in our initial genotyping sample set. In order to increase our chances of finding the Hopscotch in our association population we selected seeds from sites within the population where genotyped individuals were homozygous or heterozygous for the insertion. We chose between 10–13 seeds from each of 23 sampling sites. We treated seeds with Captan fungicide (Southern Agricultural Insecticides Inc., Palmetto, Florida, USA) and germinated them in petri dishes with filter paper. Following germination, 206 successful germinations were planted into one-gallon pots with potting soil and randomly spaced one foot apart on greenhouse benches. Plants were watered three times a day with an automatic drip containing 10-20-10 fertilizer, which was supplemented with hand watering on extremely hot and dry days.

Starting on day 15, we measured tillering index as the ratio of the sum of tiller lengths to the height of the plant (Briggs et al., 2007). Following initial measurements, we phenotyped plants for tillering index every 5 days through day 40, and then on day 50 and day 60. On day 65 we measured culm diameter between the third and fourth nodes of each plant. Following phenotyping we extracted DNA from all plants using a modified SDS extraction protocol. We genotyped individuals for the Hopscotch insertion following the PCR protocols listed above.

Tillering index data for each genotypic class did not meet the criteria for a repeated measures ANOVA, so we transformed the data with a Box–Cox transformation (λ = 0) in the Car Package for R (Fox & Weisberg, 2011) to improve the normality and homogeneity of variance among genotype classes. We analyzed relationships between genotype and tillering index and tiller number using a repeated measures ANOVA through a general linear model function implemented in SAS v.9.3 (SAS Institute Inc., Cary, North Carolina, USA). Additionally, in order to compare any association between Hopscotch genotype and tillering and associations at other presumably unrelated traits, we performed an ANOVA between culm diameter and genotype using the same general linear model in SAS. Culm diameter is not believed to be correlated with tillering index or variation at tb1 and is used as our independent trait for phenotyping analyses. SAS code used for analysis is available at http://dx.doi.org/10.6084/m9.figshare.1166630.

Genotyping for the Hopscotch insertion

The genotype at the Hopscotch insertion was confirmed with two PCRs for 837 individuals of the 1,100 screened (Table S1 and Table S2). Among the 247 maize landrace accessions genotyped, all but eight were homozygous for the presence of the insertion. Within our parviglumis and mexicana samples we found the Hopscotch insertion segregating in 37 (n = 86) and four (n = 17) populations, respectively, and at highest frequency within populations in the states of Jalisco, Colima, and Michoacán in central-western Mexico (Fig. 2). Using our Hopscotch genotyping, we calculated differentiation between populations (F_ST) and subspecies (F_CT) for populations in which we sampled sixteen or more chromosomes. We found that F_CT = 0, and levels of F_ST among populations within each subspecies (0.22) and among all populations (0.23) (Table 1) are similar to genome-wide estimates from previous studies Pyhäjärvi et al., 2013. Although we found large variation in Hopscotch allele frequency among our populations, BayEnv analysis did not indicate a correlation between the Hopscotch insertion and environmental variables (all Bayes Factors <1).

Sequencing upstream regions of the tb1 ORF

To investigate patterns of sequence diversity and linkage disequilibrium (LD) in the tb1 region and any evidence of selection on this locus, we sequenced two small (<1 kb) regions upstream of the tb1 ORF in four populations from the Jalisco region. After alignment and singleton checking we recovered 48 and 40 segregating sites for the 5′ UTR region (Region 1) and the 66 kb upstream region (Region 2), respectively. For Region 1, Ejutla A has the highest values of haplotype diversity and ${\hat{\theta }}_{\pi }$, while Ejutla B and La Mesa have comparable values of these summary statistics, and San Lorenzo has much lower values. Additionally, Tajima’s D is significantly negative in the two Ejutla populations and La Mesa, but is closer to zero in San Lorenzo (Table 2). For Region 2, haplotype diversity and ${\hat{\theta }}_{\pi }$, are similar for Ejutla A and Ejutla B, while La Mesa and San Lorenzo have slightly lower values for these statistics (Table 2). Tajima’s D is positive in all populations except La Mesa, where a slightly negative value suggests a slight excess of low frequency variants (Table 2). Pairwise values of F_ST within population pairs Ejutla A/Ejutla B and San Lorenzo/La Mesa are close to zero for both sequenced regions as well as for the Hopscotch, while they are high for other population pairs (Table 1).

Evidence of introgression around the tb1 region

We investigated the possibility of introgression as an explanation for the frequency of the Hopscotch allele in populations of teosinte using previously collected SNP data from Pyhäjärvi et al. (2013). The highest frequency of the Hopscotch insertion in teosinte was found in parviglumis sympatric with cultivated maize. Our initial hypothesis was that the high frequency of the Hopscotch element in these populations could be attributed to introgression from maize into teosinte. To investigate this possibility, we examined overall patterns of linkage disequilibrium across chromosome 1 and specifically in the tb1 region. If the Hopscotch is found in these populations due to recent introgression from maize, we would expect to find large blocks of linked markers near this element. We find no evidence of elevated linkage disequilibrium between the Hopscotch and SNPs surrounding the tb1 region in our resequenced populations (Fig. 3), and r² in the tb1 region does not differ significantly between populations with (average r² of 0.085) and without (average r² = 0.082) the Hopscotch insertion. In fact, average r² is lower in the tb1 region (r² = 0.056) than across the rest of chromosome 1 (r² = 0.083; Table 3).

Neighbor joining trees of our sequence data and data from the teosinte inbred lines (TILs; data from Maize HapMapV2, Chia et al., 2012) do not reveal any clear clustering pattern with respect to population or Hopscotch genotype (Fig. S2, see Supplemental Information with the online version of this article); individuals within our sample that have the Hopscotch insertion do not group with the teosinte inbred lines or domesticated maize that have the Hopscotch insertion. The lack of clustering of Hopscotch genotypes in our NJ tree as well as the lack of LD around tb1 do not support the hypothesis that the Hopscotch insertion in these populations of parviglumis is the result of recent introgression. However, to further explore this hypothesis we performed a STRUCTURE analysis using Illumina MaizeSNP50 data from four of our parviglumis populations (EjuA, EjuB, MSA, and SLO) (Pyhäjärvi et al., 2013) and the maize 282 diversity panel (Cook et al., 2012). The linkage model implemented in STRUCTURE can be used to identify ancestry of blocks of linked variants which would arise as the result of recent admixture between populations. If the Hopscotch insertion is present in populations of parviglumis as a result of recent admixture with domesticated maize, we would expect the insertion and linked variants in surrounding sites to be assigned to the “maize” cluster in our STRUCTURE runs, not the “teosinte” cluster. In all runs, assignment to maize in the tb1 region across all four parviglumis populations is low (average 0.017) and much below the chromosome-wide average (0.20; Table 4 and Fig. 4).

Phenotyping of Zea mays ssp. parviglumis

To assess the contribution of tb1 to phenotypic variation in tillering in a natural population, we grew plants from seed sampled from the San Lorenzo population of parviglumis, which had a high mean frequency (0.44) of the Hopscotch insertion based on our initial genotyping. We measured tiller number and tillering index, the ratio of the sum of tiller lengths to plant height, for 206 plants from within the San Lorenzo population, and genotyped plants for the Hopscotch insertion. We also measured culm diameter, a phenotype that differs between maize and teosinte but has not been shown to be affected by the Hopscotch insertion (Briggs et al., 2007). Culm diameter is meant to be an independent trait against which we can compare patterns of tillering index x Hopscotch genotype data. If tillering index in parviglumis is affected by the Hopscotch insertion, the expectation is that patterns of tillering index data will have a significant correlation with Hopscotch genotype, whereas we should find no significant correlation between culm diameter and Hopscotch genotype. Phenotypic data are available at http://dx.doi.org/10.6084/m9.figshare.776926. Our plantings produced 82 homozygotes for the Hopscotch insertion at tb1, 104 heterozygotes, and 20 homozygotes lacking the insertion; these numbers do not deviate from expectations of Hardy-Weinberg equilibrium. After performing a repeated measures ANOVA between our transformed tillering index data and Hopscotch genotype, we find no significant correlation between genotype at the Hopscotch insertion and tillering index (Fig. 5), tiller number, or culm diameter. Only on day 40 did we observe a weak but statistically insignificant (r² = 0.02, p = 0.0848) correlation between tillering index and the Hopscotch genotype, although in the opposite direction of that expected, with homozygotes for the insertion showing a higher tillering index.