Temperature response of wheat affects final height and the timing of stem elongation under field conditions

Temperature response of stem elongation in wheat grown under field conditions indicated that temperature response is highly heritable and linked to the flowering pathway.


Introduction
Temperature is a major abiotic factor affecting plant growth and development. As a consequence of global warming, wheat production could decrease by 6% for each degree Celsius of global temperature increase (Asseng et al., 2015). While heat stress during critical stages can drastically reduce yield (Gibson and Paulsen, 1999;Farooq et al., 2011), warm temperatures can decrease yield by accelerating development and thereby shortening critical periods for yield formation (Fischer, 1985;Slafer and Rawson, 1994). Despite the clear effect of temperature on growth and phenology, little is known about the genotype-specific response pattern to varying temperature conditions during crop development and its genetic control. We therefore aimed to quantify the genotype-specific temperature responsiveness of European winter wheat during the stem elongation (SE) phase.
SE is a critical phase for yield formation in wheat. It occurs between the phenological stages of terminal spikelet initiation and anthesis (Slafer et al., 2015). The start of SE coincides with the transition from vegetative to reproductive development, when the apex meristem differentiates from producing leaf primordia to producing spikelet primordia (Trevaskis et al., 2007;Kamran et al., 2014). During SE, florets are initiated at the spikelets until booting (Kirby, 1988;Slafer et al., 2015). An increased duration of SE increases the number of fertile florets due to longer spike growth and higher dry matter partitioning to the spike (González et al., 2003). This in turn increases the number of grains per spike and therefore yield (Fischer, 1985). Modifying the timing of the critical phenological stages (transition to early reproductive phase and flowering), and thus SE duration, has been proposed as a way to increase wheat yield (Slafer et al., 1996;Miralles and Slafer, 2007;Whitechurch et al., 2007) or at least to mitigate adverse climate change effects on yield, for example by enhancing earliness to escape heat during flowering (Chapman et al., 2012;Hernandez-Ochoa et al., 2019). The recent warming trend causes a faster advancement in phenology. For example, flowering time occurred earlier in Germany throughout the past decade, which is attributable to both increased temperature and selection for early flowering (Rezaei et al., 2018).
Final height is also an important yield determinant. During the 'Green Revolution' wheat yields increased by the introduction of reduced height (Rht) genes. The resulting dwarf and semi-dwarf varieties benefit from improved resource allocation from the stem to the spike and reduced lodging, allowing more intensive nitrogen application (Hedden, 2003). Gibberellin (GA)-insensitive Rht genes (Rht-A1, Rht-B1, and Rht-D1) were shown to limit cell wall extensibility which decreases growth rates (Keyes et al., 1989) without affecting development (Youssefian et al., 1992). Moreover, the allele Rht-B1c (Wu et al., 2011) and the GA-sensitive Rht12 dwarfing gene (Chen et al., 2013) delay heading.
The main abiotic factors affecting the timing of floral initiation and flowering are temperature and photoperiod, with temperature affecting both vernalization and general rate of development (Slafer et al., 2015). These developmental transitions are controlled by major genes involved in the flowering pathway, namely vernalization (Vrn), photoperiod (Ppd), and earliness per se (Eps) genes (Slafer et al., 2015). The Ppd and Vrn genes define photoperiod and vernalization requirements which jointly enable the transition to generative development and define time to flowering. On the other hand, Eps genes fine-tune the timing of floral transition and flowering, after vernalization and photoperiod requirements are fulfilled (Zikhali and Griffiths, 2015). While vernalization and photoperiod response are well known, the role of temperature per se remains less clear. Temperature affects all developmental phases, and warmer ambient temperatures generally accelerate growth and development in crops Rawson, 1994, 1995a,c;Atkinson and Porter, 1996;Fischer, 2011;Slafer et al., 2015). However, it is unclear if temperature response governs growth rate and development independently. If so, the question remains as to whether there is enough genetic variability in temperature response to be used in a breeding context (Parent and Tardieu, 2012).
Genotypic variation for growth response to temperature was reported for wheat leaf elongation rate (Nagelmüller et al., 2016), as well as for canopy cover growth (Grieder et al., 2015). Kiss et al. (2017) reported significant genotype×temperature interactions in the timing of SE as well as temperaturedependent differences in the expression of Vrn and Ppd genes under controlled conditions. Under field conditions, the response of stem elongation to temperature has not yet been investigated at high temporal resolution.
In recent years, new high-throughput phenotyping technologies have enabled the monitoring of plant height with high accuracy and frequency in the field (Bendig et al., 2013;Friedli et al., 2016;Holman et al., 2016;Aasen and Bareth, 2018;Hund et al., 2019). We have previously demonstrated that the ETH field phenotyping platform (FIP; Kirchgessner et al., 2016) can be used to accurately track the development of canopy height in a large set of wheat genotypes using terrestrial laser scanning (Kronenberg et al., 2017). Considerable genotypic variation was detected for the start and end of SE which correlated positively with final canopy height (Kronenberg et al., 2017).
While many temperature-independent factors affecting plant height are known, the influences of temperature-dependent elongation and timing of the elongation phase is less clear. We hypothesize that apart from temperature-independent factors, there is a genotype-specific response to ambient temperature which affects growth as well as the timing of developmental stages. To address this, we aimed to dissect the process towards final height into the following components: (i) temperatureindependent elongation; (ii) temperature-dependent elongation; and (iii) the duration of the elongation phase determined by the start and end of the process. To achieve this, we present a method to assess and measure these three processes under field conditions by means of high-frequency, high-throughput phenotyping of canopy height development. The resulting data were combined with genetic markers to identify quantitative trait loci (QTLs) controlling the aforementioned processes.

Experimental set-up, phenotyping procedures, and extracted traits
Field experiments were conducted in the FIP at the ETH research station in Lindau-Eschikon, Switzerland (47.449°N, 8.682°E, 520 m a.s.l.; soil type: eutric cambisol). We used a set of ~330 winter wheat genotypes (335-352 depending on the experiment) comprising current European elite cultivars (GABI wheat; Kollers et al., 2013), supplemented with 30 Swiss varieties. These were monitored over three growing seasons in 2015, 2016, and 2017. Briefly, the field experiments were conducted in an augmented design with two replications per genotype using micro plots with a size of 1.4×1.1 m. In the 2017 growing season, the experiment was repeated again, with minor changes in genotypic composition. This resulted in 328 genotypes present across all three experiments. Details about the experimental set-up for the growing seasons 2015 and 2016 are described in Kronenberg et al. (2017).
Canopy height was measured twice weekly from the beginning of shooting (BBCH 31) using a light detection and ranging (LIDAR) scanner (FARO R Focus3D S 120; Faro Technologies Inc., Lake Mary, USA) mounted on the FIP . Height measurements were concluded when no further height increase was observed in any of the genotypes. Canopy height data were extracted from the LIDAR data as described in Kronenberg et al. (2017).
The start, end, and duration of SE as well as final canopy height (FH) were extracted from the height data following Kronenberg et al. (2017): FH was defined as the point after which no increase in height was observed in several consecutive measurements. Normalized canopy height was calculated as the percentage of FH on each day of measurement for every plot and then linearly interpolated between measurement points. Growing degree-days after sowing until 15% FH (GDD 15 ) and 95% FH (GDD 95 ) were used as proxy traits for the start and end of SE, respectively. SE duration was recorded in thermal time (GDD SE ) as well as in calendar days (time SE ), as the difference between GDD 95 and GDD 15 (Kronenberg et al., 2017). Heading date was recorded as GDD (heading GDD ) when 50% of the spikes were fully emerged from the flag leaf sheath (BBCH 59; Lancashire et al., 1991). Heading data for 2015 could not be evaluated due to insufficient data availability. Therefore, a third year of heading data was gathered in 2018, when the experiment was repeated again as described in Anderegg et al. (2020).
In order to investigate short-term growth response to temperature, average daily stem elongation rates (SERs) were calculated for each plot as the difference (∆) in canopy height (CH) between consecultive time points (t): Extracting growth response to temperature Temperature response was modelled by regressing average daily SER against average temperature of the respective interval for each plot within the respective year following where T is the ambient temperature, a is the coefficient of the linear regression (i.e. growth response to ambient temperature; slp SER~T ), and ε denotes the residual error. b Tcrit is the model intercept at the temperature at which the correlation between intercept (int SER~T ) and slope is zero (see below). Per definition, the intercept of a linear model would be calculated at T=0 °C-far outside the range of observed temperatures. In the observed data, the intercept at T=0 °C correlated strongly negatively with the slope (Fig. 1A) and thus did not add much additional information concerning the performance of the evaluated genotypes. Likewise, an intercept at 20 °C, at the upper range of the observed data, correlated strongly positively with the slope (Fig. 1A). Grieder et al. (2015) performed a similar analysis for the canopy cover development during winter and found a similar, strongly negative correlation between temperature response (slope) and growth at 0 °C (intercept). We sequentially calculated the intercept at 0.01 °C increments between 1 °C and 22 °C for each plot within a year. Subsequently, we calculated the Pearson correlation coefficient between the common slope of these models and each of the different intercepts ( Fig. 1A). Based on this sequence, we empirically determined the critical temperature value (T crit ) at which the correlation between slope and intercept was zero (Fig. 1A). Hence, T crit is defined as the temperature at which intercept and slope are independent. Due to this independence, the value of the intercept at T crit can be interpreted as the intrinsic growth component independent of temperature response herein referred to as 'vigour'. Following this, two genotypes can show the same vigour but differ markedly in temperature response (Fig. 1B), have the same temperature response but differ in vigour (Fig. 1C), or differ for both temperature response and vigour (Fig. 1D).

Statistical analysis
The data were processed stepwise as follows: (i) correction for design factors and spatial trends; (ii) application of linear models to each plot to determine growth response to temperature; and (iii) prediction of adjusted means and calculation of the heritability for all traits across years. The spatial correction of canopy height per measurement time point was done using the R package SpATS (Rodríguez-Álvarez et al., 2018) following: where Y is the phenotypic value for the a plot in the jth row in the kth column planted with the ith genotype, f(r,c) is a smoothed bivariate surface defined over rows and ranges of a virtual grid, G i is the effect of the ith genotype (i=1, …, n; n=335-352 depending on the year), R j is the effect of the jth row in the virtual grid, and C k is the effect of the kth range. With the number of genotypes, the number of ranges/rows varied across years (17/21 in 2015, 18/20 in 2016, and 18/21 in 2017, respectively). The replications were arranged diagonally in this grid with a gap of five rows and ranges between them. Thus, for example in 2017, replication 1 was at rows 1-21 and ranges 1-18; replication 2 was at rows 24-41 and ranges 27-47 of the virtual grid. The function f(r,c) describing the bivariate surface can be decomposed in a nested-type ANOVA structure as described by Rodríguez-Álvarez et al. (2018). The number of spline points was set to two-thirds of the total number of rows and ranges in the virtual grid, respectively.
From this model, the predicted genotypic best linear unbiased estimates (BLUEs) plus residual error were kept as spatially corrected plot values. Thus, these new plot values were corrected for spatial effect as well as for the random row and range effects, and used for the subsequent dynamic model. For a visualization of the applied SpATS correction, see Supplementary Fig. S1 at JXB online.
Genotypic BLUEs across the three seasons were calculated for all traits using the R package asreml-4 (Butler, 2018) following: where y pij is the spatialy corrected plot value of the respective trait (FH, GDD 15 , GDD 95 , GDD SE , time SE , int SER~T , or slp SER~T ), μ is the overall mean, G i the fixed effect of the genotypes common in all three years (i=1,…,328), Y j is the fixed effect of the year (j=2015,…,2017), GY ij is the random genotype-by-year interaction, and ε pij is the residual error. In order to estimate best linear unbiased predictors (BLUPs) and heritability (H 2 ) across years, G i in Equation 3 was set as a random term and heritability was calculated following (Falconer and Mackay, 1996) using: where H 2 is the broad sense heritability, σ 2 G is the genotypic variance, σ 2 GY is the genotype×year interaction variance, and σ 2 ε is the residual variance. For heading data, only one replicate per year was available. Plotcorrected values were extracted using SpATS (Rodríguez-Álvarez et al., 2018), and heritability across 3 years was calculated by omitting the GY term in Equation 4 and dividing the residual variance by three, based on the three available year-site replications.
Genotypic BLUEs across 3 years were used for subsequent correlation analysis and genome-wide association study (GWAS). All statistical analyses were performed in the R environment (R Core Team, 2015).In order to investigate the relationship between FH, temperature response, and vigour, and to test for confounding Rht or Ppd effects on temperature response, FH was modelled using the linear model where y i is FH of the ith genotype, µ is the model intercept, β 1-6 are the main effect estimates of x=slp SER~T , int SER~T , GDD SE , Rht-B1, Rht-D1, or Ppd-D1, respectively. β 1,2 -β 5,6 are all two-way interaction effects (n=15) and ε i is the residual error. Genotypic data for Rht-B1, Rht-D1, and Ppd-D1 alleles were available for 301 genotypes obtained from Kollers et al. (2014). There, genotyping of the Rht-1 alleles was performed using PCR markers (Ellis et al., 2002), while Ppd-D1 alleles were genotyped by the presence or absence of a 2 kb insertion using specific primers (Beales et al., 2007;Kollers et al., 2014) Association study GWAS was performed on the different traits to compare the phenotypic correlations with the underlying genetic architecture of the traits. As a positive control, FH data from Germany and France reported by Zanke et al. (2014b) were also compared and analysed.
Genotyping data were made previously by the GABI wheat consortium represented by the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK; Zanke et al., 2014a) using the 90K illumina SNP-chip (Cavanagh et al., 2013;Wang et al., 2014). Monomorphic single nucleotide polymorphisms (SNPs) were discarded. The remaining markers were mapped to the IWGSC reference genome (International Wheat Genome Sequencing Consortium, 2018) by BLASTN search using an E-value threshold <1e -30 . The genome position with the lowest E-value was assigned as the respective marker location. Markers that could not be unequivocally positioned were dropped. After filtering SNPs with a minor allele frequency and missing genotype rate <0.05, a total of 13 450 SNP markers and 315 genotypes remained in the set. The reference genome position of Rht, Ppd, Vrn, and putative Eps genes was determined with BLASTN search as described above using published GenBank sequences (Supplementary Table S1).
To mitigate against multiple testing, relatedness, and population structure, three different methods were used to calculate marker-trait associations (MTAs) between phenotypic BLUPs and SNP markers. (i) We used a mixed linear model (MLM) including principal components among marker alleles as fixed effects and kinship as random effect to account for population structure (Zhang et al., 2010). This approach was chosen to stringently prevent type I errors. The MLM GWAS was performed using the R Package GAPIT (v.2, Tang et al., 2016). Kinship was estimated according to VanRaden (2008). (ii) In a generalized linear model (GLM) framework implemented in PLINK (Purcell et al., 2007), association analysis was performed using SNP haplotype blocks consisting of adjacent SNP triplets. Using haplotype blocks takes the surrounding region of a given SNP into account, thus increasing the power to detect rare variants (Purcell et al., 2007). (iii) Finally, the FarmCPU method  was used, which is also implemented in GAPIT. FarmCPU tests individual markers with multiple associated markers as covariates in a fixed effect model. Associated markers are iteratively used in a random effect model to estimate kinship. Confounding between testing markers and kinship is thus removed while controlling type I error, leading to increased power . For all methods, a Bonferroni correction was applied to the pointwise significance threshold of α=0.05, to avoid false positives. Hence, only markers above -log 10 (P-value) >5.43 were considered significant.
Linkage disequilibrium (LD) among markers was estimated using the squared correlation coefficient (r 2 ) calculated with the R package SNPrelate (Zheng et al., 2012). A threshold of r 2 =0.2 (Gaut and Long, 2003) was applied to calculate the chromosome-specific distance threshold of LD decay. Putative candidate genes were identified by searching the IWGSC annotation of the reference genome (International Wheat Genome Sequencing Consortium, 2018) for genes associated with growth and development within the LD distance threshold around the respective MTA.

Phenotypic results
We measured the canopy height of 710-756 plots per year, containing 335-352 wheat genotypes, for three consecutive years. In each season, measurements were made between 17 and 22 times during SE. Plot-based growth rates within single years extracted from these data indicate a clear relationship between growth and temperature for the period of SE, as depicted in Fig. 2. Towards the end of the measurement period in June, there was a larger deviation, which was also reflected in the quality of plot-based linear model fits of SER versus temperature (see Equation 2), summarized in Supplementary Fig. S2. For the 2015 and especially the 2016 experiment, R 2 values were low and except for the 2017 experiment, the parameter estimates were not statistically significant ( Supplementary Fig. S2a). Inspection of the best and worst model fits, however, shows that failure of fitting the model for single plots was levelled out by the replications within genotypes ( Supplementary Fig. S2b). The weak model fits therefore did not affect the genotype ranking of adjusted means across replications. ANOVA revealed significant (P<0.001) genotypic effects for both slp SER~T and int SER~T across 3 years. Both traits showed high heritabilities across years (H 2 =0.81 for slp SER~T and H 2 =0.77 for int SER~T ; Table 1). Using the BLUEs of slp SER~T , int SER~T , and temperature sum for SE (GDD SE ), FH could be predicted with high accuracy across different years (0.85≤R 2 ≤0.89) by training a linear model on the BLUEs of one year and predicting it on the BLUEs of another independent year. In order to account for possible confounding Rht and Ppd effects, the allelic status of Rht-B1, Rht-D1, and Ppd-D1 was included as contrasts in the model (Table 2). Of the 301 genotypes with available data, 7% and 58% carried the dwarfing alleles Rht-B1b and Rht-D1b, respectively, and 13% carried the photoperiod-insensitive allele Ppd-D1a. Training the model on the 3 year BLUEs resulted in a prediction accuracy of single years between R 2 =0.94 and R 2 =0.95 (Fig. 3). Type II ANOVA revealed significant effects for slp SER~T , int SER~T , GDD SE , Rht-B1, and Rht-D1. A significant (P<0.05) interaction effect was found between Rht-D1 and Ppd-D1. Furthermore, weak interactions (P<0.1) were found for Rht-B1:Ppd-D1, int SER~T :Rht-B1, and int SER~T :Rht-D1 (Table 2). High heritabilities across 3 years (0.54≤H 2 ≤0.98; Table 1) were also found for the other traits: FH, start of SE, end of SE, SE duration, and heading. All traits showed moderate genotype×year interaction effects which were smaller (except for SE duration) than the genotypic effects across years (Table 1).

Phenology, temperature response, and final height were positively correlated
To evaluate the relationships between the traits measured, Pearson correlation coefficients were calculated for each trait pair. If not indicated otherwise, the reported correlations were highly significant (P<0.001) Positive correlations were found among GDD 15 , GDD 95 , and FH (0.36≤r≤0.64, Fig. 4), indicating that taller genotypes were generally later in their development towards FH. Temperature response (slp SER~T ) and vigour (int SER~T ) also showed a strong, positive relationship to FH (r=0.85 and r=0.65, respectively). However, only temperature response correlated with GDD 15 and GDD 95 (r=0.63 and r=59, respectively), whereas vigour did not (r≤0.26, Supplementary Fig. S3).

Linkage disequilibrium and population structure
Prior to MTA analysis, we evaluated population structure and LD. Principal component analysis (PCA) of the marker genotypes revealed no distinct substructure in the investigated population. The biplot of the first two principal components showed no apparent clusters, with the first component explaining 8% and the second component explaining 3.3% of the variation in the population (Supplementary Fig. S4). This is consistent with prior work using the same population (Kollers et al., 2013;Yates et al., 2019). On average across all chromosomes, LD decayed below an r 2 of 0.2 at a distance of 9 Mb. There was, however, considerable variation in this threshold among the single chromosomes (Supplementary Table S2).

Association study
Genome-wide association results differed markedly depending on the applied model. Using an MLM with kinship matrix and PCA as covariates resulted in no significant MTA for any trait ( Supplementary Fig. S5). In contrast, the GLM using the haplotype method on temperature response yielded 2958 significant MTAs for α<0.05 and 1852 MTAs for α<0.001, respectively ( Supplementary Fig. S6). However, investigation of the respective QQ-plots showed large P-value inflation in the haplotype method whereas the P-values were slightly deflated when using the MLM approach (Supplementary Figs S5, S6).  Table 2). The model was applied on the 3 year BLUEs of all genotypes with available allelic data of the respective genes (n=301).  In contrast, with FamCPU, the QQ-plots (Fig. 5) showed no P-value inflation, except for some markers. This pattern is expected, if population structure is appropriately controlled. Therefore, FarmCPU was chosen to be the most appropriate method for the given data, despite identifying fewer significant MTAs. As a positive control, we compared our FH data and associated markers with data from Zanke et al. (2014b) who used the same population and SNP chip in field experiments in France and Germany. FH correlated strongly between the two studies (r=0.95), which is in accordance with the high heritability of the trait. In this study, we found 10 significant MTAs for FH (Table 3; Fig. 5). Zanke et al. (2014b) reported 280 significant MTAs for FH across several environments. Of these, only marker RAC875_rep_c105718_585 on chromosome 4D overlapped with the MTAs found in this study. However, by considering flanking markers, we found that of the remaining nine significant MTAs for FH, four were in LD with MTAs found by Zanke et al. (2014b;Supplementary Table S3). The significant MTA found for FH in this study are near known genes controlling FH. For example, Tdurum_contig64772_417 is 4 Mb upstream of Rht-B1 and RAC875_rep_c105718_585 is 7 Mb downstream of Rht-D1 on their respective group 4 chromosomes.

Temperature response loci are independent of vigour loci
For slp SER~T , we detected one significant (LOD=5.75) MTA on chromosome 1B (wsnp_Ex_c1597_3045682) and two almost significant (LOD=5.38 LOD=5.01) MTAs on chromosomes 4B (CAP7_c10839_300) and 5D (IAAV7104), respectively (Fig. 5). All associated markers for slp SER~T yielded small but significant allelic effects ranging from -0.061 mm °C -1 d -1 to -0.051 mm °C -1 d -1 (Table 3). The GWAS for int SER~T yielded four significant MTAs on chromosomes 2B, 4B, 4D, and 5D, respectively (Table 3; Fig. 5). Start and end of SE yielded four MTAs each, and heading yielded eight MTAs (Table 3; Fig. 5). Comparing the GWAS results for temperature response, vigour, FH, GDD 15 , GDD 95 , and heading revealed no common QTLs between slp SER~T and any other trait. Only one marker (Excalibur_ c74858_243) was significantly associated with both GDD 15 and GDD 95 , as well as heading. The lack of overlap of MTAs between temperature response, vigour, and timing of critical stages indicates that they are genetically independent. However, there is a genetic connection between vigour and FH on the one hand and between the start and end of SE and heading on the other.
To identify potential causative genes underlying the QTLs, we searched the reference genome annotation around the respective QTL intervals. For temperature response, we found an increased presence of genes or gene homologues involved in the flowering pathway, namely EARLY FLOWERING 3, FRIGIDA, and CONSTANS (Table 4). Around the QTLs associated with vigour, the annotation showed genes associated with growth (i.e. GRAS, CLAVATA, BSU1, and ARGONAUTE) as well as developmental progress (i.e. Tesmin/ TSO1-like CXC domain, BEL1, and AGAMOUS) (Table 5). Importantly, we found GAI-like protein 1 6 Mb upstream of marker Kukri_rep_c68594_530, which we identified as Rht-D1 by blasting the Rht-D1 sequence (GeneBank ID AJ242531.1) against the annotated reference genome. Genes putatively underlying the QTLs for heading, GDD 15 , and GDD 95 are listed in Supplementary Tables S4-S6. As expected, genes or gene homologues associated with the flowering pathway were found in the vicinity of the MTAs for heading. The common QTLs for heading, GDD 15 , and GDD 95 on chromosome 5B (Excalibur_c74858_243) were found to be 6.6 Mb upstream of FLOWERING LOCUS T (Supplementary Tables S4-S6). Other flowering-associated genes found near the heading QTLs were CONSTANS, FRIGIDA, and a FLOWERING LOCUS C-associated gene (Supplementary Table S4).
Moreover, a number of putative response regulators as well as genes putatively involved in light control of development (i.e. FAR1-RELATED SEQUENCE; Lin and Wang, 2004) were found near the heading QTLs (Supplementary Table S4). The remaining QTLs for GDD 95 were near genes associated with developmental progress and flowering, such as AGAMOUS, MEI2-like 1, HAPPLESS 2, and BEL1 (Supplementary Table  S5). Genes near the remaining QTLs for GDD 15 were associated with developmental progress (i.e. FLOWERING LOCUS T, BEL1, TERMINAL EAR1-like, and FAR1-RELATED SEQUENCE) as well as growth (i.e. CLAVATA and DELLA; Supplementary Table S6).

Vigour, temperature response, and the timing of SE affect final height
The phenotypic correlations show a strong connection between temperature response, vigour, and FH as well as weaker connections between GDD 15 , GDD 95 , and FH. In order to examine this interdependency on a genetic level, we used a linear model to predict FH with the SNP alleles of the QTLs for slp SER~T , int SER~T , GDD 15 , and GDD 95 as predictors. The model was able to predict FH with an accuracy R 2 =0.5; however, clusters in the data showed clear effects of Rht-D1 and Ppd-D1 alleles (Fig. 6A). Adding Rht-B1, Rht-D1, and Ppd-D1 alleles as predictors increased the prediction accuracy to R 2 =0.71 (Fig. 6B). There were significant contributions by QTLs of all three traits; however, their effects were small compared with the obvious effects of Rht-B1, Rht-D1, and Ppd-D1 (Table 6). Including all two-way interaction effects among the QTLs, Rht-1 and Ppd-1 increased the prediction accuracy to R 2 =0.87 (Fig. 6C), indicating a fine-tuning effect of temperature response, vigour, and timing of SE on FH.

Discussion
In this study, we present a method to measure temperature response during stem elongation of wheat using highthroughput phenotyping of canopy height in the field. We found a genotype-specific response of wheat to change in ambient temperature which was correlated with the timing of the developmental key stages. We decomposed this growth dynamic into a genotype-specific vigour component and temperature response component using regression models. We further related these parameters to plant height and the timing of developmental key stages. Linear regression models were used to describe wheat growth response to temperature for leaf elongation (Nagelmüller et al., 2016), canopy cover (Grieder et al., 2015), as well as SER (Slafer and Rawson, 1995a). Others proposed the use of a more complex, Arrhenius type of peak function to account for decreasing growth rates at supra-optimal temperatures (Parent and Tardieu, 2012). However such models are mainly applicable when the temperatures experienced by the crop exceed the temperature optimum. Wheat has its temperature optimum at ~27 °C (Parent and Tardieu, 2012). Temperatures in the measured growth intervals during SE did not exceed 25 °C and, given the temporal resolution of the data, a simple linear model is justified (Parent et al., 2019).
The results of the correlation analysis show a clear connection between FH and temperature response (slp SER~T ) as well as between FH and vigour (int SER~T ). This is consistent with our hypothesis that FH can be described as a function of temperature-independent growth processes and as a function of temperature response during SE. Importantly, among all components, the temperature response was a significant driver of FH and also had a strong influence on the timing. Temperature response delayed the beginning of SE, leading to a later start and end of the whole phase. This finding might appear counter-intuitive: given the assumption that plants develop faster under higher ambient temperatures, a more responsive genotype should develop faster compared with a less responsive one. Slafer and Rawson (1995b) reported an accelerated development towards floral transition with increasing temperatures up to 19 °C, whereas higher temperatures slowed development. In that respect, a more responsive genotype would experience a stronger delay of floral transition under warm temperatures.
In terms of their correlation to FH, the effects of the timing of start and end of SE are less distinct. FH was more a function of faster growth than of a longer duration of growth, especially since genotypes with a strong temperature response had a shorter duration of SE. However, the timing of the start and end of SE was linked with temperature response. Based on this result and the correlations, it would appear that temperature response influences FH directly as well as indirectly by mediating the start and end of SE. Surprisingly, we found no correlation between heading and FH despite the positive correlation of both traits with GDD 95 . A correlation between heading date and FH would therefore be expected. Previous studies reported pleiotropic effects between plant height and heading time (Griffiths et al., 2010;Mo et al., 2018).
The question of whether these trait correlations are due to pleiotropic effects will substantially impact the breeding strategy (Chen and Lübberstedt, 2010). If the relationship between phenology, FH, and temperature response were to be to a large degree pleiotropic, these traits could not be independently selected. Alternative explanations are linkage and population structure. The GABI wheat panel is made of wheat varieties from different regions of Europe. As the examined traits are major drivers of adaptation to the different regions of Europe, we anticipate a very strong selection for both temperature response and timing of critical stages. Even if there is no apparent population structure at neutral markers, there may be a strong population structure at selected loci with a strong effect on local adaptation. Our phenotypic results showed a significant interaction effect between Ppd-D1 and Rht-D1 on FH, indicating either a co-selection or a pleiotropic effect of Ppd-D1. Pleiotropic effects between height and flowering time are known for maize and rice. For example, the DWARF8 gene of maize encoding a DELLA protein is associated with height and flowering time (Lawit et al., 2010) and strongly associated with climate adaptation (Camus-Kulandaivelu et al.,  2006). The rice GHD7 locus has a strong effect on number of days to heading, number of grains per panicle, plant height, and stem growth (Xue et al., 2008). In wheat, the dwarf gene Rht-12 was shown to have a delaying effect on heading (Worland et al., 1994;Chen et al., 2013) as well as an additive interaction effect with Ppd-D1 on plant height (Chen et al., 2018). Furthermore, it was shown that the tall Rht-D1a and the photoperiod-sensitive Ppd-D1b allele positively affect leaf area and spike length throughout SE (Guo et al., 2018). To further examine the relationship among the different traits, we consider the following GWAS analysis using stringent correction of population structure. The GWAS results indicate an independent genetic control of FH, temperature response, and the timing of critical stages, whereas vigour and FH as well as heading time, and start and end of SE appear to be partly linked. Yet, FH could be predicted with surprising accuracy using the QTLs for temperature response, vigour, and start and end of SE, which reflects the correlations found in the phenotypic data.
Previous studies investigating the control of developmental key stages in wheat with respect to temperature generally adopted the concept that after fulfilment of photoperiod and vernalization, Eps genes act as fine-tuning factors independent of environmental stimuli (Kamran et al., 2014;Zikhali and Griffiths, 2015). Increasing temperature, apart from vernalization, is thought to generally quicken growth and development independent of the cultivar (Slafer and Rawson, 1995b;Porter and Gawith, 1999;Slafer et al., 2015). A genotype-specific temperature effect on the duration of different phases was not considered (Takahashi and Yasuda, 1971;Slafer and Rawson 1995c). It was, however, reported that photoperiod effects vary depending on temperature (Slafer and Rawson, 1995c). Under long days, Hemming et al. (2012) reported faster development and fewer fertile florets under high compared with low temperatures. Temperature-dependent effects were also found for different Eps QTLs (Slafer and Rawson, 1995c;Gororo et al., 2001). It has previously been suggested that Eps effects could be associated with interaction effects between genotype and temperature fluctuations (Slafer and Rawson, 1995c;van Beem et al., 2005).
The mechanisms of ambient temperature sensing and of its effects on growth and development are not yet well understood (Sanchez-Bermejo and Balasubramanian, 2016). However, important findings regarding ambient temperature effects on flowering time as well as on hypocotyl elongation have come from Arabidopsis thaliana (Wigge, 2013). With respect to these two traits, Sanchez-Bermejo and Balasubramanian (2016) reported distinct genotypic differences in temperature sensitivity. According to their results, the flowering pathway genes FRIGIDA (FRI), FLOWERING LOCUS C (FLC), and FLOWERING LOCUS T (FT) are major candidate genes for ambient temperature-mediated differences in flowering time (Sanchez-Bermejo and Balasubramanian, 2016). In the present study, we found FRI homologues near two of the three QTLs for temperature response. FRI and FLC act as the main vernalization genes in A. thaliana (Johanson et al., 2000;Amasino and Michaels, 2010). In wheat, these genes are not yet well described. However, FLC orthologues were found to act as flowering repressors regulated by vernalization in monocots (Sharma et al., 2017).
The most promising candidate gene for temperature response found near the QTLs on chromosome 1B is EARLY FLOWERING 3 (ELF3). In Arabidopsis, ELF3 was found to be a core part of the circadian clock involved in ambient temperature response (Thines and Harmon, 2010). In barley, ELF3 was shown to be involved in the control of temperature-dependent expression of flowering time genes (Ejaz and von Korff, 2017). A mutant ELF3 accelerated floral development under high ambient temperatures while maintaining the number of seeds (Ejaz and von Korff, 2017). Furthermore, ELF3 has been reported as a candidate gene for Eps1 in Triticum monococcum (Alvarez et al., 2016) as well as in wheat (Zikhali et al., 2016). A recent study in wheat showed an interaction between Eps-D1 and ambient temperature which corresponded to different expression of ELF3 (Ochagavía et al., 2019). In this study, we directly measured growth response to temperature during SE and found a significant MTA near ELF3 on chromosome 1B. Following Ochagavía et al. (2019), this indicates that growth response to temperature is connected to Eps-B1 which is a homologue to Eps-D1. Furthermore, the temperature×Eps-D1 interaction effects on heading reported by Ochagavía et al. (2019) are in agreement with the correlations found among growth response to temperature, GDD 15 , GDD 95 , and heading in the present study.
One important aspect we could not address in this study is the interaction of genotype-specific temperature response with vernalization and photoperiod (Slafer and Rawson, 1995c;Gol et al., 2017;Kiss et al., 2017). Due to the climate conditions in Switzerland, we expect fulfilment of vernalization requirement in all genotypes. However, due to the broad geographic origin of the investigated genotypes, the relationship between temperature response and the timing of SE might be confounded by different photoperiod requirements. Nevertheless, the correlations between earliness and temperature response are in agreement with Ochagavia et al. (2019). It also remains unclear whether and to what extent temperature response varies across different developmental phases and how temperature response relates to other environmental stimuli such as vapour pressure deficit or radiation. Nevertheless, the results of this study present valuable information towards a better understanding of temperature response in wheat and may be of great importance for breeding. Temperature response could provide a breeding avenue for local adaptation as well as the control of plant height.
With the recent advancements in unmanned aerial vehicle (UAV)-based phenotyping techniques, the growth of canopy cover and canopy height can be measured using image segmentation and structure from motion approaches (Bareth et al., 2016;Aasen and Bareth, 2018;Roth et al., 2018). Thus, temperature response can be investigated during the development of the vegetative canopy cover (Grieder et al., 2015) and during the generative height development as demonstrated here. It can also be assessed in indoor platforms (e.g. Parent and Tardieu, 2012) and the field using a leaf length tracker (Nagelmüller et al., 2016) measuring short-term responses of leaf growth to diurnal changes in temperature. Combining this information may greatly improve our understanding about the genetic variation in growth response to temperature.
Together, the results of this study indicate that temperature response may be exploitable as a breeding trait to adjust phenology towards specific environments, through either phenotypic or marker-assisted selection. Furthermore, a better understanding of temperature response may enhance the capability of crop models to predict crop performance under future climate change scenarios.

Conclusion and outlook
Modern phenotyping platforms hold great promise to map the genetic factors driving the response of developmental processes to environmental stimuli. To the best of our knowledge, this is the first experiment dissecting the SE process into its underlying components: temperature-dependent elongation, temperature-independent vigour, and duration of elongation. The independent loci detected for these traits suggest that it is possible to select them independently. The detected loci may be used to fine-tune height and the beginning and end of SE as they explain a substantial part of the overall genotypic variation. With increases in automation, growth processes may be monitored in the field on a daily basis or even multiple times per day. This will increase the precision in assessing genotype responses to the fluctuation in meteorological conditions and will allow quantification of the relationship of these responses to yield. Remote sensing by means of UAVs in combination with photogrammetric algorithms will allow measurement of these traits in breeding nurseries. We believe that this is paving the way for a more informed selection to climate adaptation within individual growing seasons.

Supplementary data
The following supplementary data are available at JXB online. Fig. S1. Correction of canopy height for spatial as well as random row and range effects.   .  Table S1. Genes of interest related to floral transition and flowering. Table S2. Chromosome-wise distance thresholds for LD decay <r 2 =0.2 Table S3. Corresponding marker-trait associations for final canopy height with respect to Zanke et al. (2016). Table S4. Selected putative candidate genes for heading GDD from the IWGSC reference genome annotation. Table S5. Selected putative candidate genes for GDD 95 from the IWGSC reference genome annotation. Table S6. Selected putative candidate genes for GDD 15 from the IWGSC reference genome annotation. Table S7. Three-year BLUEs of the investigated traits FH, heading GDD , GDD 15 , GDD 95 , GDD SE , time SE , slp SER~T , and int SER~T.