Determining Phenological Patterns Associated with the Onset of Senescence in a Wheat MAGIC Mapping Population

The appropriate timing of developmental transitions is critical for adapting many crops to their local climatic conditions. Therefore, understanding the genetic basis of different aspects of phenology could be useful in highlighting mechanisms underpinning adaptation, with implications in breeding for climate change. For bread wheat (Triticum aestivum), the transition from vegetative to reproductive growth, the start and rate of leaf senescence and the relative timing of different stages of flowering and grain filling all contribute to plant performance. In this study we screened under Smart house conditions a large, multi-founder “NIAB elite MAGIC” wheat population, to evaluate the genetic elements that influence the timing of developmental stages in European elite varieties. This panel of recombinant inbred lines was derived from eight parents that are or recently have been grown commercially in the UK and Northern Europe. We undertook a detailed temporal phenotypic analysis under Smart house conditions of the population and its parents, to try to identify known or novel Quantitative Trait Loci associated with variation in the timing of key phenological stages in senescence. This analysis resulted in the detection of QTL interactions with novel traits such the time between “half of ear emergence above flag leaf ligule” and the onset of senescence at the flag leaf as well as traits associated with plant morphology such as stem height. In addition, strong correlations between several traits and the onset of senescence of the flag leaf were identified. This work establishes the value of systematically phenotyping genetically unstructured populations to reveal the genetic architecture underlying morphological variation in commercial wheat.

Introduction

Wheat is a pillar of global food security, providing 20% of protein and calories consumed worldwide and up to 50% in developing countries. It is the main food staple in Central Asia, West Asia and North Africa, which have the world's highest per capita wheat consumption (Valluru et al., 2015). Global wheat production is at risk due to climate change, population growth, changing food preferences and the plant health challenges associated with its widespread cultivation. In order to maintain optimal production and profitability, wheat producers and processors must prepare for and adapt to these challenges. The current emphasis on food security has focused research attention on two avenues to improve wheat yield (Valluru et al., 2015): (1) increasing photosynthetic capacity and efficiency (Reynolds et al., 2009); and (2) increasing partitioning of assimilates to the developing spike and grain.

The timing of key developmental transitions is critical for many crops, but is particularly important in the temperate small grain cereals. For example, the transition from vegetative to reproductive growth can have major effects on biomass accumulation and harvest index that profoundly affect either the locales in which a variety can be profitably grown, or its ultimate use. Thus, crops destined for grain production should transition early, relative to the length of the growing season, to allow ripening, avoid stress, and achieve a high harvest index of grain to total biomass. Forage, biofuel or dual purpose crops could usefully transition later to allow greater total biomass accumulation, but this has to be tempered with the likelihood of deleterious stress/weather events. Flowering time, therefore, has been a key selection target since the beginning of domestication (Izawa, 2007), initially inadvertently but since modern breeding began, very directly.

Understanding the extent and basis of other aspects of phenological variation may be useful in breeding for yield potential and stress adaptation. In wheat, leaves can contribute up to 40% of the nitrogen incorporated by the grains on the fifteenth day after anthesis (Simpson et al., 1983). Therefore, lifespan of the leaves (hence yield) is a trade-off with N remobilization. Delayed leaf senescence (as in the stay-green effect), which maintains active photosynthesis for a longer period, can increase grain yields under certain circumstances (Gregersen et al., 2008). Conversely, accelerated senescence leads to low carbon (C) but high N remobilization, indicating plasticity of C and N remobilization during development, perhaps correlated with senescence. A better understanding the genetic and environmental factors affecting these processes would help optimize C and N remobilization to the actively developing grains under different growth/stress conditions.

While several patterns of senescence have been proposed (Thomas and Howarth, 2000), an ideal senescence phenotype in wheat, and in cereals in general, still needs to be identified (Gregersen et al., 2008), perhaps because of strong and variable environmental effects. In monocarpic crops such as wheat the initiation of senescence typically leads to a massive remobilization of phloem-mobile nutrients from the senescing plant parts to developing sinks, such as seeds or grains (Figure 1; Gregersen et al., 2008; Distelfeld et al., 2014). Pathogen infection also interacts with developmental processes in a complex way and symptoms of senescence often accompany the progression of disease, although senescence can also be delayed in response to pathogen infection (Häffner et al., 2015).

FIGURE 1

Figure 1. Whole leaf senescence flowchart.

A major constraint to progress in breeding for high yield varieties is the access to appropriate and consistent selection environments. The selection environment plays a key role in the efficiency of the selection process. Since environmental variables are almost impossible to control under field conditions, the identification of specific genetic factors associated to crop yield becomes more challenging (Bentley et al., 2013).

Modern controlled environment (CE) growth control and/or recording of variable environmental parameters such as temperature and watering allow for elimination or reduction of uncontrollable influences. Having control and access to an experiment's environmental parameters allows for reproducibility, and decreases the levels of uncertainty as it is easier to reverse-engineer an experiment in order to identify—or at least justify—the causes of a given phenotype, and minimizes the amount of replication per subject due to the low variability of the environment. Therefore, CE phenotyping offers closely defined conditions compared to the relatively homogeneous but less controllable growing conditions in a field plot.

To understand the genetic control of phenology and the onset of senescence in wheat, we screened a core set of the NIAB elite MAGIC wheat population (Mackay et al., 2014) across time. The eight founders of this MAGIC population were selected in partnership with UK wheat breeders to sample trait variation and germplasm important to current UK breeding programmes (Bentley et al., 2013). MAGIC populations combine high levels of genetic diversity, recombination and homozygosity (Mackay et al., 2014) to create a panel of recombinant inbred lines (RILs). A well-designed MAGIC population captures and immortalizes the variation released by intercrossing, thereby providing a stable well-defined population to be shared across sites and used across years.

In this study, plants from the elite MAGIC core set and its parents were scored throughout their life cycle to capture traits, including decimal growth stages, biomass and plant height. We discuss the use of a subset of MAGIC lines in the Smarthouse as a proof of concept that the combination of MAGIC + Smarthouse phenotyping should be repeated on a grander scale.

Materials and Methods

Plant Material

A subset of the NIAB Elite eight-founder MAGIC population described in Mackay et al. (2014) was used for all phenotypic screening. The complete population consists of approximately 1000 recombinant inbred lines (RILs) generated from three cycles of recombination between eight elite United Kingdom wheat varieties (Alchemy, Brompton, Claire, Hereward, Rialto, Robigus, Soissons, Xi-19) followed by five rounds of selfing to derive RILs. Further information about the population, including pedigree, genotype, and existing phenotype data can be found at www.niab.com/MAGIC.

The core set used in this study was selected to represent all funnels of the 210 8-way crosses within the population. Two funnels were not represented due to limited seed availability making a total of 208 RILs in the core set.

Glasshouse Cultivation

Plants were grown between mid-January 2015 and mid April 2015 in The National Plant Phenomics Centre facilities in Aberystwyth, UK. The eight parents of the MAGIC population and four additional elite varieties (Avalon, Santiago, Cadenza, and Zircon) were grown with the 208 RILs (see Table S1) under well watered conditions, with two replicates per genotype. Two seeds were sown in 8 × 8 cm pots of Levington F2 compost. After germination (approximately on the 30/10/2014) the seedlings were thinned to one per pot and transferred to a controlled environment room for vernalization (5°C, 16 h daylength) for 9 weeks. Following vernalization plants were transferred to 15 × 15 × 20 cm pots of M2 compost. Field capacity and dry matter content of the compost was determined. Plants were transferred to the growth chamber where each pot was placed into a cart on a conveyor system. Pots were weighed and watered automatically to 75% gravimetric water content daily. Growth conditions were 14 h daylength using 600W sodium lamps to supplement (350 μM/m²/sec) natural lighting, with the temperature settings of 18°C (day) and 15°C (night). Plant hygiene was monitored by visual inspection throughout the experiment, with an appropriate prophylactic and responsive spraying regime. Once the ears started to ripen, plants were removed from the system and allowed to finish ripening naturally with reduced watering. Reduced watering only occurred after all plants had passed Flag leaf senescence.

Phenotyping

Plants were manually scored for developmental stages according to the Zadoks scale (Zadoks et al., 1974) three times per week, and scored as days after sowing (DAS) when plants reached growth stage 39 (GS39; flag leaf fully emerged), GS55 (ear 50% emerged), GS65 (50% anthesis), and the onset of flag leaf senescence (FLS). At the end of the experiment, plants were harvested and above ground biomass (PW), tiller number (TN), plant height (PH), stem height (SH), top internode length (TIL), first/second/third ear length (FEL, SEL, TEL), first/other ear weight (FEW/OEW), and first flag leaf length (FFLL) were scored. The number of days between GS39 and GS55 (d1), GS55 and GS65 (d2) and GS55 and FLS (d3) were also determined. A multiple linear regression model (MLRM) was fitted to identify predictors of FLS among all the traits used in the analysis. A list of traits, abbreviations and Crop ontology terms (Shrestha et al., 2012) are provided in Table 1.

TABLE 1

Table 1. Plant trait descriptions.

Genotyping

The lines were genotyped using the Illumina Infinium iSelect 80,000 SNP wheat array (“80K array,” http://www.illumina.com/), described in Wang et al. (2014). 20,639 SNP markers were scorable and polymorphic, of which 18,601 were successfully mapped in the MAGIC population (Mackay et al., 2014; Gardner et al., 2016). Linkage map generated by mpMap is reported in Gardner et al. (2016).

Plant Stress

The young plants showed symptoms that included chlorosis and necrosis (Figure S1). The chlorotic symptoms consisted of yellow areas surrounding lesions on the leaf blades. The necrotic symptoms comprised brown spots, lens-shaped lesions, surrounded by yellow borders. Although symptoms were controlled by routine spraying (Priori Xtra, Syngenta), we speculated that it constituted an undiagnosed disease (possibly Septoria) and the degree of infection was scored manually using the seedling infection type (IT) score shown in Table S1. Plant visual stress symptoms were scored first at GS31–39 and for the second time around GS70–GS80, and the average was calculated (SM). These qualitative IT scores were converted to a numerical scale for statistical analysis.

Statistical and Quantitative Trait Locus Analysis

Statistical analyses were performed in the R environment using Core Team (2013). Quantitative Trait Locus (QTL) analysis was performed using the R package HAPPY for multi-parental populations analyses (Mott et al., 2000). The genetic analysis of multi-parental populations requires a haplotype-based approach because single marker association or interval mapping can fail to detect a QTL if the causative alleles are not dispersed among the founders with the same strain distribution pattern as the linked markers (Mott et al., 2000).

QTL Mapping

HAPPY's analysis is essentially two stage; ancestral haplotype reconstruction using dynamic programming, followed by QTL testing by linear regression:

• Assume that at a QTL, a pair of chromosomes originating from the progenitor strains, labeled s, t contribute an unknown amount T_st to the phenotype. In the special case where the contribution from each chromosome is additive at the locus then T_st = T_s + T_t, say

• A test for a QTL is equivalent to testing for differences between the T's.

• A dynamic-programming algorithm is used to compute the probability F_iLst that a given individual i has the ancestral alleles s, t at locus labeled L, conditional upon all the genotype data for the individual. Then the expected phenotype is

$\begin{array}{l}y={\displaystyle \sum _{st}}{T}_{st}{F}_{iLst},\end{array}$

and the T's are estimated by a linear regression of the observed phenotypes on these expected values across all individuals, followed by an analysis of variance to test whether the progenitor estimates differ significantly.

• The method's power depends on the ability to distinguish ancestral haplotypes across the interval.

• All inference is based on regression of the phenotypes on the probabilities of descent from the founder loci, F_nst.

The models are presented here in the linear model framework (i.e., least-squares estimation, with ANOVA F-tests).

For an additive QTL, the parameters are the strain effect sizes; for a full interaction model there is a parameter for every possible strain combination. Then the one-QTL model is E(y) = X_Lt_L.

There are S(S − 1)/2 + S parameters (where S is the number of strains) to be estimated in a full model allowing for interactions between the alleles within the locus, and S − 1 parameters in an additive model. For the full model, the i, j'th element of the design matrix X is related to the strain probabilities thus:

$\begin{array}{l}{X}_{Lij}=F_{iLst},\end{array}$

where

$\begin{array}{l}\left(s,t\right)=in(s+S\left(t-1\right)),t+S(s-1)\end{array}$

and for the additive model

$\begin{array}{l}{X}_{Lij}={\displaystyle \sum _s}{F}_{iLsj}\end{array}$

We used an additive model, where the contribution of each allele at the locus are assumed to act additively.

Furthermore, when mapping QTLs in structured populations the evidence for the existence of a QTL has to be considered in the context of other QTLs, which might explain some of the same component of variation. Population structure can produce long-range correlations between genotypes and hence ghost QTL, although the LD analysis suggests that the MAGIC population is relatively immune to this phenomenon. Although the MAGIC population is relatively unstructured, and therefore can be analyzed one locus at a time, in order to ensure the evidence for a given QTL was not confounded with that for others, statistical significance was assessed based on permuting the phenotypes (1000 times) between individuals, repeating the model fit, and finding the top-scoring marker interval. The empirical distribution of the max −logP values was then used to assess statistical significance. This technique is useful for non-normally distributed phenotypes and for estimating region-wide significance levels. We used −logP = 4 as a threshold in the multiple QTL modeling to test for association and FDR = 0.05 to identify significantly differential markers. The dashed line in QTL plots corresponds to an FDR rate of 0.05 and is calculated using the qvalue package (Storey and Tibshirani, 2003). The p-value corresponding to a q-value of 0.05 is determined by interpolation. When there are no q-values less than 0.05, the dashed line is omitted.

In addition, a multiple linear regression model (MLRM) was fitted to identify predictors of FLS among all the traits used in the analysis. We selected FLS because the trait is used as indicator of crop yield and biomass accumulation (Gan, 2014).

Principal Component Analyses (PCA) over normalized trait data was carried out to identify patterns between traits and genotypes. Biplots were used to show information on samples in a graphical manner (Kempton, 1984). PCA of marker data was carried out separately to test for population structure.

QTL Validation

We used the R package mpMap to confirm QTL mapping results (Huang and George, 2011) and to analyse the effect of including marker covariates. QTL analysis is performed using interval mapping, then selected marker covariates are included in the linear model in a forward selection process.

Data Processing

Data were pre-processed using standard methods. Data corresponding to one replicate of the MAGIC line MEL 086-1 (note, all MAGIC lines are named with prefix “MEL”) and another from one replicate of the elite line Avalon were removed due to seed infection. A small number of outliers (data points with suspicious values) were checked and, where possible, corrected. Missing values (2.19%) were imputed using multivariate imputation by chained equation (MICE) (van Buuren and Groothuis-Oudshoorn, 2011). Briefly, MICE operates under the assumption that given the variables used in the imputation procedure, the missing data are Missing At Random (MAR), which means that the probability that a value is missing depends only on observed values and not on unobserved values (Schafer and Graham, 2002). MICE creates a number of datasets by imputing missing values. That is, one missing value in original dataset is replaced by m plausible imputed values. We set m = 5 as the number of imputations. These values take imputation uncertainty into consideration. Statistics of interest are estimated from each dataset and then combined into a final one and replicates were averaged (Zhang, 2016).

Results

In this experiment, RILs, MAGIC parents (illustrated in Figure 2) and 4 other elite genotypes were grown to maturity over a single time span and within a single glasshouse chamber with controlled watering and supplementary lighting and heating. Phenotype data were curated, and missing values (accounting for 2.19% of data points) imputed. Figure S2 shows a comparison between original data (red dots) and imputed data (blue dots), which suggested a high similarity between the two distributions as indicated by the overlapping dots.

FIGURE 2

Figure 2. MAGIC parents at senescence.

Analysis of Traits

The distributions of the traits are shown in Figure S3. Most traits showed similar distributions with the exception of the discrete trait SM, which was skewed (see Table 1 for trait description). SM is a discrete trait and the skewedness of the plot reflects that most RILs' scores were in the 0–2 range. Figure 3 shows the frequency distributions for GS55, FLS, SH and d3. Pair-wise correlation analysis between all traits (Figure S3) identified strong correlations between FLS and GS39 (0.79), GS55 (0.73) and GS65 (0.69) and d3 (0.70); between FEL, SEL, and TEL (>0.86) between TIL and SH (0.76), PW and OEW (0.87), and between TEW and PW (0.87).

FIGURE 3

Figure 3. Frequency distributions of GS55, FLS, SH, and d3 from all RILs.

To determine if there was variation in duration between key developmental stages that was not simply a result of variation in overall developmental progression, we examined the time in days taken to progress from GS39 to GS55 (d1), from GS55 to GS65 (d2) and from GS55 to FLS (d3). Our MLRM identified d3 as a strong predictor (P < 0.05) of FLS, indicating that the time lapse between GS55 and FLS is a good candidate to predict FLS (Figure 4). In addition to d3, our MLRM identified other important predictors of FLS (P < 0.05). For example, the size of the flag leaf on the primary shoot, FFLL, was significantly (P < 0.05) associated with the timing of senescence (FLS). To demonstrate this result, dot size in Figure 4 was used to represent an additional feature in the plots. In the case of Figure 4, FFLL (represented by dot size) was longer in RILs that senesced earlier. Among the MAGIC founders, Brompton, Hereward and Rialto senesced after Xi-19 and the two elite controls Zircon and Cadenza. The latter also had the shortest duration between GS39 and GS55. Previously, Mackay et al. (2011) reported Cadenza as the most environmentally sensitive variety detected in 8 years in Recommended List trials showing a linear increase in yield with increasing summer rainfall. This supports the observations from the trait analysis that progression through different developmental processes e.g., flowering vs. senescence, is controlled independently. Figure 5 shows d3 in relation to FLS (the contract of d2 to FLS is shown in Figure S5). Plants which senesced earlier took less time between GS55 and GS65 (d2) and between GS55 and FLS (d3).

FIGURE 4

Figure 4. Variation between d1 and FLS. Dot size represents FFLL/100 (divided by 100 to be able to show the dots in the plot). Blue dots correspond to RILs. Other colored dots correspond to MAGIC parents.

FIGURE 5

Figure 5. Variation between d3 and FLS, where d3 = days between GS55 and FLS. Dot size represents disease score/10. Blues dots correspond to RILs and other colored dots to MAGIC parents.

To further evaluate the relationship between these traits, PCA was conducted over trait data scaled to have unit variance. Results of the analysis are shown on the biplot in Figure 6. The plot shows that PC1 and PC2 account for 48% of the total variance of the traits. Also, four clearly defined trait groups (anti clockwise) can be seen in the plot, the first one containing FLS, GS39, GS55 and GS65, d1, d3, and TN, the second group contained SEL, FEL and TEL, the third group contained OEW, PW, TIL, SH, FEW, and FFL and the fourth group contained HI and SM. Since (1) the smaller the angle between the trait vectors, the higher the correlation (2) trait values are smaller toward the middle of the plot and higher toward the edge, we can deduce that SM is negatively correlated to group one which is confirmed by the results from the correlation shown in Figure S4. The same argument could also be used between group one and group three, which is also confirmed by Figure S4.

FIGURE 6

Figure 6. (Left panel) PCA score biplot for 18 different traits. (Right panel) Biplot corresponding to 8 MAGIC parent, 4 elite lines and 208 RILs. MAGIC parents are labeled only to facilitate interpretation.

Looking at the relationship between traits and MAGIC parents, Brompton and Hereward are positively correlated, tending to be slow to FLS (late to senesce) as indicated by their proximity to the FLS vector. In contrast, Soissons and Xi-19 show an opposite effect, rapidly reaching FLS as indicated by their location at the opposite side of FLS. Through this analysis we can also see that, in general, most RILs and parents have a similar overall phenome, as represented by their location close to the center of the plot. We can also see there are a number of divergent phenotypes, such as the one at the bottom of the plot (Figure 6, right panel) which corresponds to MAGIC line MEL 091-1a or the one at the top which correspond to MEL 089-1a. Looking closer, MEL 089-1a is proximal to PW and SH while MEL 091-1a is further away, which indicates these two lines contrast strongly for these particular traits. A picture of both lines taken on 27/04/2015 was added to the plot to facilitate interpretation. The plants show clearly contrasting differences in height and biomass.

Plant Stress Analysis

A low level of chlorotic and necrotic lesions was observed on the leaves early during the growth period. Symptoms (SM) were scored independently by two people at two time points and scores were averaged. When comparing symptom scores against the onset of FLS, we found that the more severely affected plants started senescence earlier than those plants that were mildly affected (Figure 5). Our MLRM also identified SM as a predictor of FLS (P < 0.05).

Figure 6 confirms the negative correlation between FLS and SM as indicated by its opposite location from FLS. This correlation is consistent with that of Mycosphaerella graminocola, where infection induces senescence by manipulating signaling pathways in plants (Mengiste, 2012). However, it should be noted that the precise identity of the putative pathogen could not be confirmed.

Quantitative Trait Loci

We evaluated whether trait variation could be ascribed to underlying genetic variation. The lines were genotyped using the Illumina Infinium iSelect 80,000 SNP wheat array (“80K array,” http://www.illumina.com/), described in Wang et al. (2014). 20,639 SNP markers were scorable and polymorphic, of which 18,601 were successfully mapped in the MAGIC population (Mackay et al., 2014; Gardner et al., 2016); linkage map for this population was produced using mpMap and reported in Gardner et al. (2016). First, we checked for signs of population structure. To do this, the marker based relationship matrix (A) was calculated using the R package rrBLUp (Endelman, 2011), then a PCA analysis by eigenvalue decomposition of A was calculated. Results are shown in Figure S7. This shows that the first PC accounts for less than 4% of the total spectrum. This confirms the expected absence of population structure. Genome mosaics corresponding to MEL 15-2, MEL 091-1a, and MEL 209-1 are shown in Figure S8. These decompose the lines' genomes into mosaics of founder haplotypes. The lines appear to be a random mix of the founders, which indicates an absence of gross population structure.

After confirming the absence of population structure, we used HAPPY (Mott et al., 2000) to test for association between each phenotype and the predicted founder haplotypes at each locus in the genome. We used the P-value threshold < 10⁻⁴ to call QTLs [−logP = 4, corresponding to a false discovery rate (FDR) = 0.05]. This analysis identified loci associated with two phenological traits, GS39 and GS55, and a number of traits such as SEL, SH, TIL, TEW, and HI. For GS39, three significant QTLs were found on chromosome 5A, at 201.36, 212.52, and 224.64 cM, corresponding to the markers BS00009369_51, BS00021942_51 and wsnp_Ex_c5978_10478584, respectively (Table 2, Figure 7). For GS55, three QTLs were found on chromosome 5A, at 201.36, 216.05 and 227.66 cM, which correspond to BS00009369_51, wsnp_Ex_rep_c66689_65011117, and Excalibur_c7729_144, respectively (Figure 8). In both cases, these three close peaks are likely to represent a single QTL. To confirm this hypothesis, we performed composite interval analysis using 5A as covariate and identified a single clear and strong marker on 5A (Figures S11A–E).

TABLE 2

Table 2. QTLs mapped for different traits.

FIGURE 7

Figure 7. QTLs markers for the trait GS39. Dashed lines indicate threshold (FDR < 0.05, logP4).

FIGURE 8

Figure 8. QTLs markers for the trait GS55. Dashed lines indicate threshold (FDR < 0.05, logP4).

For SEL, four QTLs were identified on chromosome 5B at 60.65, 66.38, 90.8, 92.31 cM, corresponding to wsnp_Ex_c6548_11355524, BS00001101_51, wsnp_Ku_c2185_4218722, and RAC875_c19099_434, respectively (Figure 9).

FIGURE 9

Figure 9. QTLs markers for the trait SEL. Dashed lines indicate threshold (FDR < 0.05, logP4).

For SH, one QTL was identified on chromosome 4D at 32.24 cM which corresponded to the marker RAC875_c1673_193 (Figure 10). A QTL for TIL was also identified in the same location (Figure 11). QTLs for OEW and PW were also identified on 4D at 26.97 cM corresponding to RAC875_c6922_291 (Figure S10A) and 4D at 32.24 cM corresponding to RAC875_c6922_291 (Figure S10B), respectively. We identified one QTL for TEW on chromosome 4D at 26.97 cM, corresponding to RAC875_c6922_291 (Figure S10D). All of these QTLs co-located with the semi-dwarfing gene Rht-D1 (Rht2) (Ellis et al., 2002). The Rht (Reduced height) genes Rht-B1 (Rht1) or Rht-D1 (Rht2) are present in many high-yielding, semi-dwarf varieties, where they offer simple genetic control of high harvest index and resistance to lodging (Flintham et al., 1997).

FIGURE 10

Figure 10. QTLs markers for the trait SH. Dashed lines indicate threshold (FDR < 0.05, logP4).

FIGURE 11

Figure 11. QTLs markers for the trait TIL. Dashed lines indicate threshold (FDR < 0.05, logP4).

For HI, 1 weak peak QTL on chromosome 2D at 55.4 cM was identified which corresponded to RAC875_c6922_291 (Figure S10E). Figure S9 shows a contrast between FLS and HI across all the MAGIC and elite lines. Cadenza and Soissons have some of the highest HIs and the shortest time to senesce. Table 2 also shows that Soissons has the highest contribution (0.6) to that particular marker.

Discussion

This study screened a core set of lines derived from the NIAB wheat MAGIC population under Smarthouse conditions as strategy to understand the physical and genetic relationship between different phenological traits. This strategy resulted in the detection of QTL interactions with novel traits suggesting that the methodology should be taken further in the future.

Pair-wise correlation between all traits identified high positive (≥0.69) correlations between GS39, GS55, GS65 and senescence at the flag leaf (FLS), and between TEW and PW. There was also a negative correlation between FLS and length of the first flag leaf (FFLL), indicating that the shorter the flag leaf the more delayed was the start of senescence. Short flag leaves provide less nutrient assimilation therefore the plant has to compensate by either living longer, or producing a large number of tillers. Consistent with this idea, PC analysis also indicated a correlation between these traits. Suggestions that delayed leaf senescence leads to increased yield have been thrown into doubt (Borrill et al., 2015) but it may contribute under certain conditions. It will be interesting to see whether the correlation between leaf length and senescence is maintained under other environmental conditions. Genetically unstructured populations such as the MAGIC collection will be ideal to test whether experimental manipulation of the size of the sink (grain mass) can further modulate flag leaf senescence.

In order to see if there was any clustering of individuals according to phenotype, a projection plot of the MAGIC lines onto the first two PCs was generated (Figure S6). The map shows lines grouped around 3 clusters of traits where Group 1 contained Cadenza, Zircon, Xi-19, Claire, Alchemy, Soissons, and Robigus; Group 2 contained Brompton and Hereward and Group 3 contained Avalon, showed similar trait profiles. Group 1 was early to senescence and group 2 later. Group 3 contained smaller plants as indicated by their opposite location to the FEL, SEL, and TEL vectors.

FLS was also negatively correlated with disease resistance (SM), indicating that highly susceptible plants were more likely to trigger senescence early. Support for this perspective may come from the observation that Xi-19 had the earliest FLS of all the parents and controls (Figure S9) and the joint highest disease score. In the field, Xi-19 flowers considerably later than the earliest flowering parent, Soissons, which carries the Ppd-D1a allele for early flowering (Scarth et al., 1985).

After confirming absence of population structure with PCA of the kinship matrix, QTL mapping identified a very strong marker (−logP > 8.00) on chromosome 5A associated with GS39 and GS55, which we believe is likely to correspond to the vernalization gene VRN-A1. This gene plays an important role in the vernalization process in diploid (Triticum monococcum) and polyploid wheat (Triticum aestivum) (Loukoianov et al., 2005; Kiss et al., 2014). However, using HAPPY, no significant QTL was found on chromosome 2D, the location of the Ppd-D1 locus, for GS39 or GS55. In the field, the presence of the Ppd-D1 allele in Soissons results in this line flowering 7–14 days earlier than the other MAGIC founder lines. These contrasting results between field and CER for Ppd-D1 and Vrn-A1 associated QTL suggest that the plants in this experiment might have experienced reduced vernalization as a result of a lack of cold treatment. However, this idea was discarded because inspection of CER records did not show any temperature discrepancies during vernalization. Another possibility is that plants displayed disease-like symptoms at an early stage. We also noticed that Cadenza, one of the 4 elites and a genotype that does not need vernalization, was one of the first to senesce, had the shortest duration between GS39 and GS55 and was more disease susceptible than the similarly early-flowering Soissons (Figures 4, 5). The fact that Cadenza has no vernalization requirement, might suggest that indeed plants we not fully vernalized or encountered a de-vernalizing effect. Whatever the cause, the results of the experiment appear to have been strongly affected by a vernalization issue. This may explain some of the “anomalous” behavior of Cadenza, Xi-19 and Zircon, all of which do not require vernalization.

Further insight is provided by QTL validation analyses using mpMap. With no covariates included in the QTL model, the mpMap interval mapping approach produces very similar results to HAPPY. However, many more QTLs (−logP > 10) are detected using a model with 10 covariates in mpMap, as can be seen in Figure S11. For GS39 and GS55, it can be seen that although the 5A QTL is still the highest peak, the Ppd-D1 marker is significant and detected as the third (GS39) or 2nd (GS55) highest QTL. For FLS, 5A is also the most prominent QTL, but there is no evidence for a QTL around the Ppd-D1 locus. This supports the observations from the trait analysis that progression through different developmental processes e.g., flowering vs. senescence, is controlled independently. Furthermore, QTL detected using mpMap for the length of the interval (d3) from flowering to senescence (Figure S11C) show a distinct pattern from both GS55 and FLS, although some loci are in common (e.g., 4D). As well as phenological traits, markers associated to morphological traits were also identified. For example, a strong marker on chromosome 4D was associated with shoot height (SH) and shoot number (TIL), as well as ear weight (OEW) and above ground biomass (PW). In all these cases, the QTL interval includes the Rht-D1 (Rht2) locus. The Rht-D1b allele at this locus causes a semi-dwarfing phenotype in wheat, is strongly correlated with a reduction in height and several other morphological traits and is segregating in the MAGIC population (dwarfing alleles are present in all parent lines except Robigus and Soissons). Interestingly, this locus also shows up in the highest QTL interval for the d3 developmental interval in the mpMap covariate analysis (Figure S11D), suggesting that there may be independent effect on timing of progression through development. The Rht-D1b allele has a premature stop codon resulting in reduced sensitivity to gibberellic acid, which has been associated with reduced plant height and earlier heading date (Wilhelm et al., 2013). Our analysis indicates there are differential effects on the duration of other developmental processes not directly related to height or flowering per se.

In addition to FLS, GS55, and GS39, the multiple covariate analysis also identified a strong (−logP > 18) peak on chromosome 7B (Figure S11E). Chromosome 7 has been previously associated to Septoria leaf blotch in an analysis of wheat-barley disomic addition lines. The highest level of resistance to infection by S. tritici was found in the H. vulgare chromosome addition line 7 followed by 4 and 6 (Rubiales et al., 2001).

Another interesting result is related to HI, for which a weak peak QTL was identified on the 2D chromosome, the location of Ppd-D1. In the field, the presence of the Ppd-D1 allele in Soissons results in flowering 7–14 days earlier than the other MAGIC founder lines. In our analysis, Soissons have some of the highest HIs and the shortest time to senesce but it also has the highest contribution (0.6) to that particular maker.

This study provides the first systematic phenological characterization of a wheat MAGIC population under controlled environment conditions. With careful developmental staging and end of life measurements of the MAGIC core set we were able to identify previously detected QTL loci on chromosomes 5A and 4D associated with the onset of senescence at the flag leaf. This powerful multi-founder population captures much of the genetic variation present in elite cultivars and a more detailed knowledge of fine-scale developmental and physiological patterns can be exploited for fine-tuning wheat's response to the environment. We have shown that the combination of MAGIC + Smarthouse can help extend the current understanding of developmental plasticity in elite wheat varieties with potential application for responding to the adaptation challenges facing agriculture in a changing climate.

Author Contributions

Conceived and designed the study: AB and AC; analyzed the data: AC; assisted with QTL analysis: RM, KG, and IM; provided genetic data: KG and IM; provided scoring data: FC; wrote the paper: AC; provided comments and corrected the manuscript: All authors.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer FM and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

Acknowledgments

Access to the National Plant Phenomics Centre was provided by a National Capability for Crop Phenotyping grant (BBSRC ref number BB/J004464/1). The creation of the NIAB MAGIC population was supported by BB/E007201/1. We are grateful to the team of National Plant Phenomics Centre for carrying out the experiments, particularly Julie Pruvost. Thanks are extended to Prof. Luis A. J. Mur (Aberystwyth University, UK) and Dr. Flavio M. Santana (Embrapa Wheat) for critical discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01540

References

Bentley, A., Jensen, E., Mackay, I., Hönicka, H., Fladung, M., Hori, K., et al. (2013). Genomics and Breeding for Climate-Resilient Crops, Vol 2. Berlin; Heidelberg: Springer-Verlag.

Borrill, P., Fahy, B., Smith, A. M., and Uauy, C. (2015). Wheat grain filling is limited by grain filling capacity rather than the duration of flag leaf photosynthesis: a case study using NAM RNAi plants. PLoS ONE 10:e0134947. doi: 10.1371/journal.pone.0134947

PubMed Abstract | CrossRef Full Text | Google Scholar

Core Team (2013). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available online at: http://www.R-project.org/

Distelfeld, A., Avni, R., and Fischer, A. M. (2014). Senescence, nutrient remobilization, and yield in wheat and barley. J. Exp. Bot. 65, 3783–3798. doi: 10.1093/jxb/ert477

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, M., Spielmeyer, W., Gale, K., Rebetzke, G., and Richards, R. (2002). “Perfect” markers for the Rht-B1b and Rht-D1b dwarfing genes in wheat. Theor. Appl. Genet. 105, 1038–1042. doi: 10.1007/s00122-002-1048-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Endelman, J. B. (2011). Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4, 250–255. doi: 10.3835/plantgenome2011.08.0024

CrossRef Full Text | Google Scholar

Flintham, J. E., Börner, A., Worland, A. J., and Gale, M. D. (1997). Optimizing wheat grain yield: effects of Rht (gibberellin-insensitive) dwarfing genes. J. Agric. Sci. 128, 11–25. doi: 10.1017/S0021859696003942

CrossRef Full Text | Google Scholar

Gan, S.-S. (2014). Leaf senescence as an important target for improving crop production. Adv. Crop Sci. Tech. 2:e116. doi: 10.4172/2329-8863.1000e116

CrossRef Full Text | Google Scholar

Gardner, K. A., Wittern, L. M., and Mackay, I. J. (2016). A highly recombined, high-density, eight-founder wheat MAGIC map reveals extensive segregation distortion and genomic locations of introgression segments. Plant Biotechnol. J. 14, 1406–1417. doi: 10.1111/pbi.12504

PubMed Abstract | CrossRef Full Text | Google Scholar

Gregersen, P. L., Holm, P. B., and Krupinska, K. (2008). Leaf senescence and nutrient remobilisation in barley and wheat. Plant Biol. 10, 37–49. doi: 10.1111/j.1438-8677.2008.00114.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Häffner, E., Konietzki, S., and Diederichsen, E. (2015). Keeping control: the role of senescence and development in plant pathogenesis and defense. Plants 4:449. doi: 10.3390/plants4030449

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, B. E., and George, A. W. (2011). R/mpMap: a computational platform for the genetic analysis of multiparent recombinant inbred lines. Bioinformatics 27, 727–729. doi: 10.1093/bioinformatics/btq719

PubMed Abstract | CrossRef Full Text | Google Scholar

Izawa, T. (2007). Adaptation of flowering-time by natural and artificial selection in Arabidopsis and rice. J. Exp. Bot. 58, 3091–3097. doi: 10.1093/jxb/erm159

PubMed Abstract | CrossRef Full Text | Google Scholar

Kempton, R. A. (1984). The use of biplots in interpreting variety by environment interactions. J. Agric. Sci. 103, 123–125. doi: 10.1017/s0021859600043392

CrossRef Full Text

Kiss, T., Balla, K., Veisz, O., Láng, L., Bedõ, Z., Griffiths, S., et al. (2014). Allele frequencies in the VRN-A1, VRN-B1 and VRN-D1 vernalisation response and PPD-B1 and PPD-D1 photoperiod sensitivity genes, and their effects on heading in a diverse set of wheat cultivars (Triticum aestivum L.). Mol. Breed. 34, 297–310. doi: 10.1007/s11032-014-0034-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, C., Jaiswal, P., Hebbard, C., Avraham, S., Buckler, E. S., Casstevens, T., et al. (2008). Gramene: a growing plant comparative genomics resource. Nucleic Acids Res. 36, D947–D953. doi: 10.1093/nar/gkm968

PubMed Abstract | CrossRef Full Text | Google Scholar

Loukoianov, A., Yan, L., Blechl, A., Sanchez, A., and Dubcovsky, J. (2005). Regulation of VRN-1 vernalisation genes in normal and transgenic polyploid wheat. Plant Physiol. 138, 2364–2373. doi: 10.1104/pp.105.064287

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, I., Horwell, A., Garner, J., White, J., McKee, J., and Philpott, H. (2011). Reanalyses of the historical series of UK variety trials to quantify the contributions of genetic and environmental factors to trends and variability in yield over time. Theor. Appl. Genet. 122, 225–238. doi: 10.1007/s00122-010-1438-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, I. J., Bansept-Basler, P., Barber, T., Bentley, A. R., Cockram, J., Gosman, N., et al. (2014). An eight-parent multiparent advanced generation inter-cross population for winter-sown wheat: creation, properties, and validation. G3 4, 1603–1610. doi: 10.1534/g3.114.012963

PubMed Abstract | CrossRef Full Text | Google Scholar

Mengiste, T. (2012). Plant immunity to necrotrophs. Ann. Rev. Phytopathol. 50, 267–294. doi: 10.1146/annurev-phyto-081211-172955

PubMed Abstract | CrossRef Full Text | Google Scholar

Mott, R., Talbot, C. J., Turri, M. G., Collins, A. C., and Flint, J. (2000). A method for fine mapping quantitative trait loci in outbred animal stocks. Proc. Natl. Acad. Sci. U.S.A. 97, 12649–12654. doi: 10.1073/pnas.230304397

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynolds, M., Foulkes, M. J., Slafer, G. A., Berry, P., Parry, M. A. J., Snape, J. W., et al. (2009). Raising yield potential in wheat. J. Exp. Bot. 60, 1899–1918. doi: 10.1093/jxb/erp016

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubiales, D., Moral, A., and Martín, A. (2001). Chromosome location of resistance to septoria leaf blotch and common bunt in wheat-barley addition lines. Euphytica 122, 369–372. doi: 10.1023/A:1012952819255

CrossRef Full Text | Google Scholar

Scarth, R., Kirby, E. J. M., and Law, C. N. (1985). Effects of the photoperiod genes Ppd1 and Ppd2 on growth and development of the shoot apex in wheat. Ann. Bot. 55, 351–359.

Google Scholar

Schafer, J. L., and Graham, J. W. (2002). Missing data: our view of the state of the art. Psychol. Methods 7, 147–177. doi: 10.1037/1082-989X.7.2.147

PubMed Abstract | CrossRef Full Text | Google Scholar

Shrestha, R., Matteis, L., Skofic, M., Portugal, A., McLaren, G., Hyman, G., et al. (2012). Bridging the phenotypic and genetic data useful for integrated breeding through a data annotation using the Crop Ontology developed by the crop communities of practice. Front. Physiol. 3:326. doi: 10.3389/fphys.2012.00326

PubMed Abstract | CrossRef Full Text | Google Scholar

Simpson, R. J., Lambers, H., and Dalling, M. J. (1983). Nitrogen redistribution during grain growth in wheat (Triticum aestivum L.): IV. Development of a quantitative model of the translocation of nitrogen to the grain. Plant Physiol. 71, 7–14. doi: 10.1104/pp.71.1.7

PubMed Abstract | CrossRef Full Text | Google Scholar

Storey, J. D., and Tibshirani, R. (2003). Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. U.S.A. 100, 9440–9445. doi: 10.1073/pnas.1530509100

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, H., and Howarth, C. J. (2000). Five ways to stay green. J. Exp. Bot. 51, 329–337. doi: 10.1093/jexbot/51.suppl_1.329

PubMed Abstract | CrossRef Full Text | Google Scholar

Valluru, R., Reynolds, M. P., and Lafarge, T. (2015). Food security through translational biology between wheat and rice. Food Energy Security 4, 203–218. doi: 10.1002/fes3.71

CrossRef Full Text | Google Scholar

van Buuren, S., and Groothuis-Oudshoorn, K. (2011). mice: multivariate imputation by chained equations in R. J. Stat. Softw. 45:67.

Google Scholar

Wang, S., Wong, D., Forrest, K., Allen, A., Chao, S., Huang, B. E., et al. (2014). Characterization of polyploid wheat genomic diversity using a high-density 90 000 single nucleotide polymorphism array. Plant Biotechnol. J. 12, 787–796. doi: 10.1111/pbi.12183

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, E. P., Boulton, M. I., Al-Kaff, N., Balfourier, F., Bordes, J., Greenland, A. J., et al. (2013). Rht-1 and Ppd-D1 associations with height, GA sensitivity, and days to heading in a worldwide bread wheat collection. Theor. Appl. Genet. 126, 2233–2243. doi: 10.1007/s00122-013-2130-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zadoks, J. C., Chang, T. T., and Konzak, C. F. (1974). A decimal code for the growth stages of cereals. Weed Res. 14, 415–421. doi: 10.1111/j.1365-3180.1974.tb01084.x

CrossRef Full Text | Google Scholar

Zhang, Z. (2016). Multiple imputation with multivariate imputation by chained equation (MICE) package. Ann. Trans. Med. 4:30.

PubMed Abstract | Google Scholar