Allelic Variation in Developmental Genes and Effects on Winter Wheat Heading Date in the U.S. Great Plains

Heading date in wheat (Triticum aestivum L.) and other small grain cereals is affected by the vernalization and photoperiod pathways. The reduced-height loci also have an effect on growth and development. Heading date, which occurs just prior to anthesis, was evaluated in a population of 299 hard winter wheat entries representative of the U.S. Great Plains region, grown in nine environments during 2011–2012 and 2012–2013. The germplasm was evaluated for candidate genes at vernalization (Vrn-A1, Vrn-B1, and Vrn-D1), photoperiod (Ppd-A1, Ppd-B1 and Ppd-D1), and reduced-height (Rht-B1 and Rht-D1) loci using polymerase chain reaction (PCR) and Kompetitive Allele Specific PCR (KASP) assays. Our objectives were to determine allelic variants known to affect flowering time, assess the effect of allelic variants on heading date, and investigate changes in the geographic and temporal distribution of alleles and haplotypes. Our analyses enhanced understanding of the roles developmental genes have on the timing of heading date in wheat under varying environmental conditions, which could be used by breeding programs to improve breeding strategies under current and future climate scenarios. The significant main effects and two-way interactions between the candidate genes explained an average of 44% of variability in heading date at each environment. Among the loci we evaluated, most of the variation in heading date was explained by Ppd-D1, Ppd-B1, and their interaction. The prevalence of the photoperiod sensitive alleles Ppd-A1b, Ppd-B1b, and Ppd-D1b has gradually decreased in U.S. Great Plains germplasm over the past century. There is also geographic variation for photoperiod sensitive and reduced-height alleles, with germplasm from breeding programs in the northern Great Plains having greater incidences of the photoperiod sensitive alleles and lower incidence of the semi-dwarf alleles than germplasm from breeding programs in the central or southern plains.

Introduction by epistatic interactions among genes or QTL [19] and may be affected by gibberellin insensitivity through the reduced-height loci [20].
The reduced-height loci are known to affect growth and development in winter wheat [21]. Reduced-height genes encode DELLA proteins, which are growth repressors that are degraded by a process that involves gibberellin [22]. The reduced-height alleles Rht-B1b and Rht-D1b each have a single SNP that causes a premature stop codon, resulting in reduced sensitivity to gibberellic acid and therefore, semi-dwarf stature. Gibberellin insensitivity has been associated with earlier heading date in wheat, and Rht-B1b and Rht-D1b are the predominant semi-dwarf alleles in wheat worldwide [20]. Rht-B1b has been reported as the leading cause of semi-dwarf stature among winter wheat in the U.S. Great Plains region [23]. Guedira et al. [23] found the reduced-height alleles Rht-B1b and Rht-D1b to be widespread in U.S. germplasm, and Rht8c to be less common but linked with the photoperiod insensitive allele Ppd-D1a.
The allelic diversity of candidate loci known to affect heading date has been characterized in several worldwide [20,24] and regional [25][26][27] collections of spring and winter wheat. Analyses of core collections of germplasm have indicated substantial variation in the geographic distribution of vernalization and photoperiod alleles. Variation in the presence and distribution of alleles and haplotypes has been found to vary among continents e.g., [24] and countries [26], with some haplotypes under-represented or absent from particular geographic regions. A more complete understanding of the genetic controls-including the allelic variants and effects of single genes, and distribution of favorable multi-locus genotypes-is important for plant breeders to prepare for future climate scenarios including those that are more variable or extreme than today's conditions [28][29].
We examined allelic variation, distribution, and effects of vernalization, photoperiod, and reduced-height loci in a collection of hard winter wheat germplasm representative of the U.S. Great Plains. This is the first article to determine allelic diversity at the photoperiod loci in winter wheat germplasm from the U.S. Great Plains. The objectives of the experiment were 1. To characterize allelic variants at loci known to affect flowering time, including vernalization, photoperiod, and reduced-height genes; 2. To assess the effect of allelic variants on heading dates in nine U.S. Great Plains environments; and 3. To investigate temporal and geographic distribution of alleles and haplotypes among germplasm developed during different periods or from different regions of the U.S. where hard winter wheat is grown.

Germplasm
The germplasm was a collection of 299 winter wheat genotypes (or "entries") comprising the Triticeae Coordinated Agricultural Project (http://www.triticeaecap.org) hard winter wheat association mapping panel (S1 Table). The entries included modern and historic cultivars and experimental breeding lines. Public breeding programs contributed 270 entries, private breeding programs contributed 27, and two were historic cultivars developed before 1900. The publically developed germplasm came from Nebraska (55 entries

Environments and Experimental Design
The entries were evaluated in nine field trials conducted across the U.S. Great Plains region during 2011-2012 and 2012-2013, as described previously [30] and in Table 1. Each environment was a combination of location, year, and moisture treatment, with three possible moisture regimes: full irrigation, partial irrigation, and rainfed (no supplemental irrigation). Environments grown in 2011-2012 include rainfed in Bushland, TX (Bu12R), full irrigation in Greeley, CO (Gr12F), partial irrigation in Greeley, CO (Gr12P), rainfed in Ithaca, NE (It12R), and Manhattan, KS (Ma12). The 2012-2013 environments include rainfed in Ardmore, OK (Ar13R), Fort Collins, CO (Fo13), rainfed in Hays, KS (Ha13R), and rainfed in Ithaca, NE (It13R). Different environments at the same location during the same year were side-by-side treatments with a moisture differential. There were separate rainfed and irrigated treatments at Fo13 and Ma12, but irrigation was not applied until after heading date so there were no significant treatment effects on heading date and the moisture treatments were treated as replications.
Four environments (Bu12R, Gr12P, Gr12F, Ha13R) were unreplicated and arranged in an augmented row-column design with two check varieties. The experimental entries were unreplicated except for 'Wichita' (CI 11952), which was included in the panel twice. The two check varieties were replicated 15 times each and systematically placed throughout the field. Fo13 had a similar experimental design but included two replications. The check varieties at Bu12R, Fo13, Gr12P, Gr12F, and Ha13R were 'Hatcher' [31] and 'Settler CL' [32], and these varieties were also included as experimental entries in the trial. Irrigation was applied at Gr12P and Gr12F using drip irrigation. Gr12P was irrigated less frequently and at a reduced volume than Gr12F. Irrigation totaled 101.6 mm at Gr12P and 335.3 mm at Gr12F. It12R and It13R used a similar experimental design but plots were arranged as 15 incomplete blocks, with one plot within each block planted to each of the check varieties Settler CL and 'Jagger' [33]. Both It12R and It13R included four replications, and these trials are described in detail by Guttieri et al. [34]. Ma12 had four replications of a modified row-column design arranged as incomplete blocks using a single, locally adapted check variety, 'Everest' (http://kswheatalliance.org/ varieties/everest/). The Ar13R trial was arranged as a randomized complete block design with two replications. Harvested plot area ranged from 2.2 to 4.6 m 2 .

Phenotypic Evaluation
Crop development was determined for each field plot using Zadoks' scale [35]. Heading date (stage 59) was recorded when the spike had fully emerged from the flag leaf sheath in approximately 50% of tillers. Days to heading were recorded as the number of days from January 1 to heading date. Growing degree-days from 1 January to heading were determined previously, and found to be strongly correlated with days from January 1 to heading [30]. However, number of days to heading was used in these analyses because estimates of allelic effects in days are more widely understood than estimates in°C days.

Genetic Evaluation
Genotypes of nine loci associated with flowering time were obtained primarily from KASP analysis, as described below. Those results were supplemented with polymerase chain reaction (PCR)-based analysis from a previous published report [21] and unpublished data from the Byrne lab at Colorado State University and the USDA-ARS Small Grains Genotyping lab in Manhattan, KS.
DNA preparation for genetic analyses. Genomic DNA was extracted at Colorado State University from leaf tissue from single seedlings using the phenol-chloroform method modified slightly from Riede and Anderson [36].
The spring allele Vrn-A1b was distinguished from Vrn-A1a and the vrn-A1 winter allele using codominant marker wMAS000034 [4]. Two additional markers were evaluated to determine copy number variation (CNV) of the winter allele vrn-A1: Vrn-A1_Exon4_C/T and Vrn-A1_Exon7_G/A [9]. The C allele at Vrn-A1_exon4 is associated with two or fewer copies, and the G allele at Vrn-A1_exon7 is associated with a single copy, so entries were classified as having three or more copies, two copies, or one copy of vrn-A1. A PCR assay for Vrn-A1 [45] was also run that detected the same SNP polymorphism and distinguished two winter allelesrenamed vrn-A1w and vrn-A1v by Eagles et al. [46]-associated with winter dormancy release and freezing tolerance. Thus, the winter alleles detected at Vrn-A1 using PCR correspond with CNV detected using KASP.
Three KASP markers used to distinguish spring alleles at Vrn-B1. TaVrn-B1_D-I, wMAS000037, and Vrn-B1_C are codominant markers that detect the Vrn-B1a, Vrn-B1b [47], and Vrn-B1c [48] spring alleles, respectively. An additional codominant marker, TaVrn-B1_1752, detected an A/G polymorphism in intron 1 of the vrn-B1 gene associated with Averaged across treatments' indicates that separate side-by-side rainfed and full-irrigation treatments at this location did not differ significantly for average heading date, so were treated as two replications. 2 Rainfed treatment was harvested on 18 July 2013 and fully irrigated treatment was harvested on 22 July 2013. differences in vernalization requirement duration [49]. The spring and winter alleles at Vrn-D1 were distinguished using a single dominant marker, wMAS000039 [43].
Photoperiod insensitive allele Ppd-A1a.1 was assayed with the marker TaPpd-A1prodel, which detects a deletion characteristic of the insensitive allele [50]. Alleles at Ppd-B1 were distinguished using two markers: wMAS000027, which detects the 'Chinese Spring'-type insensitive allele with a truncated copy and TaPpdBJ003, which identifies the 'Sonora 64'-type insensitive allele based on the presence of an intercopy junction [16]. Photoperiod sensitive and insensitive alleles at Ppd-D1 were distinguished using a single codominant marker, wMAS000024 that detects a deletion upstream of the coding region responsible for the photoperiod insensitive phenotype [15].
Single markers were also used to detect point mutations at the reduced-height loci Rht-B1 and Rht-D1 [51]. Mutants at Rht-B1 were genotyped using wMAS000001 that detected the causative SNP for semi-dwarf stature. Likewise, a single marker, wMAS000002 was used to detect the polymorphism at Rht-D1 that is associated with semi-dwarf stature. Additional information about KASP assays having wMAS designations is available at http://www. cerealsdb.uk.net/cerealgenomics/CerealsDB/kasp_download.php. Table 2. Description of photoperiod (Ppd), reduced-height (Rht), and vernalization (Vrn) loci, alleles, and phenotypes evaluated in this study.

Statistical Analyses
Best linear unbiased predictors were calculated separately for each environment based on field design and spatial trends using SAS 9.3 (SAS Institute, Inc., Cary, NC). Six different spatial correlation models (row-column, spherical, exponential, power, anisotropic power, and Matérn, [52] were tested for environments with an augmented row-column design (Bu12R, Fo13, Gr12P, Gr12F, Ha13). The best model was selected based on the AIC fit statistic. Replications were treated as a random effect for environments with multiple replications (Ar13R, Fo13, It12R, It13R, Ma12). Further statistical analyses were performed using the R version 3.1.3 [53]. Combined analyses of the effects of alleles at single genes were evaluated on all entries with homozygous allele calls using the 'car' package [54]. The ANOVA model terms consisted of the environment, gene, and when significant, interaction between the gene and environment. All terms were fit as fixed effects. Entries with heterozygous calls were treated as missing. When a gene had more than two alleles, pairwise comparisons were run using the 'lsmeans' package [55] to test Table 3. Description of multi-locus genotypes at photoperiod and reduced-height loci. The photoperiod (Ppd) alleles are 'a' insensitive and 'b' sensitive. The reduced-height (Rht) alleles are 'a' tall wild type and 'b' semi-dwarf. 'Het' is heterozygous at the locus. Of 299 total entries, 285 entries have complete genotypic data at all five loci, and 280 have complete data with all homozygous allele calls. differences between each pair of alleles. Individual environments were analyzed using the linear model (lm) function in the 'stats' package [53] to evaluate the proportion of variability in heading date explained by one or more gene, and to test for interactions between genes. Little's missing completely at random test (MCAR) was conducted using the 'BaylorEdPsych' package [56] to ensure that genotypic data were missing at random. The MCAR test was evaluated on all entries with at least one genotyped allele (297 entries total) and all loci with varying allele calls (Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, Rht-D1, vrn-A1 (winter alleles and CNV), and vrn-B1 (winter alleles)); the data were found to be missing at random (Χ 2 = 5.40, P = 1.00). The MCAR test was also evaluated on a subset of 293 entries with one or more genotype at Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, or Rht-D1, which were also found to be missing data at random (Χ 2 = 3.98, P = 0.91). Missing genotypic data are believed to be from poor quality or low quantity of DNA for marker assays, and not due to variation in heading date.
Models that tested effects of multiple genes included all entries with complete genotypic data at those loci. For model comparison a subset of 280 entries with complete genotypic data and all homozygous calls at the photoperiod and reduced-height loci were analyzed (Table 3). First, a full model that included the main effects of all five genes and all ten two-way interactions were fit to each data set. In the combined analysis, the main effect of the environment, two-way interactions between the environment and each gene, and ten three-way interactions between the environment and pairs of genes were also included in the model. Then, model selection was applied in R using the 'MuMIn' package [57] using the 'dredge' function, and the best-fit model was identified based on lowest AIC value. In some cases the best-fit model included non-significant terms, but removing these terms detracted from model fit. The relative importance of main effects in the best-fit model was determined through model averaging, using Akaike weights weighted across all evaluated models. The percent variation accounted for by each model term was calculated as the percent of sums of squares for a given term relative to the total sums of squares or to the sums of squares for all genetic terms in the model.

Allelic Diversity of Candidate Genes
The 299 winter wheat entries were genotyped at 11 candidate genes, but sample size varied among the loci due to different amounts of missing data at each gene. Genotypic data were missing from up to 14 entries per locus. The genotypes at each candidate gene for every entry are provided in S1 Table. Winter growth habit is determined by recessive alleles at all three Vrn1 genes. Spring alleles at Vrn-A1, Vrn-B1, or Vrn-D1 were not detected for any entry, validating our assumption that this germplasm consists exclusively of winter wheat entries. There was variation among winter alleles at vrn-A1 and vrn-B1. The SNP in vrn-A1 associated with copy number [9] varied, with 264 (89%) of entries predicted to have three or more copies, four entries (1%) having two copies, and 29 entries (10%) having a single copy. Copy number variation corresponds with the winter alleles at Vrn-A1 described by Eagles et al. [46], such that the strong winter Wichita allele (Vrn-A1w) is associated with three or more copies of the gene and the 'Veery' [58] allele (Vrn-A1v) with two or fewer copies. Increased copy number at vrn-A1 has been associated with greater vernalization requirements, resulting in later flowering when the vernalization requirement is only partially fulfilled [16].
Two winter alleles at Vrn-B1 were previously described to affect heading date following a short vernalization duration, but with no effect when the vernalization requirement is fully satisfied [49]. The winter allele characteristic of 'AGS2000' [59] has an A at position 1752 and was associated with lower vernalization requirements and earlier heading than the 'NC-Neuse' [60] allele (C at position 1752). Low variation at this locus was observed in this germplasm. Most of the germplasm (291 entries, 99%) carried vrn-B1-Neuse; only 'TAM 401' [61] and TX05A001822 had vrn-B1-AGS2000.
The KASP assay identified the photoperiod insensitive allele at Ppd-D1. The PCR assay for Ppd-D1 further characterized allelic variation. Photoperiod insensitivity at Ppd-D1 is caused by a 2089 bp deletion upstream of the coding region [15]. In presence of the deletion, a 288 bp band is produced that corresponds with Ppd-D1a. The presence of the photoperiod sensitive Ppd-D1b allele is detected by amplifying a 414 bp fragment within the deletion region with primer pairs Ppd-D1_F and Ppd-D1_R1. Ppd-D1a was identified in 88 entries (29%), and the sensitive allele Ppd-D1b was found in 204 entries (68%). Guo et al. [40] described two small (24 and 15 bp) insertions within the intact 2089 bp region characteristic of an alternate photoperiod sensitive allele that results in a PCR fragment size of 453 bp. The alternate photoperiod insensitive allele was detected in seven entries (2%).
Among the 285 entries with complete allelic data at all five photoperiod and reduced-height loci, the most common allelic combination at the two reduced-height genes was Rht-B1b/Rht-D1a (203 entries, 71%), followed by Rht-B1a/Rht-D1a (63 entries, 22%), and Rht-B1a/Rht-D1b (14 entries, 5%, Table 3). The remaining entries were heterozygous at one or both reducedheight genes. No entries had mutations at both Rht-B1 and Rht-D1. These results are in agreement with those presented by Guedira et al. [23] in a similar set of germplasm.

Changes in Diversity of Photoperiod Alleles Over Time
The prevalence of the photoperiod sensitive alleles Ppd-A1b, Ppd-B1b, and Ppd-D1b has decreased gradually among U.S. Great Plains hard winter wheat germplasm over the past century ( Table 4). The first appearance of Ppd-A1a in this germplasm collection occurred with the derivation of 'TAM 301' [63] in 1991 (S1 Table). Ppd-A1a is rare in this germplasm but the prevalence of this allele has increased over time. Among entries derived after 1999, 4% were found to carry Ppd-A1a (Table 4). The photoperiod sensitive allele Ppd-B1b is found in 89% of entries derived prior to 1960, 67% of entries derived between 1960 and 1979, and 58% of those derived between 1980 and 1999. Among germplasm derived after 1999, Ppd-B1a is the more common allele and is found in 52% of entries. The prevalence of Ppd-D1b is similar among the 19 entries derived before 1960 (84%) as the 33 entries derived between 1960 and 1979 (85%). However, the percentage of entries carrying Ppd-D1b drops to 70% among entries derived between 1980 and 1999, and is further reduced to 64% among entries derived in or after 2000.

Effect of Photoperiod Alleles on Heading Date
Heading date had a grand mean of 130.8 days among all 299 entries and nine environments. In the combined analysis Ppd-A1 had a significant effect (P < 0.001) on heading date across environments. Entries with the photoperiod sensitive allele Ppd-A1b reached heading an average of 2.0 days later than those with Ppd-A1a (S3 Table). The effect of Ppd-A1 was not significantly different among environments, indicating a lack of significant genotype-by-environment (G×E) interaction at Ppd-A1. However, in analyses of individual environments the effect was significant (P < 0.05) only at Ar13R, where the effect size was 5.0 days. Only six entries had Ppd-A1a and the lack of significant effects in most environments is likely influenced by low diversity at this locus.
The effect of Ppd-B1 on heading date was significant in the combined analysis (P < 0.001). Ppd-B1b was associated with 3.0 days later heading across environments, but there was significant G×E interaction. Ppd-B1 had a significant effect (P < 0.01, S3 Table) in all individual environments. The effect of Ppd-B1b ranged from a minimum of 0.8 days later heading date in Fo13 (latitude of 40.65°N, Table 1) to a maximum of 5.2 days later in Ar13R (latitude of 34.18°N). There was a significant effect (P < 0.001) of Ppd-D1 on heading date in the combined analysis. The average effect of Ppd-D1b was 3.2 days later heading, which was larger than the effects of Ppd-B1b or Ppd-A1b. The effect of Ppd-D1 was significant (P < 0.05) in six environments, and the effect size ranged from a minimum of 0.9 day in Fo13 to a maximum of 5.2 days in Ar13R and It12R (S3 Table).
Models including Ppd-D1, Ppd-B1, and the interaction between Ppd-D1 and Ppd-B1 explained up to 65% of the variability in heading date (S4 Table). The interaction between Ppd-B1 and Ppd-D1 was significant in most environments (excluding Ha13R), but when the interaction was included the main effects of the loci generally became non-significant (S4 Table).
The interaction was such that entries with both photoperiod sensitive alleles reached heading much later than those with a single photoperiod sensitive allele or both photoperiod insensitive alleles. The effect of carrying both photoperiod sensitive alleles was greater than 6 days in four environments: Ar13R (8.6 days), It12R (8.0 days), Bu12R (6.7 days), and Ma12 (6.4 days). The effect of the significant interaction between Ppd-B1 and Ppd-D1 is illustrated for three Colorado environments (Gr12P, Gr12F, and Fo13) in Fig 1A-1C. The effect of interaction between Ppd-D1b and Ppd-B1b in the Colorado environments ranged from 1.5 to 4.1 days (S4 Table).

Effect of Semi-Dwarf Alleles on Heading Date
In the combined analysis, the semi-dwarf allele Rht-B1b had a significant (P < 0.001) effect on heading date and significant (P < 0.001) G×E interaction. Rht-B1b was associated with an average effect of 2.4 days earlier heading, but the effect ranged from 1.1 days in Ha13R to 4.2 days in Ar13R (S3 Table). There was a low level of genetic diversity at Rht-D1 and the locus did not have a significant effect on heading date when it was the only genetic term in the model.
There were 289 entries with homozygous allele calls at Rht-B1 and Rht-D1. When both loci were included in the model, the main effects of both Rht-B1 and Rht-D1 were significant (P < 0.001, Table 5). There was significant G×E interaction for Rht-B1 but not Rht-D1. Across environments Rht-D1b had an effect of 2.9 days earlier heading. The effect of Rht-B1b was significant (P < 0.05) in Ar13R, Bu12R, Gr12P, and It12R and ranged from a minimum of 2.6 days earlier heading in Gr12P to a maximum of 4.7 days earlier in Ar13R (Table 6). Interaction between Rht-B1 and Rht-D1 could not be tested because not all four allelic combinations were present in this germplasm.

Effect of Vernalization Alleles on Heading Date
No entry had alleles for spring growth habit detected at Vrn-A1, Vrn-B1, or Vrn-D1. The effect of different vrn-B1 winter alleles was significant (P <0.01) in the combined analysis, with the 'Neuse' allele heading 2.2 days later (S3 Table). However, only two entries carried the vrn-B1 'AGS2000' allele. The G×E term was not significant for vrn-B1.
Variation observed for the winter allele vrn-A1 had a significant effect (P < 0.001) on heading date that varied among environments (S3 Table). The differential effect of three or more vrn-A1 copies relative to a single copy ranged from a minimum of 0.5 days later heading in Gr12P to a maximum of 4.7 days later heading in Ar13R. The effect of three or more copies of vrn-A1 relative to two copies also had a significant (P < 0.01) effect in the combined analysis, such that entries with two copies reached heading 1.9 days earlier. This effect had significant G×E interaction, and in individual environments the effects of three or more copies relative to two copies were only significant in Ar13R, where the effect size was 5.5 days (S3 Table).
Copy number variation at Vrn-A1 has been previously shown to have a positive association with vernalization requirements and survival under freezing conditions [46]. The strong effect of vrn-A1 in Ar13R suggests winter temperatures at this environment may not have fully satisfied the vernalization requirements of all entries. Ar13R experienced a mild winter with daily highs always above 0°C. The effect of allelic variation at Vrn-A1 did not have a significant effect at It12R, which experienced a cold winter with 11 days that had highs below 0°C and 65 days with lows below 0°C.

Multivariate Analyses of Alleles on Heading Date
Combined analyses across environments. The significance of Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, Rht-D1, and their interactions were tested in a multiple-gene model using a subset of 280 entries with complete genotypic data and homozygous calls at all loci (Table 3). For the combined analysis the full model included main effects of the environment, main effects of all five genes, all 10 two-way interactions between genes, and the interaction of each genetic term with the environment. Model selection was used to identify the best-fit model based on lowest AIC value.
The best-fit model to explain heading date across environments included effects of the environment; major effects of all five candidate loci (Ppd-D1, Ppd-B1, Ppd-A1, Rht-B1, and Rht-D1); the interaction between Ppd-D1 and Ppd-B1; the interaction between Ppd-B1 and Rht-B1; two-way interactions between the environment and Ppd-D1, Ppd-B1, and Rht-B1; and a three-  way interaction between the environment, Ppd-D1, and Ppd-B1 (Table 7). This model explained 96.2% of the phenotypic variation in heading date, with most of the variation (92.2% of sums of squares) explained by the environment, and a smaller amount (4.2%) explained by genetic terms. Of the variation explained by genetic terms, 30.1% was due to the main effect of Ppd-B1, 29.6% was due to Ppd-D1, and 13.6% was due to their interaction. Smaller amounts of the genetic variation were explained by the two-way interactions between Ppd-D1 (9.1%) or Ppd-B1 (8.7%) and the environment, the three-way interaction between Ppd-D1, Ppd-B1, and the environment (4.1%), the main effect of Rht-B1 (2.7%), and the interaction between Rht-B1 and Ppd-B1 (1.2%). Less than 1% of the genetic variation was due to Ppd-A1, Rht-D1, the twoway interaction between Ppd-B1 and Ppd-A1, or the two-way interaction between Rht-B1 and the environment. The relative importance of each locus in explaining variation in heading date was also evaluated through model averaging using Akaike weights. The Akaike weights were 1.00 for Ppd-B1, Ppd-D1, and Rht-B1; 0.96 for Ppd-A1; and 0.71 for Rht-D1. Individual analyses for each environment. Multi-gene models were fit to each environment separately using the same model selection approach as in the combined analysis. The full model included the main effects and all two-way interactions between five genes: Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, and Rht-D1. The best-fit models varied among environments and are summarized in Table 7, and not all genetic terms are significant. The effects of Ppd-D1, Ppd-B1, and the interaction between Ppd-D1 and Ppd-B1 were included in the best-fit model for each environment. The interaction between Ppd-D1 and Ppd-B1 had a maximum effect size of 8.5 days in Ar13R. Rht-B1 was also included for every environment except Gr12F. None of the best-fit models included the main effect of, or interactions withPpd-A1 or Rht-D1. Some bestfit models also included interactions between Ppd-B1 and Rht-B1 or Ppd-D1 and Rht-B1. The combined allelic effects explained an average of 44% of the variability in heading date, but ranged from a minimum of 14.7% in Fo13 to a maximum of 69.5% in Ma12 (Table 8). Table 6. Allelic effects (number of days) of Rht-B1b in a model that estimated days after 1 January to heading in 285 winter wheat entries grown in nine environments. The environments are described in Table 1. The model terms were fit separately for each environment. The model effects included environment, Rht-B1, Rht-D1, and Rht-B1-by-environment interaction. The intercept (Int) describes the number of days from 1 January to heading in each environment before the allelic effects are applied. The allelic effect (number of days) at each locus is added to the Int value. The allelic effect of Rht-D1b was -2.68 days and did not have significant interaction with environment.   Table 8. Allelic effects (number of days) of gene-based terms included in the best-fit model for winter wheat heading date in each of nine environments, and the proportion of variability (R 2 ) in heading date explained by all terms in each model. The environments are described in Table 1. The model terms were fit separately for each environment. The intercept (Int) describes the number of days from 1 January to heading in each environment before the allelic effects are applied. The allelic effect (number of days) at each locus is added to the Int value. temperatures than northern locations, and Montana and North Dakota have greater changes in day length throughout the year than the central and southern plains, due to their northern latitudes.

Env
In a combined analysis across environments there were significant differences (P < 0.001) in heading date among entries from different states of origin. However, not all pairs of states had significantly different heading dates. The entries were divided into three broad regions within the U.S. Great Plains: northern plains (Montana, North Dakota, and South Dakota; 39 entries), southern plains (Texas and Oklahoma; 105 entries), and central plains (Colorado, Kansas, and Nebraska; 119 entries, Fig 2A). There were significant differences (P < 0.01) in heading dates among entries originating from each pair of regions. Average heading dates were earliest among entries derived in the southern plains (128.8 ± 14.3 days after 1 January). Heading occurred about 2.8 days later in entries derived in the central plains (131.6 ± 13.1 days) compared to entries from the southern plains. Latest heading dates were observed for germplasm from the northern plains (134.9 ± 12.2 days).
All three photoperiod insensitive alleles (Ppd-A1a, Ppd-B1a, Ppd-D1a) were present at higher levels in germplasm from the southern plains than those from the central or northern plains (Fig 2B). The photoperiod insensitive allele Ppd-A1a was uncommon in this germplasm. Ppd-A1a was only present in six (2%) of the 263 entries with known geographic origins, all of which originated in the southern plains. Ppd-B1a was much more evenly distributed-present in 40% of entries overall-but varied gradually among the three regions. Ppd-B1a was found at highest proportions (57%) in germplasm from the southern plains, and at much lower levels in the central (34%) or northern plains (13%). Similar patterns were observed for Ppd-D1a, which was found in 44% of entries from the southern plains, 19% of entries from the central plains, and only 10% of entries from the northern plains.
The semi-dwarf alleles at the reduced-height loci were more common in germplasm derived in the southern plains than those originating from the northern or central plains. Rht-B1b was found in 71% of entries overall, including 83% of entries from the southern plains, 63% of entries from the central plains, and only 51% of entries from the northern plains ( Fig 2B). The semi-dwarf allele Rht-D1b was rare in the panel, and present in only 11 (4%) entries. However, the distribution of this allele was much higher in the southern plains (8 entries, 8%) than the northern plains (1 entry, 3%) or central plains (2 entries, 2%).

Discussion
Our first objective was to characterize allelic variants at loci known to affect developmental timing of wheat through the vernalization, photoperiod sensitivity, and reduced-height pathways. Marker analysis revealed the incidence of semi-dwarf alleles Rht-B1b and Rht-D1b in this germplasm was comparable to that reported by Guedira et al. [23] in a similar set of germplasm. The photoperiod insensitive alleles Ppd-B1a and Ppd-D1a were present at similar levels as detected by Kiss et al. [24] in a worldwide collection of winter wheat. We found Ppd-A1a to be rare in this germplasm, and present at slightly lower levels than reported among European accessions [26]. Further allelic variation could exist in our panel as additional allelic variants [40,50,64] or copy number variants [25,65] not detected by our assays.
Guo et al. [40] described four modern and two ancient Ppd-D1 haplotypes present in a worldwide collection of wheat, and detected all four modern haplotypes in accessions from the U.S. and Canada. Our analyses could only differentiate two modern and one ancient haplotype, which suggests additional haplotypes could be present. Seven of our entries had an alternate photoperiod sensitivity allele at Ppd-D1b, described by Guo et al. [40] as haplotype VI. The alternate allele did not have a significant effect on heading date across environments, probably due to the small sample size. The entries with haplotype VI are 'OK Rising' (PI 656382), 'Thunder CL' [66], TAM 303, KS00F5-20-3, 'Overley' (PI 634974), 'Chisholm' [67], and 'Custer' (OK88767-11). Haplotype VI is thought to be an ancient deletion and is associated with Aegilops tauschii accessions and synthetic hexaploid wheats [40]. All seven entries we characterized as haplotype VI have pedigree contributions from Aegilops tauschii or Eastern European Effects of Major Genes on Wheat Heading Date germplasm. The Eastern European parentage might also include Aegilops tauschii, as exotic germplasm and interspecific crosses were introduced to this region during the 1970s and 1980s [68].
Additional alleles have also been described at Ppd-A1 [50]. The photoperiod insensitive allele Ppd-A1a was expected to be associated with earlier heading, but did not have a significant effect on heading date on this germplasm. We only characterized six entries as carrying Ppd-A1a, all of which originated from the southern plains. This region has less variation in day length than the other U.S. Great Plains regions, so photoperiod sensitivity is expected to have less of an effect, and photoperiod insensitivity to be more common, than at more northern latitudes. Guedira et al. [49] found the photoperiod insensitive allele Ppd-A1a.1-first identified in Japanese germplasm-to be common among winter wheat varieties originating in the eastern U.S. [50]. It is likely that the effect of Ppd-A1a on heading date would be significant among germplasm with greater diversity at this locus, such as a panel that included more entries from southern or eastern regions of the United States.
Our second objective was to estimate the allelic effects of known developmental genes on heading date. We found the timing of heading date in winter wheat from the U.S. Great Plains to be more strongly affected by photoperiod loci than the vernalization or reduced-height loci. Alleles of these five developmental genes (Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, and Rht-D1) had the greatest genetic effects on variation in heading date (Table 7). However, terms included in the best-fit model varied slightly among environments (Table 8). Additionally, the effects of Ppd-A1 and Rht-D1 were not significant in any individual environment, and have relatively low Akaike weights, which indicates a lower level of likelihood these loci are affecting heading date. Differential effects among environments suggest the model might be improved if the environmental effects were deconstructed into multiple environmental variables (such as effects of temperature, moisture, day length, etc.). Furthermore, additional variation in number of days after 1 January to heading might be improved through genotyping and inclusion of other loci known to affect developmental timing, such as Vrn-B3 and Vrn-D3 [11,19,45].
In our study, Rht-B1 was found to have a significant effect on heading date in most environments (Table 8). While it is unclear how reduced sensitivity to GA affects heading date, Rht-B1 has been previously shown to have a small effect on heading date [20]. Wilhelm et al. [69] suggest one possibility is the role of other genes that are tightly linked with Rht-B1, such as Teosinte branched 1 (TaTb1), which is associated with tillering and fertility. Recently, Li. et al. [64] reported numerous allelic variants at Rht1 (Rht-A1, Rht-B1, and Rht-D1), which were largely unaccounted for in our analyses. Among these variants were novel alleles that increased plant height relative to the tall alleles Rht-B1a or Rht-D1a.
Most of the genetic effects were due to the main effects and interaction of Ppd-B1 and Ppd-D1, and their three-way interaction with the environment. The interaction between Ppd-B1 and Ppd-D1 was such that entries carrying both photoperiod sensitive alleles reached heading much later than those carrying a single photoperiod sensitive allele or both insensitive alleles (Fig 1, S4 Table, Table 8). The presence of both photoperiod sensitive alleles delayed heading date by an average of 4.8 days across all environments, and a maximum of 8.5 days in Ar13R (Table 8). While the genetic effects on heading date are small relative to the environmental effects (Table 7), even minor variation in the timing of heading date-variation of several days -may allow for more precise targeting of wheat varieties to different regions, especially if climatic patterns shift within established growing regions [70], or wheat cultivation expands to new regions [71].
The allelic effects of Ppd-B1 and Ppd-D1 varied among environments, suggesting differential expression under varying environmental conditions such as temperature, or moisture. The greatest significant effects-more than five days at each locus-were observed at Ar13R (S3 Table). The effect of the interaction between these loci was 8.6 days at Ar13R (S4 Table). Variation in the effect size and level of significance does not follow a discernable trend among environments; however, the G×E term is likely affected by different sources of environmental variability, such as climatic conditions, management practices including planting and harvest dates, or biotic and abiotic stresses. For instance, heavy rainfall during May and June, and warm maximum daily temperatures throughout the winter, spring, and summer have influenced the large allelic effects at Ar13R. Interaction between Ppd-D1 and Ppd-B1 has been previously characterized [72], is known to vary with CNV at Ppd-B1 [46], and is known to have much larger effects in spring than winter wheat [24].
Allelic variation associated with CNV at Vrn-A1 had a significant effect on heading date in this germplasm (S3 Table). Low copy number at Vrn-A1 is associated with earlier flowering following a short vernalization period [16]. Copy number variation at Vrn-A1 was inferred by SNPs on exons four and seven. The SNP on exon four is associated with having two or fewer copies of the gene, and was used by Eagles et al. [46] to distinguish Vrn-A1v associated with earlier heading from Vrn-A1w associated with greater freezing tolerance. We found the effect of CNV at Vrn-A1 to vary among environments with contrasting winter temperatures. The effect of CNV at vrn-A1 was largest (4.7 days) at Ar13R (S3 Table), and that environment experienced a mild winter that had daily highs above 0°C and only three days with lows below 0°C. The second largest significant effect was seen at Bu12R, which is also located in a mild, southern environment. By contrast, variation at Vrn-A1 did not have a significant effect at It12R, which experienced a cold winter with 65 days-including 29 continuous days-that had lows below 0°C.
Our third objective was to investigate temporal and geographic patterns of allelic distribution in the U.S. Great Plains. Changes in allele frequencies at the photoperiod genes indicate selection for photoperiod insensitivity during the last century of winter wheat breeding (Table 4). This could be the direct result of selecting for improved adaptability based on flowering time [73]. Most strikingly, the proportion of entries carrying the Ppd-B1b sensitive allele has gradually declined over time, from 89% of entries derived before 1960 to only 48% of germplasm derived after 2000 (Table 4). An increase in the number of varieties carrying the photoperiod insensitive allele suggests a contribution to greater adaptation to specific environments.
Selection against photoperiod sensitivity is also apparent in the distribution of alleles from breeding programs in different sub-regions of the U.S. Great Plains. We saw the incidence of photoperiod sensitive alleles decrease from northern to southern latitudes (Table 8). Various environmental variables, such as temperature and day length, vary with latitude. Extreme temperature events such as freezing temperatures and heat stress can have very pronounced and detrimental effects on yield, especially during vulnerable developmental stages. Therefore, targeting heading date could reduce the chance that extreme temperatures or stress damage the developing reproductive structures.

Conclusions
In this study we used a panel of 299 hard winter wheat entries representative of modern and historic U.S. Great Plains germplasm to evaluate allelic diversity and effects of vernalization, photoperiod, and reduced-height loci on the timing of heading date. We found most of the genetic effects of heading date to be explained by Ppd-B1, Ppd-D1, and the interaction between these loci. A smaller amount of variation was explained by Rht-B1. Across nine environments, the interaction between Rht-B1 and Ppd-B1 also had a small but significant effect on heading date. Both photoperiod sensitive and insensitive alleles were common for Ppd-B1 and Ppd-D1, and an alternate photoperiod sensitivity allele associated with ancestral wheats was detected at Ppd-D1. There was limited allelic diversity at Ppd-A1 and Rht-D1, and these loci did not have a significant effect on heading date. The presence of photoperiod sensitive alleles Ppd-A1b, Ppd-B1b, and Ppd-D1b has been decreasing over the past century of wheat breeding, and these alleles are less common in the southern sub-region of the U.S. Great Plains than either the central or northern plains. Our analyses enhance the understanding of the effects of developmental genes on winter wheat under varying environmental conditions. This information can potentially be used to improve breeding strategies for current and future climate scenarios.
Supporting Information S1 Table. Description of U.S. winter wheat entries evaluated, including entry name, year derived, state breeding program and region of the U.S. Great Plains in which derived, and allele calls. The table is sorted first alphabetically by region, then temporally by year, and finally alphabetically by entry name. The genotyped loci include photoperiod (Ppd) genes where alleles are 'a' insensitive and 'b' sensitive, reduced-height (Rht) genes where alleles are 'a' tall and 'b' semi-dwarf, and vernalization (Vrn) genes. Vrn-A1, Vrn-B1, and Vrn-D1 allele calls are 'W' for winter growth habit. Variation in winter alleles at vernalization loci include copy number variation (CNV) for vrn-A1, Wichita-type ('w') or Veery-type ('v') alleles for vrn-A1, and Neuse-type ('N') or AGS2000-type ('A') alleles for vrn-B1. Missing genotypic data are indicated with a dash (-). (DOCX) S2 Table. Description of KASP markers used to genotype 299 U.S. Great Plains hard winter wheat entries. KASP detected allelic variants at Vrn-A1, Vrn-B1, Vrn-D1, Ppd-A1, Ppd-B1, Ppd-D1, Rht-B1, and Rht-D1. (DOCX) S3 Table. Allelic effects (number of days) of photoperiod, reduced-height, and vernalization loci on heading date in each of nine environments. The environments are described in Table 1. The model terms were fit separately for each locus. The intercept (Int) describes the number of days from1 January to heading in each environment before the allelic effect is applied. The allelic effect at each locus is added to the Int value. Allelic effects were fit separately for each environment when there was significant genotype-by-environment (G×E) interaction at that locus in the combined analysis across all environments. Allelic effects from the combined analyses are reported for Ppd-A1 and vrn-B1 because significant G×E interaction was not observed at these loci. The grand mean heading date across all germplasm and environments was 131.8 ± 0.3 days. (DOCX) S4 Table. Allelic effects (number of days) of photoperiod sensitive alleles Ppd-D1b and Ppd-B1b. Allelic effects and interaction of Ppd-D1b and Ppd-B1b on winter wheat heading date in each of nine environments, and proportion of variability (R 2 ) in heading date explained by all terms in each model. The environments are described in Table 1. The model terms were fit separately for each environment. The intercept (Int) describes the number of days from 1 January to heading in each environment before the allelic effects are applied. The allelic effect (number of days) at each locus is added to the Int value. (DOCX)