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.

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.

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.

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).

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 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).

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

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: XLij=FiLst, where (s,t)=in(s+S(t-1)),t+S(s-1) and for the additive model XLij=∑sFiLsj 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.

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 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).

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).

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).

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.

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.

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.

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).

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).

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).

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.

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.