Quantification of the effects of VRN1 and Ppd-D1 to predict spring wheat (Triticum aestivum) heading time across diverse environments.

Heading time is a major determinant of the adaptation of wheat to different environments, and is critical in minimizing risks of frost, heat, and drought on reproductive development. Given that major developmental genes are known in wheat, a process-based model, APSIM, was modified to incorporate gene effects into estimation of heading time, while minimizing degradation in the predictive capability of the model. Model parameters describing environment responses were replaced with functions of the number of winter and photoperiod (PPD)-sensitive alleles at the three VRN1 loci and the Ppd-D1 locus, respectively. Two years of vernalization and PPD trials of 210 lines (spring wheats) at a single location were used to estimate the effects of the VRN1 and Ppd-D1 alleles, with validation against 190 trials (~4400 observations) across the Australian wheatbelt. Compared with spring genotypes, winter genotypes for Vrn-A1 (i.e. with two winter alleles) had a delay of 76.8 degree days (°Cd) in time to heading, which was double the effect of the Vrn-B1 or Vrn-D1 winter genotypes. Of the three VRN1 loci, winter alleles at Vrn-B1 had the strongest interaction with PPD, delaying heading time by 99.0 °Cd under long days. The gene-based model had root mean square error of 3.2 and 4.3 d for calibration and validation datasets, respectively. Virtual genotypes were created to examine heading time in comparison with frost and heat events and showed that new longer-season varieties could be heading later (with potential increased yield) when sown early in season. This gene-based model allows breeders to consider how to target gene combinations to current and future production environments using parameters determined from a small set of phenotyping treatments.

Introduction

The timings of heading (head emergence from leaf whorl) and flowering (first anther burst on spikes) are primary determinants of wheat adaptation to diverse environments to avoid frost, heat, and terminal drought stresses, and are influenced by three groups of genes responding to average temperature [earliness per se (EPS)], cold temperature [vernalization (VRN)], and day length [photoperiod (PPD)] (Rousset ). Although many crop models have been parameterized to predict wheat phenology (Jones ; Keating ), wheat breeding and agronomy research would benefit from robust models that quantitatively predict the interaction of genotype and environment based on known genes (White ).

Conventional process-based ecophysiological models simulate wheat development and growth through consideration of genotype, environment, and management effects (G×E×M) (Jones ; Keating ). Phenological stages are simulated in response to temperature (development rate and VRN requirement), PPD (day length), and the interaction effects of these inputs. To simulate the phenology of different genotypes, parameters of these models (Bertin ) are typically calibrated and validated using field trials from multiple sites and years (White ). However, parameter estimation of process-based models is time-consuming and expensive when considering large numbers of genotypes (Bertin ), so new genotypes are slow to be characterized into models.

Genes controlling development have been identified and cloned in multiple species, and researchers have linked parameters of process-based models to quantitative trait loci (QTL) or known genes to try to explain variation in genotype responses (gene-based models) (Bertin ). Yin , 2000) introduced a gene-based model to predict barley flowering time, specific leaf area, and yield from detected QTLs in a bi-parental population. Similar QTL approaches have been used to model maize leaf elongation (Chenu ), soybean phenology (Stewart ), peach fruit quality (Quilot ), barley phenology (Yin ), rice phenology (Nakagawa ), and Brassica oleracea phenology (Uptmoor ). White and Hoogenboom (1996) and Hoogenboom and White (2003) developed a gene-based model, GeneGro, to predict phenology, growth habit, and seed size of soybean (e.g. Ppd, Hr). In this approach, homozygous genotypes with two dominant alleles are counted as 1 and recessive genotypes as 0. The integer number of dominant alleles is summed across loci for a specific trait and then fitted linearly with line-specific parameter values of process-based models. Similar approaches have been extended to soybean development and yield (Messina ) and to wheat phenology (White ). Simulations of gene networks associated with flowering time have been undertaken by Chapman for sorghum and by Welch and Salazar for Arabidopsis, considering the major gene pathways for flowering time. Wenden also presented a computational model based on signal and genes to predict pea flowering time.

The accuracy of genetically parameterized development models is typically reduced when genotypic (line-specific) parameters are replaced by generic (gene-specific) parameters (White ; Uptmoor ). Line-specific parameters vary with genotype, but gene-specific parameters are similar among genotypes, and even among species (Bertin ). A potential source of error is firstly to fit line-specific parameters of process-based models, and then re-estimate the new parameters associated with these QTL/genes. Furthermore, most genetically parameterized models consider only additive effects (Messina ; White ). However, even the known genes in wheat interact in complex ways. In order to flower, wheat plants (including completely ‘spring’ wheats) generally require a VRN event, followed by exposure to a lengthening PPD. The VRN requirement of wheat is determined mainly by three homoeologous VRN1 genes, Vrn-A1, Vrn-B1, and Vrn-D1 (Rousset ). In winter wheats, under long days, PHOTOPERIOD1 (PPD1) can promote flowering time via effects of a CCT domain promoter, CO2 (a wheat homologue of constans), releasing the expression of VRN3 and then VRN1 (Li ).

The aim of this study was to build a gene-based model to predict spring wheat heading time from allele combinations of VRN1 (Vrn-A1, Vrn-B1, and Vrn-D1) and PPD1 (Ppd-D1 only) genes, using data from field experiments executed at one location. Gene effects were estimated from data collected under VRN and PPD treatments, and an optimization in one step was used to estimate the final model parameters. Observations of heading time from trials across the Australian wheatbelt were used to validate the gene-based model. As an example of applying the model, for a range of representative locations, virtual genotypes were created to indicate the gene combinations that would minimize exposure to heading time frost and heat risk across the Australian wheatbelt.

Materials and methods

Sources of data

To derive parameters for calibration of the model, we used two years of trials in one location with treatments of VRN and PPD applied to 210 wheat lines that had been genotyped for known genes. The model was validated on multiple datasets comprising a large range of sowing times and latitudes across the Australian wheatbelt for 172 out of the 210 wheat lines.

Calibration dataset

The VRN and PPD trials comprising four treatments were conducted in 2008 and 2009 in South Perth, Western Australia (31.99°S, 115.88°E). The VRN treatments were either a natural VRN (V1) or a pre-imbibed seed VRN for 8 weeks at 4 °C (V2). The PPD treatments were natural day length (from 9 to 11h, P1) or extended day length to midnight (day length of 17–18h, P2; Fig. 1). The day length was extended with OSRAM L36W/850 (3300 luminous flux) lights mounted 1.5 m above the plants.

Fig. 1.

Daily minimum (MinT) and maximum (MaxT) temperature and day length for natural and extended treatments during growth season. The day length was extended to midnight using lights (3300 luminous flux).

Each plot consisted of a single row of seeds sown in 17cm rows on 11 June in both years, so as to emerge just prior to the shortest day of the year (i.e. 20–21 June) when day length at the site is 8.8h. Heading times were recorded when 50% of heads were fully emerged from the flag leaf. There were three replicates for 36 lines, but only one replicate for the remaining lines. The lines included all major spring wheat varieties grown in Australia over the last 25 years, as well as 86 pre-release experimental lines supplied by wheat breeding companies.

The hourly temperature at 2 m screen height was recorded by on-site weather stations. There were two data loggers for treatments of natural and extended PPD in 2008, respectively, but only one data logger in 2009. There were no differences in temperature between treatments of natural and extended PPD (data not shown).

Validation datasets

The validation datasets had been collected previously by multiple researchers from 79 sites across the Australian wheatbelt with a wide range of environments (from 17.2°S to 38.1°S; from 114.7 E to 152.3°E; Fig. 2) over 7 years (from 2005 to 2011). These trials were designed to evaluate adaptation of wheat cultivars and the impact of sowing time on heading, and comprised 68 different sowing dates extending from 30 March to 11 September. The validation datasets contained 172 of the lines that had been grown in the calibration datasets. The heading time had been determined as the date when 50% of plants in a plot had reached heading (Zadoks stage 55).

Fig. 2.

The locations of calibration (blue) and validation (red) datasets with circle size indicating the number of observations in 190 trials across the Australian wheatbelt (grey shaded region).

Hourly weather records for Katanning, Northam, Roseworthy, and Gatton were obtained from on-site loggers (2 m above ground; TinyTag Gemini loggers). Weather records for other sites were obtained from the nearest Bureau of Meteorology sites (http://www.longpaddock.qld.gov.au/silo/). The maximum distance to the closest weather station was 31 km, and the median distance was 7 km, with little difference in elevation. The topography in the areas of these trials is typically flat with consequently little effect on daily temperature values. See Supplementary Materials and Methods at JXB online for additional detail about validation datasets.

Genotyping of the VRN1 and Ppd-D1 genes

At the time of sowing in the 2008 South Perth trials, DNA was extracted from ten seedlings for the determination of homozygotes at four loci in each of the 210 lines. Homozygotes for Vrn-A1, based on the variation identified in the promoter region of the Vrn-A1 locus, were classified using the protocol described by Yan . Homozygotes of both Vrn-B1 and Vrn-D1 were identified as described by Fu . Homozygotes of Ppd-D1 were identified using the protocol described by Beales . The resolution of PCR products was found to be vastly improved by the addition of 1% polyvinylpyrrolidone to the extraction buffer to purify the DNA sample prior to PCR. In the designation of the alleles observed, we adopted the allele classification utilized by Eagles , with a and v for spring and winter alleles/homozygotes of VRN1 genes, respectively, and a and b for insensitive and sensitive alleles/homozygotes of Ppd-D1, respectively.

In contrast to spring or facultative wheats, true winter wheats carry winter VRN alleles at VRN2 and/or VRN3 loci (Yan ). Based on pedigree, all of the lines in our experiments were expected to carry spring alleles at the other VRN loci, VRN2 and VRN3, which are understood to be downstream in the flowering pathway of the VRN1 loci (Yan ), i.e. all of the lines grown here were ‘spring’ wheats with no effect of VRN2 or VRN3 loci.

The original phenology model of APSIM-Wheat (APSIM-Wheat-O)

APSIM (Agricultural Production Systems Simulator) is a daily time-step cropping system model that simulates soil water, residue, and nutrient dynamics, and the growth and development of more than 30 crops (Keating , http://www.apsim.info).

In the release version 7.4 of APSIM (designated ‘APSIM-Wheat-O’ here), the growing season between sowing and harvest is divided into nine phases, each of a fixed thermal time (TT) (target TT). The commencement of each phase is determined by the accumulation of TT, except from sowing to germination, which is driven by soil water content and sowing depth. On the day when the accumulated TT from the start of a phase exceeds the target TT of that phase, the phase is completed, and the ‘excess’ TT on that day is credited towards the target in the next phase.

The accumulated TT () towards the target TT is a sum of daily TT (ΔTT) modified by genetic and environmental factors:

ΔTT is calculated from daily mean temperature using three cardinal temperatures (0, 26, and 34°C for base, optimum and maximum temperature, respectively): when daily mean temperature is <26 °C,ΔTT is equal to the daily mean. Genetic effects include a VRN factor (F ) and PPD factor (F ) calculated daily (see below) and equal to 1 before emergence and after floral initiation. Environmental factors include soil water (F ), nitrogen (F ) and phosphorus stresses (F ). In the current release (APSIM 7.4, as at 1 December 2012), these environmental factors are set to a default of 1.0 and have no influence on wheat phenology.

Vernalization

Daily VRN in APSIM-Wheat-O (ΔV) is simulated from daily mean crown temperature (T ), daily maximum air temperature (T ) and minimum air temperature (T ) with crown temperature being daily mean temperature adjusted by snow depth (Ritchie ).

Devernalization (ΔV ) can occur if daily T is above 30 °C and the total VRN (V) is less than 10:

The total VRN (V) is calculated as:

The VRN factor (F ) is calculated from plant emergence to floral initiation and is updated daily.

where R is the sensitivity to VRN.

Photoperiod

The PPD factor (F ) is calculated by

where L is the day length plus civil twilight (h) (i.e. the centre of the Sun’s disc is 6° below the horizon), and R is the sensitivity to PPD.

Both R and R need to be determined for each new genotype added to APSIM. For the purposes of this paper, the entire phenology algorithm from APSIM was converted to a program in the R language (R Development Core Team, 2008), and modifications of the algorithm were compared using R scripts implemented across a high-performance computer system.

A modified phenology model of APSIM-Wheat (APSIM-Wheat-M)

Prior to the introduction of the gene-based parameters, two changes were made to the calculation of effects of VRN and PPD.

By definition, the minimum function in Equation 1 of APSIM-Wheat-O considers only the maximum effect of either VRN or PPD on the accumulation of daily TT. However, a multiplicative function to accommodate observed interactions between VRN and PPD effects has been used in other models, e.g. ARCWHEAT 1 (Weir ) and Sirius (Jamieson ). In order to capture the observed effects of interactions, Equation 1 was changed to:

Slafer and Rawson (1994) indicated that the developmental rate of wheat is sensitive to PPD until flowering and so the second change was to extend the effect of PPD to be applied between the floral initiation and flowering stage.

A gene-based model of APSIM-Wheat (APSIM-Wheat-G)

In the gene-based model (‘APSIM-Wheat-G’, modified from APSIM-Wheat-M), the VRN (R ) and PPD (R ) sensitivities in Equations 5 and 6 were related to the number of sensitive alleles of the VRN1 and Ppd-D1 genes. Previous studies have found different effects among VRN1 and Ppd-D1 genes (González ; Loukoianov ; Eagles ; Allard ), so the different effects were allowed to vary in magnitude via a weighting function. A multiplicative model was used to simulate the interaction between VRN1 and Ppd-D1 genes on PPD sensitivity.

Each locus of the VRN1 and Ppd-D1 genes was associated with a genotype value of 0 for the spring or PPD-insensitive homozygote at that locus and 1 for the winter or PPD-sensitive homozygote at that locus. No heterozygous lines were present in the dataset. The total weighted numbers of VRN1 (N ) and Ppd-D1 (N ) genes were calculated by weighting and summing the genotype values of 0 or 1 at each VRN1 and Ppd-D1 locus:

where h , h , and h were a weighting for the effect at each VRN1 locus on VRN sensitivity, and h , h , and h were the weighted VRN1 effects on PPD requirement.

Linear functions were used to simulate the relationships between the weighted numbers of VRN1 or Ppd-D1 alleles and the sensitivity of target processes:

where k and k were slopes of lines for VRN and PPD, respectively, and b and b were intercepts of lines to indicate the unknown effects of the VRN and PPD genes, respectively.

Estimation of parameter values for APSIM-Wheat-O and APSIM-Wheat-M

In APSIM-Wheat-O, the VRN sensitivity (R ), PPD sensitivity (R ) and target TT from floral initiation to flowering (TT ) are three line-specific parameters to determine the heading time of wheat. The heading time is estimated as Zadok’s stage 55 in APSIM. Over the last 20 years or so, the phenology model had been calibrated for 103 named varieties in APSIM release 7.4, including 40 of the lines in our datasets. Perhaps due to fitting to ‘local data’, APSIM-Wheat-O gave poor predictions (R 2 <0.7) for at least nine of the 40 lines in our calibration dataset (RMSE=7.2 d), so we refitted the three line-specific parameters (R , R , and TT ) of all 210 lines for APSIM-Wheat-O (Equation 1) and APSIM-Wheat-M (Equation 7). Using both the calibration and validation datasets together to cover a wide range of environments and maximize model performance, the three line-specific parameters in APSIM-Wheat-O and APSIM-Wheat-M (610 parameters in total per model) were fitted to minimize RMSE between observed and predicted heading times.

Estimation of parameter values for APSIM-Wheat-G

Using the observed calibration data from the 2 years of trials at South Perth, the effects of the VRN1 and Ppd-D1 genes were estimated using the restricted maximum likelihood (REML) method in ASReml (version 3, http://www.vsni.co.uk/software/asreml). We first calculated the heading time differences for each line and year by subtracting the mean values in the V2P2 treatment from each of the other treatments. A mixed model (Equation 12) was used to estimate the heading time differences (fixed effects) for each multi-locus genotype (MLG , Table 1) for the differences V2P1 – V2P2, V1P2 – V2P2, and V1P1 – V2P2, separately, while the other components of years, lines, and residuals were considered random effects.

Table 1.

Frequency distribution of MLGs for calibration and validation datasetsThe MLGs are the alleles of Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1. a and v indicate homozygous genotypes for spring and winter alleles of VRN1 genes, respectively, while a and b indicate homozygous genotypes for the insensitive and sensitive alleles of Ppd-D1, respectively.

MLGs	Calibration	Validation
a-a-a-a	8	7
a-a-a-b	1	1
a-a-v-a	14	13
a-a-v-b	21	17
a-v-a-a	5	5
a-v-v-a	39	33
a-v-v-b	11	8
v-a-a-a	26	24
v-a-a-b	4	4
v-a-v-a	32	32
v-a-v-b	21	12
v-v-a-a	11	9
v-v-v-a	14	6
v-v-v-b	3	1

where H is the heading time difference for the ith MLG in the jth year and kth line, μ is the overall mean, M is the ith MLG effect (fixed), y is the jth year effect (fixed), L is the kth line effect (random) and ε are the residuals.

Considering two alleles at each of the Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1 loci results in 16 MLGs. In the pre-VRN and extended-PPD treatment (V2P2), we assumed that the VRN requirement was satisfied by 8 weeks of cold treatment at 4 °C, and that PPD sensitivity was overcome by extending day length to 17–18h. This ‘control’ treatment was used to estimate EPS for each line, i.e. the shortest possible TT to flowering. The differences between VRN treatments were used to estimate the effects of VRN1 genes (i.e. V1P2 – V2P2), and the differences between PPD treatments to estimate effects of Ppd-D1 (i.e. V2P1 – V2P2).

The final gene-based model included the four gene-specific parameters (k , k , b , and b ) and one line-specific parameter (TT ). These parameters were fitted simultaneously for 210 lines (214 parameters in total). We generated exhaustive combinations of these parameters across a wide range (from 0 to 3 for k , k , b , and b at 0.01 intervals; from 300 to 1300 degree days (°Cd) for TT at 5 °Cd intervals), and then calculated heading times for the eight calibration environments and lines (>5 billion estimates). The parameter values were selected according to the minimum RMSE between the observed and predicted heading times of all treatments in the calibration dataset.

Virtual genotypes

To account for local variability in frosts and heat (damaging at heading) and in seasonal rainfall, the range of heading times at a specific sowing date and location is an important component of adaptation when breeding wheat lines for target production regions. Virtual genotypes were created including all MLG of Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1, and a full range of values of TT to investigate the possible range of heading times given these known values. APSIM-Wheat-G was then used to simulate heading times across a wide sowing window (from 1 March to 30 September) for the years 1960 to 2009 (SILO Patched Point Dataset) for 1479 stations across the Australian wheatbelt (map from Zheng ). The median heading times were calculated to show the distribution of heading time in the Australian wheatbelt and compared with the probabilities of last frost day (<0 °C) and first heat day (>35 °C) (Zheng ).

Results

Trial observations and differences among MLGs for the VRN1 and Ppd-D1 genes

The VRN and PPD trials in South Perth for 210 lines provided 1359 observed heading times which ranged from 18 August to 19 October (63 d). The validation set of 172 lines at 79 sites (Fig. 2) with sowing times between 30 March to 11 September (165 d) comprised 4475 observed heading times (30 May to 24 November; 47–160 d after sowing). Cultivars Wyalkatchem, Janz, and Gladius were sown in about 50% of experiments.

Lines in the calibration and validation datasets included 14 of 16 (24) possible homozygous MLGs for the four genes in our datasets (Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1; Table 1). Most lines were in three MLGs: 39 lines with Vrn-A1a, Vrn-B1v, Vrn-D1v, and Ppd-D1a; 32 for Vrn-A1v, Vrn-B1a, Vrn-D1v, and Ppd-D1a; and 26 for Vrn-A1v, Vrn-B1a, Vrn-D1a, and Ppd-D1a.

Effects of the VRN1 and Ppd-D1 genes

Pooling data for 2008 and 2009, the heading times in the pre-VRN and extended-PPD treatment (V2P2) were significantly correlated with other treatments (0.62, 0.76, and 0.71 for V1P2, V2P1, and V1P1, respectively), indicating that EPS was a major determinant of heading time. The lowest correlation was found between V2P1 and V1P2 (0.42, P <0.001), indicating a weaker relationship between VRN requirement and PPD sensitivity.

The difference between V2P2 and the natural VRN and extended-PPD trial (V1P2) should estimate the effects of the VRN1 loci (and also unknown VRN genes), while the pre-VRN and natural PPD trial (V2P1) can be used to estimate the effects of Ppd-D1 (and unknown PPD genes). The gene effects of VRN1 and Ppd-D1 were predicted by including the MLGs for the VRN1 and Ppd-D1 genes in the fixed part of REML model.

Combined over the 2 years, the heading times of V2P2 (EPS) ranged from 68 to 96 DAS (i.e. TT of 901–1281 °Cd after sowing; data not shown) with a normal distribution (mean=1026±51 °Cd, Kolmogorov–Smirnov test, P <0.001).

Comparison of the treatments V1P2 and V2P2 (V1P2 – V2P2; Fig. 3 and Supplementary Table S3 at JXB online) indicated the independent effect of VRN when the PPD requirement was satisfied. Ten sets of lines were identified as having one of five MLGs for the VRN1 loci (a-a-v, a-v-v, v-a-a, v-a-v, and v-v-v; Fig. 3 Supplementary Table S3) and being either sensitive or insensitive at the Ppd-D1 locus. For these VRN1 groups, the alleles of Ppd-D1 had a minimal influence on heading times, e.g. the predicted effects of VRN1 genes on heading times were 14.9±13.3 vs 17.8±11.0 °Cd for a-a-v-a and a-a-v-b homozygous genotypes respectively; and 178.1±13.7 vs 209.9±29.5 °Cd for v-v-v-a and v-v-v-b, respectively. The lines with winter alleles at all VRN1 loci exhibited the greatest effects on heading time (e.g. 183.8±12.6 °Cd, averaged for both Ppd-D1 alleles for V1P2 – V2P2). For other MLGs of VRN1, there were strong interactions among alleles of VRN1 genes. For example, the effect of the winter genotype for Vrn-D1 was only 1.3 °Cd when present with spring genotypes of Vrn-A1 and Vrn-B1, while the effect was 25.8 °Cd with spring alleles of Vrn-A1 and winter alleles of Vrn-B1, and was 32.7 °Cd with winter alleles of Vrn-A1 and spring alleles of Vrn-B1 (Fig. 3). Similar interactions were found for Vrn-A1 and Vrn-B1. The average effect of the winter VRN1 genotype compared with the spring genotype in this comparison of VRN treatments was 76.8 °Cd for Vrn-A1, 33.8 °Cd for Vrn-B1 and 43.6 °Cd for Vrn-D1, i.e. compared with the spring genotype, the winter genotype at Vrn-A1 had almost exactly twice the effect as the other genes in terms of delaying heading time. The three VRN1 genes explained almost all variation in VRN effects (V1P2 – V2P2), as the residual VRN effect was only 15.4±17.7 °Cd for MLG that were homozygous for spring alleles at all three VRN1 loci.

Fig. 3.

Estimated heading time difference (°Cd) between V2P2 (‘control’) and other treatments (see text for description of REML method). The MLGs are for Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1 loci (see Table 2). The differences between V1P2 and V2P2 (V1P2 – V2P2) indicate VRN effects, and those between V2P1 and V2P2 (V2P1 – V2P2) indicate PPD effects. Results are shown as ±standard deviation.

The PPD effect contrast (V2P1 – V2P2) for the five VRN1 MLGs discussed above indicated that the winter alleles of VRN1 loci had strong interactions with PPD sensitivity of Ppd-D1. For example, the PPD effect was 120.5±17.8 °Cd for the a-a-v-a genotype, which was substantially less than the 162.2±14.7 °Cd for a-a-v-b. The average effects of the VRN1 winter genotypes on the interaction of the Ppd-D1 sensitivity allele were 41.7 °Cd for Vrn-A1, 99.0 °Cd for Vrn-B1, and 51.8 °Cd for Vrn-D1, i.e. the Vrn-B1 winter alleles resulted in about twice the delay in heading time compared with the other VRN1 alleles. There were apparently other unknown genes influencing the PPD response, as the PPD effect was still 121.5±23.8 °Cd for MLG homozygous for spring alleles of all three VRN1 genes and the insensitive allele of Ppd-D1.

Performance of APSIM-Wheat-O and APSIM-Wheat-M

Using both the calibration and validation datasets together, each genotype was fitted using either the APSIM-Wheat-O or APSIM-Wheat-M. The phenology parameters (R , R , and TT ) of APSIM-Wheat-O were refitted due to poor predictions (RMSE=7.2 d) for several lines using the default values supplied in APSIM version 7.4. After re-fitting, there was good agreement between observed and simulated heading times for APSIM-Wheat-O (RMSE=3.9 d, y=0.99x+1.11, R 2=0.96, P <0.001, N=5843; see Supplementary Fig. S1 at JXB online). Data fits for the APSIM-Wheat-M (with multiplicative instead of minimum function for VRN and PPD effects) were better than for APSIM-Wheat-O (RMSE=3.3 d, y=0.97x+2.73, R 2=0.97, P <0.001, N=5843; Supplementary Fig. S1).

Estimation of parameter values for the gene-based model (APSIM-Wheat-G)

According to the observation of VRN1 and Ppd-D1 gene effects (Fig. 3) and in order to keep the gene-based model simple, we ignored minor epistatic effects among the VRN1 genes. However, we did consider gene interactions between VRN1 and Ppd-D1 effects. Based on the earlier description of results from Fig. 3, the weighted effects (h ) of the different loci on VRN were set to integers 2, 1, and 1 for Vrn-A1, Vrn-B1, and Vrn-D1, respectively, in Equation 8, i.e. winter alleles at the Vrn-A1 locus had double the effect of those at the Vrn-B1 and Vrn-D1 loci. The intercept of Equation 10 was set to zero, as three homoeologous VRN1 genes explained most variation of VRN effects (V1P2 – V2P2; Fig. 3). The weighted effects (h ) due to interaction of the VRN1 and Ppd-D1 genes were set to integers 1, 2, and 1 for Vrn-A1, Vrn-B1, and Vrn-D1, respectively, in Equation 9. As b is set as 0, APSIM-Wheat-G required three gene-specific parameters (k , k , and b ) and one line-specific parameters (TT ) to predict heading time, i.e. the MLG parameters are the same for any lines that have the same combinations of VRN1 and Ppd-D1 alleles. When these parameters were fitted together using an optimization algorithm to minimize the RMSE on the South Perth calibration dataset, R and R were given by the following equations derived from Equations 8–11. The fitted parameters were 0.3 for k , 0.2 for k , and 0.9 for b .

There was good agreement between observed and simulated heading time for the four treatments in the calibration dataset (RMSE=3.2 d, y=0.89x+9.5, P <0.001, N=1359; Fig. 4). There was slightly poorer performance for V2P1 treatment (RMSE=3.7, y=0.71x+25.8, P <0.001, N=340), compared with the V2P2 treatment (RMSE=2.8 d, y=1.07x – 5.92, P <0.001, N=340). We did test APSIM-Wheat-G with the exact ratios (rather than the integer values given above) of estimated gene effects of Vrn-A1, Vrn-B1, and Vrn-D1 (Fig. 3). However, as similar results were obtained, the integer weighted effects were retained for model simplicity.

Fig. 4.

Comparison between observed and simulated heading times for calibration datasets using APSIM-Wheat-G. These trials were conducted over 2 years in South Perth, Western Australia, with four treatments: natural and pre-VRN (V1 and V2, respectively), and natural and extended PPD (P1 and P2, respectively). The total RMSE was 4.1 (y=1.06x+4.9, P <0.001, N=1359; R 2=0.86, dashed line; 1:1, solid line).

The target TTs for the length of the floral initiation to flowering stages (TT ) ranged from 455 to 1025 °Cd across a normal distribution (663±96 °Cd, Kolmogorov–Smirnov test, P <0.001; see Supplementary Fig. S2 at JXB online) with most lines (70%) having a TT between 500 and 800 °Cd.

The parameter values of F , F , and TT of APSIM-Wheat-O, APSIM-Wheat-M, and APSIM-Wheat-G are listed in Supplementary Table S1 for 124 Australian wheat varieties and other commercial spring wheat releases (data on proprietary breeding lines were excluded). Eagles presented alleles of Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1 for many Australian wheat lines, including 56 lines in our datasets, but 10 of them had different alleles in some loci (see Supplementary Table S2 at JXB online). Further research is needed to check these lines. However, Eagles indicated that it is not uncommon that commercial seed from released cultivars could contain mixtures of (homozygous) lines that vary for alleles at one or more loci.

Validation of the gene-based model (APSIM-Wheat-G)

There was a high correlation between observed and simulated heading times for validation datasets using the parameters derived only from the South Perth trial (RMSE=4.3 d, y=0.98x+0.28, P <0.001, R 2=0.96, N=4475; Fig. 5A). Meanwhile, there was also good agreement between APSIM-Wheat-M and APSIM-Wheat-G (RMSE=2.5 d, y=1.00x – 0.93, P <0.001, R 2=0.99, N=5834; Fig. 5B) with only a slight loss of accuracy in APSIM-Wheat-G. The RMSE slightly increased from 3.4 to 4.3 d when line-specific parameters (APSIM-Wheat-M) were replaced by gene-specific parameters (APSIM-Wheat-G) (Fig. 5 and Supplementary Fig. S1). When predicting heading dates of multiple genotypes within single trials, the median RMSE and R 2 were 3.8 d and 0.72, respectively (n=35 trials), and when predicting heading dates of a single genotype across multiple trials, the median RMSE and R 2 were 3.8 d and 0.95, respectively (n=141 genotypes) (see Supplementary Fig. S3 at JXB online). There were no systematic errors related to sowing date or latitude for APSIM-Wheat-G (see Supplementary Fig. S4 at JXB online).

Fig. 5.

Comparison between observed and simulated heading times for validation datasets using APSIM-Wheat- G (A), and simulated heading times between APSIM-Wheat-M and APSIM-Wheat-G for calibration and validation datasets together (B) with a linear fit indicated by dashed line, and 1:1 by a solid line.

Virtual genotypes to predict suitable genotypes for different regions and management

Virtual genotypes were created using all MLGs of Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1, and several combinations of TT to explore the potential range of heading times in this diverse germplasm. The results for the earliest and latest MLGs showed that later sowing times resulted in a smaller range of heading times for sites/years where there were multiple sowing times (Fig. 6). In the more northern latitudes of Gatton (27.55°S, 152.33°E) and Geraldton (28.78°S, 114.61°E), the widest range of heading times occurred for May sowing, while at the southern sites of Katanning (33.68°S, 117.61°E) and Roseworthy (34.53°S, 138.69°E), the widest heading date range was for April planting. Compared with the virtual genotypes, the observed data covered the full range of heading times when sowing was late but did not sample the latest heading times when sowing was early (Fig. 6). For any given trial and sowing date, the majority of lines were mid-maturing, as indicated by the darker colour (higher frequency).

Fig. 6.

Comparison between heading times (day of year, DOY) and sowing times (DOY) for a subset of selected sites and years where large numbers of lines were assessed at multiple sowing times. The lower and upper solid lines show the heading times of the virtual shortest-season lines and longest-season lines, respectively. The shortest-season lines had spring alleles of Vrn-A1, Vrn-B1, and Vrn-D1, an insensitive allele of Ppd-D1, and TT of 455 °Cd. The longest-season lines had winter alleles of Vrn-A1, Vrn-B, and Vrn-D1, a sensitive allele of Ppd-D1, and TT of 1025 °Cd. The coloured strips show the distribution of the number of wheat lines in a trial with a specific observed heading time at a specific sowing date in the validation datasets.

The model allows an interpretation of the potential heading dates for any combination of genotype, location, and sowing date. Early June is the only sowing opportunity that can occur in almost any part of the wheatbelt, and generally suits a ‘mid-season’ maturity type, e.g. with a winter allele at Vrn-A1 (v), spring alleles at Vrn-B1 and Vrn-D1 (a), a sensitive allele at Ppd-D1 (b), and an intermediate TT (800 °Cd), which is representative of the released variety Ellison. As is observed in trials and expected from the temperature climatology of the wheatbelt, Fig. 7 shows that this genotypic combination would have a heading time that gradually became later from north to south and from inland to coast. The earliest heading times would be at the start of September in the west of central Queensland, south-west of Queensland, and north-west Western Australia and in parts of central South Australia (Fig. 7). Later heading times would occur at the end of October in the south of Western Australia and south-east New South Wales. Similar maps of heading time for other virtual genotypes, including lines that are ‘not in the current germplasm’, have been generated as a reference for breeders (data not shown).

Fig. 7.

Median heading times simulated by APSIM-Wheat-G at 1479 locations across the Australian wheatbelt for 1 June sowing (1960–2009). The map shows heading times for four virtual wheat genotypes with a winter allele at Vrn-A1 (v), spring alleles at Vrn-B1 and Vrn-D1 (a), a PPD-sensitive allele at Ppd-D1 (b), and an intermediate thermal time from floral initiation to flowering (400–1000 °Cd). The similar cultivars include Baxter and Ellison.

Wheat growers aim to choose a genotype of suitable maturity to avoid frost and heat stresses in the spring. For an early sowing time (e.g. 1 May) in Moree (29.48°S, 149.84°E; Fig. 8), a grower would need a variety with at least two VRN1 winter alleles (including VRN-A1), a PPD-sensitive allele at Ppd-D1, and a TT of at least 1000 °Cd. By late June/early July, a suitable variety would need a much lower TT (400–600 °Cd, i.e. stronger expression of EPS) and fewer or no VRN1 winter alleles in order to flower early in the ‘safe’ window. For any given sowing date and genotype, the spread of heading times across years is approximately 1 month (the width of boxplots). Therefore, to minimize the risk of coincidence of frost with heading time (i.e. the entire boxplot is to the right of the frost line), a rule of thumb is that the genotype should flower such that its median heading date is approximately 2 weeks after the date of last frost. Maps of heading times have been generated for other stations as references for breeders (data not shown).

Fig. 8.

Impact of sowing time and allele combinations for VRN1 and Ppd-D1 genes and TT on heading times compared with the occurrence of extreme-temperature events (for example, Moree in New South Wales for 1960–2009). The boxplot shows the variation in heading time (x-axis) for different sowing times (y-axis: every half month from 15 April to 15 July) and for MLGs. The MLGs are the alleles of Vrn-A1, Vrn-B1, Vrn-D1, and Ppd-D1. a and v indicate homozygous genotypes for the spring and winter alleles of the VRN1 genes, respectively, while a and b indicate homozygous genotypes for the insensitive and sensitive alleles of Ppd-D1, respectively. The panels show a range of different TT from 400 to 1000 °Cd. Probabilities of last frost days (left solid line) and first heat days (right solid line) are calculated as the percentile of last frost days (<0 °C) and first heat days (>35 °C) from 1960 to 2009. The lower horizontal dashed line indicates the date of 10% risk of last frost day (x-axis) and the upper horizontal dashed line indicates the date of 30% risk of first heat day (x-axis). The low-risk period for frost and heat is highlighted in grey and helps to determine the best sowing window, on the y-axis. See Zheng for methods to calculate last frost days, first heat days, and the low-risk window.

Discussion

Quantifying the effects of the VRN1 and Ppd-D1 genes on heading time

Field trials with VRN and PPD treatments were used to estimate the effects of VRN1 and Ppd-D1 genes in a broad sample of Australian wheat lines. Assuming that the pre-VRN and extended-PPD treatment (V2P2) completely satisfied the VRN and PPD requirements, this treatment should therefore estimate only the effects of EPS. A significant correlation between V2P2 and other treatments indicated that EPS was a major determinant of heading times for winter sowing, as was found by Rousset for autumnal sowings of winter wheat.

The effects of VRN1 genes were estimated as the difference between V1P2 and V2P2 with the long PPD requirement being satisfied. Three homoeologous VRN1 genes (Vrn-A1, Vrn-B1 and Vrn-D1) explained most variation of VRN effects (V1P2 – V2P2; Fig. 3). For lines with three spring genotypes at the VRN1 loci, the VRN effects were only about 15 °Cd, which is about 1–3 d in most of the Australian wheatbelt, indicating that there were no other major genes influencing VRN in these lines. This may be related to the fact that the VRN2 and VRN3 genes are typically fixed to spring alleles in Australian wheat lines (Dion Bennett, personal communication). The results showed that the winter allele of Vrn-A1 had the strongest effect on delaying heading time, with the effects of the winter alleles of Vrn-B1 and Vrn-D1 being about half that of Vrn-A1 (Fig. 3). Vrn-A1 has been suggested elsewhere to have the strongest effects on wheat development (Trevaskis ; Loukoianov ; Allard ). While some researchers indicated stronger winter allele effects of Vrn-B1 than that of Vrn-D1, others found the weaker effects for the winter allele of Vrn-B1 (Loukoianov ; Allard ).

The effect of Ppd-D1 was estimated by the difference between V2P1 and V2P2, in which the VRN requirement was satisfied. There were obviously residual effects of PPD (Fig. 3), as only Ppd-D1 was used to evaluate PPD effects in this study and there were notable effects of the PPD treatment in lines that were insensitive at Ppd-D1. The residual effects may be related to other major PPD genes, e.g. Ppd-A1 and Ppd-B1 (González ), the FT locus (Bonnin ), or other unknown genes. All 210 lines in this study contained an apparently insensitive allele of Ppd-A1 (data not shown). Others have shown that the three homoeologous PPD1 genes (Ppd-A1, Ppd-B1, and Ppd-D1) had different effects on wheat development (González ). Our results showed that alleles of VRN1 genes had a significant influence on the effects of Ppd-D1 (Fig. 3), even when the VRN requirement was satisfied. Li suggested that PPD1 can promote the flowering time of wheat through CO2 (constans) genes and VRN3 genes under long days after the VRN requirement is satisfied. Allard showed that the rate of VRN varied with the change in PPD. In the calibration dataset, Vrn-B1 had the strongest effect in delaying heading time under long days, while Vrn-A1 and Vrn-D1 had smaller delay effects.

Predicting phenology with a gene-based model

The gene-based model was built to predict wheat phenology based on the conventional process-based model APSIM-Wheat with gene parameters (for VRN1 and Ppd-D1) and line-specific parameters of the gene-based model estimated from two sets of VRN and PPD at a single location. For 210 genotypes, APSIM-Wheat-O requires a total of 630 parameters (estimated for each genotype from multiple trials), while APSIM-Wheat-G requires 210 genotype parameters (estimated from one trial), and three gene-specific parameters (fixed values). The accuracy of development models is expected to decrease as line-specific parameters are replaced by gene-specific (‘generic’) parameters (Yin ; White ; Uptmoor ). For example, RMSE increased from 6.0 and 9.0 d to 8.6 and 9.9 d for calibration and validation datasets, respectively, when predicting wheat flowering time (White ). Uptmoor indicated that R 2 was reduced from 0.81 and 0.65 to 0.59 and 0.50, respectively, for low- and high-temperature trials, respectively, in predictions of flowering time of B. oleracea. The performance reduction of gene-based models is caused by undetected effects of minor QTL (Yin ) or unknown genes but also by poor estimation of allele effects of known QTL (Uptmoor ) or genes. Compared with existing gene- and QTL-based models, the performance of the gene-based model in this study is a significant improvement and has been validated over a wide range of environments and almost 4500 observations. The precision was only slightly less than the process-based model APSIM-Wheat (RMSE increased from 3.9 to 4.3; Fig. 5A). White used the same VRN1 and Ppd-D1 genes to predict wheat phenology and validated this at 34 diverse locations with a higher RMSE (9.5 d), but did not account for the interaction of VRN1 and Ppd-D1 loci and did not use a calibration experiment to estimate the effects of the alleles. A source of error in gene-based models is the method of estimation of gene parameter values. In the previous gene-based models, the line-specific parameters of process-based models were first fitted with optimization algorithms (normally by obtaining the minimum RMSE), and then these parameter values were re-estimated as a function of the presence/absence of alleles at QTL/genes and new parameters associated with these QTL/genes (White and Hoogenboom, 1996; Hoogenboom and White, 2003; Messina ; White ). Therefore, the systematic errors of models and uncertainty are accumulated in a two-step optimization. Additionally, as most process-based models are non-linear, there exist multiple parameter sets to optimize the same models. In this study, all gene- and line-specific parameters were fitted together in one step to reduce systematic errors related to optimization.

Only additive effects were considered in previous gene-based models of flowering time (White and Hoogenboom, 1996; Hoogenboom and White, 2003; Messina ; White ). However, many researchers have indicated that variation in the effects of VRN1 and PPD1 genes delay flowering time (Trevaskis ; González ; Loukoianov ; Allard ). The weighted gene effects were introduced into a gene-based model to simulate the different effects of homoeologous genes (Equations 8 and 9).

There is still some variation among lines for the target TT from floral initiation to flowering (TT ), which is the minimum TT requirement when the VRN requirement and PPD sensitivity are satisfied and is related to the trait EPS (Fig. 4 and Supplementary Fig. S2). Several EPS genes have been located on wheat chromosomes 1A, 2B, 3A, 4B, 4D, 5A, 6B, 6D, and 7B (Snape ; Bullrich ), and in barley, many of these genes have been identified and in some cases cloned (Comadran ; Faure ; Zakhrabekova ). In future development of the model, we propose that the line parameter TT can be replaced with gene parameters associated with EPS genes. At present, new lines can be included in our model by genotyping them for the VRN1 and PPD1 alleles and conducting a simple study with one treatment of pre-VRN and extended PPD to determine the TT . Researchers are also beginning to identify alternative VRN1 alleles that differ in effect size, compared with the current spring or winter alleles (Eagles ). Incorporation of these multi-allelic effects in the model would require calibration experiments with all four treatments, as presented here. Using experiments on winter wheats, the model should be able to be extended to include the effects of the ‘winter wheat’ loci VRN2 and VRN3. The simple interaction between the VRN1 and Ppd-D1 genes and allele effects in this study could also be captured by introducing gene regulatory networks and expression levels and underlying biological mechanisms, but this would probably reduce the accuracy of the model (Salazar ; Wenden ). Given the increased interest in the genetic manipulation of spike development time (i.e. influencing grain set per spike; Miralles ), additional experiments and model development would be needed to accurately predict the lengths of development stages that occur prior to heading, e.g. time to end of tillering, time to floral initiation. Others have found that to model these stages with precision requires a greater number of parameters, but that this might still be achieved with only about three parameters per genotype (He ).

Targeting breeding activities for current and future climates

In the Australian wheatbelt, there has been an increase in mean temperature of 0.19 °C per decade over the last 40 years (Murphy and Timbal, 2008) with predictions that warming associated with climate change is expected to be 1–2 °C by 2030. If the same maturity types were grown, this degree of warming would shorten the growing season of wheat by 4–6 weeks (Zheng ). The shorter growing season infers a reduction in the available time to accumulate radiation, CO2, and nutrients for storage of assimilate and resulting grain yield and economic return. It is thus urgent to examine adaptation of wheat genotypes and associated production system to future climate, as breeding cycles take 7–12 years to deliver new lines (Chapman ).

A key consideration for wheat growers is to minimize the risks of frost and heat stresses in the period prior to and soon after heading when the crop is most sensitive to impacts on grain number per spike. In most Australian environments, the crop should flower in a ‘safe’ window, soon after the last frost event, so that the crop also potentially avoids heat events and terminal drought stress. Sowing time depends on a ‘planting’ rain event, so growers aim to choose a genotype of suitable maturity for a given location and sowing time. By understanding and manipulating VRN, PPD, and EPS, breeders can adjust heading date to fit the seasonal temperature and rainfall pattern in current and future environments (Sharma ; Zheng ). There is an expectation that in some regions and management combinations, stronger VRN or PPD combinations may need to be combined with longer TT , especially to counter the effects of increasing seasonal temperatures.

Supplementary data

Supplementary data are available at JXB online.

Supplementary Materials and methods. Detailed information about the validation datasets.

Supplementary Fig. S1. Comparison between observed and simulated heading times (Days after sowing, DAS) with original and modified APSIM-Wheat.

Supplementary Fig. S2. The frequency distribution of target thermal time from floral initiation (FI) to flowering (FL)

Supplementary Fig. S3. The histogram graph of root mean square error and coefficient of determination.

Supplementary Fig. S4. The residuals of (predicted – observed) heading times from APSIM-Wheat-G compared with sowing dates or latitudes.

Supplementary Table S1. The alleles of Vrn-A1, Vrn-B1, Vrn-D1 and Ppd-D1 and the parameter values of APSIM-Wheat phenology model for 125 Australian wheat varieties.

Supplementary Table S2. Comparison of the variety composition of alleles of Vrn-A1, Vrn-B1, Vrn-D1 and Ppd-D1 between our dataset and Eagles .

Supplementary Table S3. Estimated heading time difference between V2P2 (‘control’) and other treatments.