<figure><img src="image_1.png" /></figure>

HAL
open science

# Identification of QTLs affecting post-anthesis heat stress responses in European bread wheat

Gaëtan Touzy, Stéphane Lafarge, Elise Redondo, Vincent Lievin, Xavier Decoopman, Jacques Le Gouis, Sébastien Praud

► To cite this version:

Gaëtan Touzy, Stéphane Lafarge, Elise Redondo, Vincent Lievin, Xavier Decoopman, et al. Identification of QTLs affecting post-anthesis heat stress responses in European bread wheat. TAG Theoretical and Applied Genetics, 2022, 135 (3), pp.947-964. 10.1007/s00122-021-04008-5. hal-03524084

HAL Id: hal-03524084

https://hal.inrae.fr/hal-03524084

Submitted on 6 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Identification of QTLs affecting post-anthesis heat stress responses in European bread wheat

Gaëtan Touzy<fnref n="1" />,<fnref n="2" /> · Stéphane Lafarge<fnref n="2" /> · Elise Redondo<fnref n="2" /> · Vincent Lievin<fnref n="2" /> · Xavier Decoopman<fnref n="2" /> · Jacques Le Gouis<fnref n="3" /> · Sébastien Praud<fnref n="2" />

Received: 30 July 2021 / Accepted: 23 November 2021 / Published online: 5 January 2022
© The Author(s) 2021

## Abstract

**Key message** The response of a large panel of European elite wheat varieties to post-anthesis heat stress is influenced by 17 QTL linked to grain weight or the stay-green phenotype.

**Abstract** Heat stress is a critical abiotic stress for winter bread wheat (*Triticum aestivum* L.) especially at the flowering and grain filling stages, limiting its growth and productivity in Europe and elsewhere. The breeding of new high-yield and stress-tolerant wheat varieties requires improved understanding of the physiological and genetic bases of heat tolerance. To identify genomic areas associated with plant and grain characteristics under heat stress, a panel of elite European wheat varieties (N = 199) was evaluated under controlled conditions in 2016 and 2017. A split-plot design was used to test the effects of high temperature for ten days after flowering. Flowering time, leaf chlorophyll content, the number of productive spikes, grain number, grain weight and grain size were measured, and the senescence process was modeled. Using genotyping data from a 280 K SNP chip, a genome-wide association study was carried out to test the main effect of each SNP and the effect of SNP × treatment interaction. Genotype × treatment interactions were mainly observed for grain traits measured on the main shoots and tillers. We identified 10 QTLs associated with the main effect of at least one trait and seven QTLs associated with the response to post-anthesis heat stress. Of these, two main QTLs associated with the heat tolerance of thousand-kernel weight were identified on chromosomes 4B and 6B. These QTLs will be useful for breeders to improve grain yield in environments where terminal heat stress is likely to occur.

## Abbreviations

GWAS Genome-wide association study

QTL Quantitative trait loci

SNP Single nucleotide polymorphism

## Introduction

As one of the most important crops in the world, wheat (*Triticum aestivum* L.) has seen a steady growth in yields since the late 1960s in many areas (Calderini and Slafer 1998). However, a reversal of this trend has been observed in several countries over the last two decades (e.g., Peltonen-Sainio et al. 2009; Brisson et al. 2010; Rife et al. 2019). For example, lower yields in France have mainly been explained by climate variability causing drought stress during booting and high temperatures during grain filling (Brisson et al. 2010). It has been estimated on a world scale that climatic variations such as precipitation and heat stress may have led to a 5.5% reduction in grain yield since 1980, which corresponds to a loss of 35 M tonnes, equivalent to the annual wheat harvest of a country like France (Lobell et al. 2011). By the end of the twenty-first century, temperatures are expected to increase across most regions in the world, and on average by 0.3–1.7 °C in the lowest greenhouse gas emission scenario, and 2.6–4.8 °C in the highest one (Stocker et al. 2013). Combining climate predictions with a wheat

Communicated by Peter Langridge.

<figure><img src="image_3.png" /></figure>

Jacques Le Gouis
jacques.le-gouis@inrae.fr

<footnote n="1">¹ Arvalis-Institut du Végétal, Biopole Clermont Limagne, 63360 Saint-Beauzire, France</footnote>

<footnote n="2">² Centre de Recherche de Chappes, Route d'Ennezat CS90216, 63720 Chappes, France</footnote>

<footnote n="3">³ UMR 1095 GDEC, INRAE, Université Clermont Auvergne, Clermont-Ferrand, France</footnote>
simulation model, Semenov and Shewry (2011) showed that not drought, but the frequency and magnitude of heat stress around flowering time and during grain filling could significantly impact wheat yield, especially for wheat varieties that have not been specifically selected for their tolerance to heat stress, such as those commonly grown in Europe.

Heat stress scenarios can be characterized according to (1) the timing or the stage at which the plant is affected by stress, (2) the duration of stress, and (3) the intensity, i.e., the extremes of temperature to which the plant is exposed. Two kinds of stress scenarios are frequently associated with yield loss. “Heat shocks” are short periods (a couple of hours) of stress at very high temperatures (> 40 °C), and “chronic heat stress” is when average daily temperatures are higher than the optimum over a longer period. Crop responses to high temperature vary widely as they depend on the stress scenario, the germplasm (e.g., European wheat does not react in the same way as stress-adapted Australian wheat) and the experimental conditions (e.g., controlled or field trials). Porter and Gawith (1999) reviewed more than 65 papers describing the response of wheat to extreme temperatures and they determined that 20.7 (± 1.4) °C was the optimal temperature at the grain filling stage, while 35.4 (± 2.0) °C was the maximum temperature above which growth stops. To evaluate the impact of heat stress on plants, numerous studies have related plant behavior to the number of days during which the maximum temperature exceeded 25 °C at different growth stages (e.g., Sofield et al. 1977; Hunt et al. 1991; Gate 1995). Gate et al. (2010) applied this threshold to a multi-environment trial network in France and showed that for each day with a maximum temperature higher than 25 °C between heading and heading + 750 degree days (°Cd), the yield loss may be as much as 0.14 t ha⁻¹.

Heat stress impacts plant development and yield through numerous mechanisms. Firstly, higher temperatures during the crop cycle accelerate some biological processes (Johnson and Thornley 1985) such that physiological stages occur earlier (Porter and Gawith 1999). While earlier completion of some growth stages could be advantageous, wheat in Europe is still impacted by heat stress events during grain filling (Gate 2007; Gate et al. 2008). There are three categories of negative physiological impacts as reviewed by Cossani and Reynolds (2012). (1) Light interception is diminished with a reduction of the leaf area index and green area duration. E.g., Xu et al. (1995) showed that a reduction of leaf chlorophyll content during a post-anthesis heat stress led to earlier senescence. (2) Radiation use efficiency (RUE) is lessened due to altered protein structure, reducing the intensity of photosynthesis, decreasing CO₂ solubility and increasing photorespiration. (3) The partitioning of total assimilates is modified. If the stress occurs around the time of anthesis, the number of grains per spike may be impacted because

of pollen sterility or grain abortion (Tashiro and Wardlaw 1990). If the stress occurs after flowering, the remobilization of water-soluble carbohydrates and starch synthesis could be impacted with consequences on grain weight (Stone and Nicolas 1994; Wheeler et al. 1996; Calderini et al. 1999). High temperature during grain cell division might also reduce the final grain size (Nicolas et al. 1984; Commuri and Jones 2001; Barnabás et al. 2007).

Senescence is a natural process during which chlorophyll is catabolized and nutrients are remobilized from the source to sink organs of the plant (Taiz and Zeiger 2002). In monocarpic species, plants which senesce later and/or more slowly than a standard reference are described as “staying green.” Stay-green plants may be able to maintain photosynthesis during grain filling (Thomas and Howarth 2000). The stay-green phenomenon, when quantified, is an indicator of tolerance under heat and drought stresses (Olivares-Villegas et al. 2007; Pinto et al. 2010; Bouffier 2014). Different techniques have been used previously to monitor senescence and extract the main parameters (e.g., start, duration, rate and end) defining the stay-green trait, such as visual scoring (Lim et al. 2007), measurement of the chlorophyll content with a SPAD meter (Borrell et al. 2000; Christopher et al. 2008) and an estimation of the normalized difference vegetative index (e.g., Lopes and Reynolds 2012; Bouffier 2014; Christopher et al. 2016).

Several genetic studies have been conducted in heat stress conditions to identify quantitative trait loci (QTL) associated with pertinent agronomical and physiological traits, reviewed by Bouffier (2014), Farooq et al. (2014), and Tricker et al. (2018). Many of these studies were performed with bi-parental spring wheat populations growing in arid environments. However, evaluating the impact of heat stress alone, without simultaneous drought stress, is challenging in field trial situations. Field studies for QTL detection by CIMMYT in Mexico (Pinto et al. 2010, 2016; Bouffier et al. 2015) circumvented this difficulty by conducting experiments with an optimal and a late sowing date while ensuring both sets of plants were well irrigated. In effect, the delayed sowing exposed plants to terminal heat stress due to high summer temperatures, but they did not suffer from any drought stress. Performing field trials to screen varieties for heat tolerance is more complicated in Europe than at the CIMMYT sites. For instance, late sowing may increase the likelihood of heat stress, but it is not certain to occur each year in Europe’s temperate climate. Moreover, testing winter wheat panels that are composed of both photoperiod-sensitive and photoperiod-insensitive genotypes may bias the results. Vernalization and photoperiod requirements may cause varieties to flower at different times, such that escape and tolerance could easily be confused under heat stress. Most studies quantifying the impact of heat stress alone have indeed used controlled environments (Stone and

<figure><img src="image_4.png" /></figure>
Nicolas 1994; Gibson and Paulsen 1999; Spiertz et al. 2006; Shirdelmoghanloo et al. 2016; Telfer et al. 2021).

The development of new improved wheat varieties would be greatly facilitated if the genetic bases of heat tolerance were deciphered. The objective of this study was to quantify the effects of high temperatures during grain filling on elite European winter wheats. Seeking QTL from European germplasm for potential introgression through a European breeding program minimizes the risk of introducing negative QTL as the genetic background is closer. Our study was conducted in controlled conditions, and each plant was exposed to a heat treatment after it had flowered without any acclimation period. Flowering time, leaf chlorophyll content, the number of productive spikes, grain number, grain weight and grain size were measured. A GWAS was then performed to identify genomic areas associated with these traits and their response to heat treatment.

# Materials and methods

## Plant material and growth conditions

A total of 199 elite European winter wheat varieties, registered in Europe from 1974 to 2010, were assembled in a diversity panel (Online Resource 1). These varieties were evaluated in the greenhouse during the spring in 2016 and 2017 at the Biogemma site in Chappes, France.

Plants were sown on 22 February 2016 and 2 March 2017 in 1 cm-deep trays containing commercial potting substrate TS 3 from Klassman (Bourgoin Jallieu, France). The seedlings were raised in the greenhouse maintained at 21 °C daytime maximum and 15 °C night-time minimum. After two weeks, seedlings were vernalized for 49 to 56 days at 8 °C with an 8-h photoperiod. Then plants were transplanted into individual 1L pots (12 cm long and 11 cm wide). The rooting medium in the pots was commercial “Fleurissement 384” from Klassman, plus a controlled-release fertilizer “Nutricote type 70” from Fertil (Boulogne Billancourt, France), with 13% N, 13% P₂O₅ and 13% K₂O at 0.48 g per pot and with minor nutrients added before transplanting. Pots were watered daily and kept in trays containing water about 1 cm deep throughout the experiment to avoid water stress. Biocontrol methods and chemical pesticides were used to control diseases and pests. Air temperature and relative humidity were continuously monitored throughout the experiment at 15-min intervals using a temperature and humidity probe. Daytime maximum temperature of 21 °C was held for 12 h from 08:00 to 20:00. Similarly, night-time minimum temperature of 15 °C was held for 8 h from 22:00 to 06:00. The transition period between the daytime maximum and night-time minimum temperature was 2 h. Relative humidity was held at 60%. The photoperiod was 16 h

and sodium-vapor lamps (adjusted to be 100 cm away from the plant canopy) provided a photosynthetic complement if natural luminosity was less than 90 W m⁻².

## Split-plot design

A split-plot design was used. The first effect (whole-plot) was temperature regime, either non-stress (NS) or stress (S). The second effect (sub-plot) was the different varieties from the panel. Each treatment × variety combination was replicated three times. Fifteen tables (100 × 200 cm) were set out in a 3 × 5 matrix (Online Resource 2), and 120 pots were arranged in an 8 × 15 matrix on each table, which resulted in a density of 60 plants m⁻². Spring wheat varieties were positioned all around the tables to minimize border effects. To compensate for the known temperature and luminosity gradients in the greenhouse, the three replicates were arranged in an East/West direction, i.e., replicate 1 on Tables 1, 4, 7, 10 and 13; replicate 2 on Tables 2, 5, 8, 11 and 14; and replicate 3 on Tables 3, 6, 9, 12 and 15. Each treatment occupied two and a half tables, e.g., for replicate 1, Tables 1, 4 and half of Table 7 were assigned to one treatment, while Tables 10, 13 and the other half of Table 7 were assigned to another treatment (Online Resource 2). The tables were split into two, defining five sub-blocks per whole-plot, where three replicate checks were placed at random (checks were the registered lines Bermude, Boregar and Premio). The relative location of temperature treatments was switched between the 2016 and 2017 experiments and different randomizations were drawn up for genotypes and check locations.

Forty thermo buttons from the Proges Plus company (Willems, France) were used, two per table, to monitor temperature every half hour during the growth cycle. The *krige* function from the package “gstat” (Pebesma 2004) was used to perform an ordinary kriging operation and assign a mean temperature value for each pot.

## Heat treatment

Three days after anthesis (DAA), plants in the S temperature regime were put into another greenhouse at a temperature of 29 °C day/23 °C night but with the same photoperiod for ten days. The plants with the same flowering date were arranged on the same table so relative positions during this time were different from the original plot. However, after the heat treatment, the plants were returned to their former location.

The quality of temperature control was assessed by thermo buttons. The average daytime temperatures (between 07:00 and 21:00) and night-time temperatures (between 21:00 and 07:00) in the 2016 experiment were, respectively, 30.4 ± 2.9 °C and 24.4 ± 2.6 °C for S and 20.0 °C ± 2.6 and 16.4 °C ± 2.9 for NS, and in the 2017 experiment were 29.0 ± 3.0 °C and 23.8 ± 3.5 °C for S and 21.7° C ± 1.7 and 17.6
± 2.3 °C for NS. Relative humidity during the day and night was 49.2 ± 9.2% and 61.0 ± 8.4% for S treatment and 74.0 ± 10.3% and 80.2 ± 8.3% for NS treatment in 2016; and was 60.0 ± 9.6% and 70.3 ± 9.0% for S treatment and 77.1 ± 8.0% and 82.8 ± 6.3% for NS treatment in 2017. Vapor pressure deficit (VPD) was calculated from temperature and humidity values. VPD during the day and night, respectively, was 2.13 ± 0.38 kPa and 1.86 ± 0.23 kPa for the S treatment and 1.71 ± 0.21 kPa and 1.48 ± 0.15 kPa for the NS treatment in 2016, and 2.40 ± 0.37 kPa and 2.07 ± 0.32 kPa for the S treatment and 1.99 ± 0.16 kPa and 1.67 ± 0.16 kPa for the NS treatment in 2017.

## Measurement of yield components

For each pot, flowering time (D.Z65) corresponding to Zadoks stage 65 (Zadoks et al. 1974) was scored daily on the main shoot. At the end of the experiments, the number of fertile spikes per plant (SPP) was counted, then using a small plot harvester grains were harvested separately from the main shoot and the other shoots. Grains were dried for one week at 35 °C to ensure a homogenous moisture content before threshing. Using a MARVIN digital seed analyzer (MARViTECH, Wittenburg, Germany), the number, width, length and area of dry grains were determined for grain from both the main shoot and the tillers. The grains were weighed (GW), and the thousand kernel weights (TKW) were calculated.

## Measurement of chlorophyll content and estimation of senescence traits

A Dualex Scientific™ sensor (ForceA, Orsay, France) was used to measure chlorophyll (CHL) content (µg cm⁻³) from light transmission as follows: $CHL = \frac{\text{Near infrared transmission} - \text{Red transmission}}{\text{Red transmission}}$. CHL was measured weekly, a minimum of five different dates for each plant from flowering time until maturity. Measurements were taken at five different points along the flag leaf blade.

Senescence was modeled for each plant using the following logistic function:

$$CHL = D + \frac{K}{1 + a \times \exp^{(r \times t)}} \quad (1)$$

where CHL is a vector of the chlorophyll content (µg cm⁻³) from the flowering date to maturity, D is the final CHL of the dead plant, K corresponds to the difference in CHL between the maximum and final values, a and r are indicators of the senescence rate, and t is the thermal time from plant anthesis. To enable proper data fits, an extra data point corresponding to the CHL value of the fully senesced plant was added 200 °C day⁻¹ after the last measurement.

Fitted logistic functions were obtained using the nls function (Bates and Chambers 1992) from the stat package (R Development Core Team 2011). A first run was performed on raw data to identify outliers. To do this, the residuals distribution was compared to a normal distribution and outliers were identified using method I of the getOutliers function (Van Der Loo 2010) from the extremevalues package (R Development Core Team, 2011). A second run was then performed without the outliers to estimate the different parameters of the logistic model. Pseudo R² were calculated to evaluate the strength of fitted models according to Efron's (1978) formula:

$$R^2 = 1 - \frac{\sum_{i=1}^{N} (CHL_i^2 - \overline{CHL_i})^2}{\sum_{i=1}^{N} [(CHL_i^2 - \overline{CHL})^2]} \quad (2)$$

where N is the number of observations in the model, CHL is a vector of the chlorophyll content measurements (µg cm⁻³) from flowering to maturity, CHL is the mean of the CHL values, and CHL is the value predicted by the model.

The logistic function in Eq. 1 provided a close fit to the experimental data (e.g., Fig. 1). The model converged on more than 95% of the plants in both experiments. A second data cleaning operation was performed by hand to ensure that the adjustment of the kinetics corresponds correctly to the senescence process; this removal of outliers was based on senescence parameters distributions, visual curve fitness and R². Overall, it was not possible to satisfactorily model senescence data for 199 plants (7.8%). The mean of R² for

<figure><img src="image_5.png" /><figcaption>Fig. 1 Example of the fitted logistic curve for modeling senescence and estimating stay-green traits. The red curve is an example of the fitted logistic curve for plant 455 in the 2016 experiment; the black open dots represent the raw data; the red crosses represent outliers; green and blue dashed lines are, respectively, the y-axis and x-axis coordinates of the start (START), the point of inflection (PI) and the end (END) of senescence; the dashed gray curve is the third derivative of the function which at its minima defines the start and end of senescence and at its maximum defines the x-axis of the point of inflection; the orange arrow indicates the maximal senescence rate; and the gray shading represents the area under the senescence curve (AUC)</figcaption></figure>
the remaining plants was 0.94 (± 0.05) and 0.98 (± 0.02) for the 2016 and 2017 experiments, respectively. All the senescence parameters presented characteristics of normal distribution except for a.

Additional senescence traits were estimated using the properties of the logistic function (Fig. 1). The x-axis and y-axis of the point of inflection were obtained as $XPI = \frac{|\ln(a)|}{r}$ and $YPI = \frac{K}{2}$. The maximal rate of senescence was obtained with the first derivative as VIT = $f'(XPI)$. The start (START) and end (END) of senescence were determined as the thermal time when the second derivative of the function attained its minimums. The chlorophyll content at these two points was noted as fSTART and fEND. We calculated the area under the senescence curve from anthesis until the completion of senescence (AUCK). It was calculated as the integral from anthesis to the end of senescence for the x-axis, and from K to 0 for the y-axis. AUCK is related to the start, the speed of senescence and the CHL maximum. We also calculated the area under the curve without the CHL maximal, named AUC. AUC is represented in Fig. 1 as the gray area, and can be considered to express the stay-green trait as it combines both the start and rate of senescence.

## Statistical analysis

The same experimental designs were used in both experiments (2016 and 2017), so we used the same baseline model to estimate genotype × treatment interaction (G × T). This model was specific to the split-plot design of our study and was adapted to a multi-year analysis and SNP by treatment interaction. The models were fitted using a mixed model written in R using the ASReml-R package (Butler et al. 2009).

The baseline model for the G × T estimate in one experiment was written as follows:

$$Y_{ijkl} = \mu + R_j + T_k + T_k \times R_j + B_l(R_j T_k) + G_i + G_i \times T_k + \varepsilon_{ijkl} \quad (3)$$

where $Y_{ijkl}$ is the phenotypic value of genotype i (i = 1:199) in the sub-block l of the treatment k of the replicate j, $\mu$ is the general mean, $R_j$ is the fixed effect of the replicate j (j = 1:3), $T_k$ is the fixed effect of treatment k with k = 1 for non-stress treatment (NS) and k = 2 for stress treatment (S), $T_k \times R_j$ is the fixed effect of treatment k within replicate j, $B_l(R_j T_k) \sim N(0, \sigma_l^2)$ is the random effect of the incomplete sub-block l (l = 1:5) within treatment k of replicate j, $G_i \sim N(0, \sigma_g^2)$ K is the random genotype effect where K is the genomic relationship matrix calculated using the formula from VanRaden (2008) (see below), $G_i \times T_k$ is the random genotype × treatment interaction effect with $G_i \times T_k \sim N(0, \sigma_{g,t}^2 \cdot K \circ I_{nT})$, where $K \circ I_{nT}$ is the Hadamar product between the K matrix and the $I_{nT}$, identity matrix and $\varepsilon_{ijkl}$ is the residual error $\varepsilon \sim N(0, \sigma_{\varepsilon}^2)$.

The random polygenic effect (K) was modeled with a kinship matrix as a variance/covariance matrix. Pairwise kinship coefficients were computed according to the first method described by VanRaden (2008), with a subset of 27,680 SNP obtained by removing highly correlated SNP. In brief, we computed $r^2$ for each pair of SNP located on the same chromosome. A dendrogram was built with this matrix of $r^2$ and one SNP was randomly selected from each cluster with a distance < 0.1 (where distance = $1-r^2$).

Broad sense heritability ($H^2$) was estimated as follows:

$$H^2 = \frac{\sigma_g^2}{\sigma_g^2 + \frac{\sigma_{g,t}^2}{t} + \frac{\sigma_{\varepsilon}^2}{tn}} \quad (4)$$

where $\sigma_g^2$, $\sigma_{g,t}^2$, and $\sigma_{\varepsilon}^2$ are genotype, genotype by treatment interaction and residual variance, and t and n are the number of treatments and replicates.

For a combined analysis across the two years, the baseline model was enhanced with a year effect to give:

$$Y_{ijklm} = \mu + A_m + T_k + T_k \times A_m + R_j(A_m) + T_k(A_m R_j) + B_l(A_m R_j T_k) + G_i + G_i \times T_k + G_i \times A_m + \varepsilon_{ijklm} \quad (5)$$

where the new terms $A_m$ and $G_i \times A_m \sim N(0, \sigma_{g,a}^2 \cdot K \circ I_{nA})$ are the effects of the mth year and of its interaction with genotype, respectively. In this model, we must stress that part of the $G_i \times T_k \times A_m$ interaction was included in the model residual, resulting in an underestimation of the specific influence of $G_i \times T_k$ and $G_i \times A_m$.

Broad sense heritability ($H^2$) was estimated as follows:

$$H = \frac{\sigma_g^2}{\sigma_g^2 + \frac{\sigma_{g,t}^2}{t} + \frac{\sigma_{g,a}^2}{a} + \frac{\sigma_{\varepsilon}^2}{atn}} \quad (6)$$

where $\sigma_g^2$, $\sigma_{g,t}^2$, and $\sigma_{\varepsilon}^2$ are genotype, genotype by treatment interaction and residual variance, respectively, and t, a and n are the number of years, treatments and replicates.

To test the stress impact and relationship between traits the genomic best linear unbiased predictions (gBLUP) for each genotype and treatment were extracted from Eqs. 3 or 5 as follows:

$$\overline{GGT}_{ik} = \hat{\mu} + \hat{G}_i + \hat{G}_i \times \overline{T}_k \quad (7)$$

Based on the $\overline{GGT}_{ik}$ the relationships between traits were estimated with the Pearson's correlation and principal component analysis (PCA).

The stress index was calculated as follows:

$$SI = 100 \times (\overline{GGT}_{.2} - \overline{GGT}_{.1}) / \overline{GGT}_{.2} \quad (8)$$
where $\overline{GGT}_{.1}$ is the average of gBLUP in stress treatment (S) and $\overline{GGT}_{.2}$ is the average of gBLUP in no-stress treatment (NS). This formula was calculated for each genotype to give the stress intensity at genotype level, and based on this index the correlation between traits was also calculated.

Finally, the baseline model was adapted to multi-year GWAS split-plot analysis as follows:

$$Y_{ijklm} = \mu + A_m + T_k + T_k(A_m) + \alpha_i + \alpha_i \times T_k + A_m(R_j) + T_k(A_m R_j) + B_l(A_m R_j T_k) + G'_i + G'_i \times T_k + G'_i \times A_m + \varepsilon_{ijklm} \quad (9)$$

where the new terms $\alpha_i + \alpha_i \times T_k$ represent the allelic fixed effect of the genotype i at SNP $\alpha_i$ and the interaction between marker $\alpha_i$ and treatment k, respectively. As part of the genetic effect is captured by the SNP, the notation for the genetic background is $G'_i$, so it is defined in the same way as the genetic effect for Eq. 5.

All the varieties were genotyped with the Affymetrix Axiom 280 K SNP array (Rimbert et al. 2018). Only the Polymorphic High-Resolution SNP were used in this analysis. All SNP were physically mapped on the genome reference sequence RefSeq V1.0 (IWGSC et al. 2018). Heterozygous loci were considered as missing data. SNP with more than 10% missing data, and SNP with a minor allele frequency (MAF) below 5% were discarded. Missing data were imputed with Beagle 4.1 (Browning and Browning 2016). Finally, 164,198 SNP were tested with Eq. 9.

To avoid proximal contamination, SNP were tested using a K matrix computed with SNP that were not located on the same chromosome as the SNP tested (Rincent et al. 2014).

The significance of random factors was tested one by one using the likelihood ratio test (LRT) (Kendall and Stuart 1979), based on log-likelihood (Lmax) differences between the complete and the reduced model without the test factor. A Wald test was performed on the complete model for testing the SNP main effect and the SNP by treatment interaction. Following the method of Gao et al. (2009), SNP were considered to be significantly associated with a trait if the -log10(P-value) exceeded a threshold of 5.26. The estimation of the part of the variance explained by the SNP effect on the treatment was calculated using the $R^2$ resulting from the regression of the $Gi \times Tk$ coefficients depending on these effects using Eq. 5.

QTL confidence interval boundaries from GWAS results were defined following the method described in Cormier et al. (2014). Briefly, SNP belonging to the same linkage disequilibrium (LD) cluster were defined as a group of quantitative trait nucleotides. Clustering was performed by averaging $r^2$ distances and the tree was cut at 1 “critical LD unit” (critical LD = 0.24, Breseghello and Sorrells 2006). QTL

boundaries were defined as the maximal and minimal map positions of the SNP within the quantitative trait nucleotides extended by the LD decay specific to the genomic area of the QTL.

## Results

### Characterization of growth conditions and modeling senescence

A panel of 199 elite European wheat varieties were grown in controlled greenhouse conditions in a split-plot design with or without a 10-day post-anthesis heat stress in the years 2016 and 2017. NS temperatures were regulated at 21 °C daytime maximum and 15 °C night-time minimum and S temperatures at 29 °C daytime maximum and 23 °C night-time minimum. According to multiple well-spaced thermobutton measurements, on average thermal time was 27.8 °C day$^{-1}$ and 27.1 °C day$^{-1}$ in the S treatment and 19 °C day$^{-1}$ and 20.6 °C day$^{-1}$ in the NS treatment, for 2016 and 2017, respectively. VPD was also calculated from the temperature and relative humidity measurements (see "Materials and methods") and though slightly higher in 2016 than in 2017, in both years of the experiment there was a clear difference in VPD between S and NS conditions. From the regular monitoring of conditions and measurement of phenotypes of plants (including the check variety) and yield components, it was possible to analyze the effect of heat stress on the panel in the absence of drought stress, because plants were well-watered throughout.

We used some of the non-invasive physical measurements to model the wheat senescence process and thus quantify the behavior of the wheat genotypes in a biologically meaningful way. Here, we modeled senescence in terms of chlorophyll content (estimated from light transmission) and thermal time using the logistic function in Eq. 1 (see "Materials and methods"). The model provided a close fit to the experimental data, as shown for one example of one genotype in one experiment in Fig. 1. The model converged on more than 95% of the plants in both years of the experiment. The second data cleaning step were found to be essential to ensure the shape of the curves corresponded to a normal biological process, for example, the total chlorophyll content at the end of senescence could have been estimated at around 40 µg cm$^{-3}$ or -10 µg cm$^{-3}$, which does not reflect the real distribution for this trait. For this reason, it was not possible to model the senescence of 199 plants (7.8%). The mean of $R^2$ for the remaining plants was 0.94 (± 0.05) and 0.98 (± 0.02) for the 2016 and 2017 experiments, respectively. Additional senescence traits were estimated using the properties of the logistic function, essentially the start, the end and the rate of senescence (Fig. 1). In practice, the stay-green
trait can be quantified as AUC (the gray area represented in Fig. 1) as this measure combines the start and rate of senescence.

From the regular monitoring of conditions and comprehensive measurement of phenotypes of plants (including the check variety) and yield components, it was possible to closely model the senescence process and analyze to what extent the genotype of wheat affected the observed and modeled phenotypes (Online Resource 3).

## Phenotypic analysis and variance decomposition of heat stress in wheat

Significant genotypic effects were observed for all traits except for senescence parameters *D* and *a* in both years, according to the variance components associated with the different effects extracted from the model in Eq. 3 for 2016 and 2017 experiments (Online Resource 3). Traits with significant genotypic effects showed highly variable heritability, ranging from 0.20 for fEND in 2017 to 0.95 for D.Z65 in 2016 and 2017.

As the heat treatment was applied to each plant 3 DAA, as expected, no significant impact was observed on flowering time in the 2016 and 2017 experiments. However, a slight effect of heat stress was observed on the number of fertile spikes per plant: the stress index (SI) was 3.14% (not significant, ns) in 2016 and 4.85% (*P*-value < 0.05) in 2017. The number of grains on the main shoot was not affected as the SI was 7.15% (0.05 < *P*-value < 0.01; Table 1), and on tillers there was no effect in 2016 with an SI of 8.26% (ns) and a strong effect in 2017 with an SI of 13.92% (*P*-value < 0.01; Online Resource 3). S effects on tillers were expected because they flower later than the main shoot, during or after the heat stress.

Considering the senescence traits, only the *x*-axis of the point of inflection (XPI) and the end of senescence (END) in 2016 were significantly affected by the S conditions with SI of 6.89% and 6.44%, respectively. Heat stress impacted the area, length and width of grain from both the main shoot and the tillers in both years. For the TKW of the grain from the main shoot, the means were 46.77 g in 2016 and 41.96 g in 2017 for the NS treatment, but 34.88 g and 32.86 g for the S treatment (Fig. 2). S reduced the TKW of the main shoot by 25.73% in 2016 and 21.88% in 2017. The performance of the varieties and their SI are available in Online Resource 4.

Variance decomposition revealed significant G×T interactions for TKW, grain area, grain width and grain length for both the main shoot and the tillers in 2016 or 2017 (Online Resource 3). The ratio of the G×T variance on the sum of the genetic and G × T variances was the highest for TKW and grain width on the main shoot. For TKW, this ratio was

22.45% in 2016 and 22.76% in 2017. For other traits, the G×T interaction varied between the two years.

## Multi-year phenotypic analysis and variance decomposition

Variance components associated with different effects were extracted from the model (Eq. 5), combining the 2016 and 2017 experiments, and are presented in Table 1. Significant genotypic effects were observed for all traits except for senescence parameters *D*, *a* and *r*, and fEND. Traits with significant genotypic effects showed highly variable heritability ranging from 0.33 for the rate of senescence (VIT) to 0.88 for the flowering date (D.Z65).

No significant effect of treatment was observed for D.Z65, SPP and senescence parameters. Only slight effects were observed for XPI, START and AUC with SI of 3.75%, 5.39% and 3.75%, respectively. The shapes of the mean senescence curves were similar but in S conditions senescence started earlier (Fig. 3). A mildly significant stress effect was observed for the number of grains on the main shoot and on tillers, the SI being 7.15% and 11.24%, respectively. For other traits such as TKW and grain weight, area, width and length, a significant effect was observed with an SI varying from 4.16% for the length of grain from tillers to 28.74% for the weight of grain from the main shoot.

The year effect was most significant for senescence traits (Table 1), with only a weak effect on other traits. Average TKW was higher in 2016 than in 2017 (Fig. 2).

Looking at the interactions, we found a significant G×Y interaction for all traits, except for senescence parameter *a*. The genetic variance of the senescence rate (VIT) shows that this trait has the strongest G×Y interaction. The grain parameters showed significant high G×T effects.

For TKW, the variance explained by G×T and G×Y interactions was similar for the main shoot at 4.24 g² and 4.39 g², and slightly lower for the tillers with 2.93 g² and 3.21 g², respectively. This strong G×Y effect could be explained by uncontrolled environmental factors such as global radiation, which was far lower in 2016 than in 2017, even when supplemented with sodium-vapor lamps.

## An overview of senescence traits

The variance components of senescence traits are detailed in Table 1. The change in heritability and the response of senescence to the temperature treatment was reported for the whole grain filling period. Figure 3 compares the senescence curve variability based on raw data per plant and based on gBLUP of the senescence parameters. For raw data, the size of the confidence interval (CI) around the mean (± one standard deviation) remained stable during senescence until around 400 °C day when the CI widened considerably
<table><caption>**Table 1** Variance components for traits measured on 199 wheat varieties grown for two years with two post-flowering temperature treatments</caption><thead><tr><th>Variable</th><th>Units</th><th>Mean</th><th>Sd</th><th>H<sup>2</sup></th><th>SI</th><th>ssY</th><th>ssT</th><th>ssYT</th><th>varG</th><th>varGT</th><th>varGY</th><th>varRes</th></tr></thead><tbody><tr><td>Flowering Date</td><td>D</td><td>175.7</td><td>6.84</td><td>0.88</td><td>-0.06</td><td>32.17</td><td>0.31</td><td>0.57</td><td>24.29***</td><td>0.14</td><td>5.03***</td><td>9.04</td></tr><tr><td>Senescence Parameters</td><td>D</td><td>5.10</td><td>1.09</td><td>0.05</td><td>1.66</td><td>118.08***</td><td>0.4</td><td>4.09</td><td>0.04</td><td>0.25**</td><td>0.25*</td><td>6.27</td></tr><tr><td></td><td>K</td><td>46.14</td><td>3.16</td><td>0.77</td><td>-0.27</td><td>15.5</td><td>2.74</td><td>40.45</td><td>7.62***</td><td>2.07E-06</td><td>1.62***</td><td>17.44</td></tr><tr><td></td><td>A</td><td>6.15E-04</td><td>4.90E-04</td><td>4.31E-03</td><td>-29.22</td><td>1.08E-04*</td><td>4.35E-06</td><td>1.99E-06</td><td>7.68E-09</td><td>1.68E-07</td><td>1.07E-07</td><td>1.95E-05</td></tr><tr><td></td><td>R</td><td>2.28E-02</td><td>4.94E-03</td><td>0.18</td><td>-1.67</td><td>1.36E-03***</td><td>4.05E-06</td><td>8.36E-05</td><td>2.16E-06</td><td>2.71E-06**</td><td>5.19E-06***</td><td>6.70E-05</td></tr><tr><td>Senescence traits</td><td>XPI</td><td>°Cd</td><td>657</td><td>82.54</td><td>0.81</td><td>3.75</td><td>164.647***</td><td>25.379</td><td>20.941</td><td>3778***</td><td>143</td><td>508***</td><td>6675</td></tr><tr><td></td><td>YPI</td><td>µg cm<sup>-3</sup></td><td>28.2</td><td>2.08</td><td>0.74</td><td>0.11</td><td>67.25***</td><td>0.03</td><td>0.46</td><td>2.09***</td><td>0.13*</td><td>0.62***</td><td>4.32</td></tr><tr><td></td><td>VIT</td><td>µg cm<sup>-3</sup>, °Cd<sup>-1</sup></td><td>-0.26</td><td>0.06</td><td>0.33</td><td>-2.28</td><td>0.09***</td><td>8.18E-04</td><td>0.01</td><td>6.47E-04**</td><td>3.69E-04***</td><td>9.31E-04***</td><td>0.01</td></tr><tr><td></td><td>START</td><td>°Cd</td><td>532</td><td>67.39</td><td>0.74</td><td>5.39</td><td>18.764</td><td>28.689*</td><td>6276</td><td>2897***</td><td>221</td><td>665***</td><td>6819</td></tr><tr><td></td><td>START</td><td>µg cm<sup>-3</sup></td><td>47.0</td><td>3.28</td><td>0.81</td><td>-0.04</td><td>64.1**</td><td>0.03</td><td>5.3</td><td>6.58***</td><td>0.06</td><td>1.57***</td><td>9.13</td></tr><tr><td></td><td>END</td><td>°Cd</td><td>781</td><td>107.61</td><td>0.74</td><td>2.62</td><td>350.664***</td><td>23.832</td><td>63.877*</td><td>4814***</td><td>226</td><td>1051***</td><td>12.938</td></tr><tr><td></td><td>END</td><td>µg cm<sup>-3</sup></td><td>9.3</td><td>1.16</td><td>0.14</td><td>0.81</td><td>121.98***</td><td>0.3</td><td>1.9</td><td>0.11</td><td>0.22**</td><td>0.27**</td><td>5.26</td></tr><tr><td></td><td>AUC</td><td>µg cm<sup>-3</sup>°Cd</td><td>26.828</td><td>4176</td><td>0.79</td><td>3.75</td><td>244.029.819***</td><td>48.250.209</td><td>13.692.149</td><td>9.649.426***</td><td>247.761</td><td>2.040.704***</td><td>16.651.766</td></tr><tr><td></td><td>AUCK</td><td>µg cm<sup>-3</sup>°Cd</td><td>34.148</td><td>5746</td><td>0.81</td><td>3.72</td><td>582.017.317***</td><td>54.314.087</td><td>29.880.570</td><td>15.481.554***</td><td>629.865*</td><td>2.995.929***</td><td>20.730.987</td></tr><tr><td>Number of Spikes</td><td>SPP</td><td>plant<sup>-1</sup></td><td>4.8</td><td>1.05</td><td>0.70</td><td>4.01</td><td>2.67</td><td>2.07</td><td>1.07</td><td>0.64***</td><td>0.06*</td><td>0.21***</td><td>1.67</td></tr><tr><td>Main Shoot</td><td>NbGrainS</td><td>39.5</td><td>8.11</td><td>0.80</td><td>7.15</td><td>184.54†</td><td>190.65†</td><td>65.76</td><td>39.5***</td><td>4.29***</td><td>5.48***</td><td>57.16</td></tr><tr><td></td><td>GW</td><td>g</td><td>1.54</td><td>0.45</td><td>0.72</td><td>28.74</td><td>0.49*</td><td>3.3***</td><td>0.19</td><td>0.06***</td><td>0.01***</td><td>0.01***</td><td>0.11</td></tr><tr><td></td><td>TKW</td><td>g</td><td>39.1</td><td>7.29</td><td>0.64</td><td>23.90</td><td>166.16*</td><td>1582.09***</td><td>141.67*</td><td>12.3***</td><td>4.24***</td><td>4.39***</td><td>32.20</td></tr><tr><td></td><td>GrainArea</td><td>mm<sup>2</sup></td><td>14.5</td><td>1.65</td><td>0.72</td><td>14.40</td><td>6.53*</td><td>70.3***</td><td>2.84</td><td>0.85***</td><td>0.2***</td><td>0.2***</td><td>1.53</td></tr><tr><td></td><td>GrainWidth</td><td>mm</td><td>3.42</td><td>0.24</td><td>0.66</td><td>10.01</td><td>0.07</td><td>1.83***</td><td>0.16*</td><td>0.01***</td><td>4.62E-03***</td><td>4.99E-03***</td><td>0.03</td></tr><tr><td></td><td>GrainLength</td><td>mm</td><td>5.84</td><td>0.29</td><td>0.85</td><td>4.53</td><td>0.26*</td><td>1.24***</td><td>0.01</td><td>0.04***</td><td>4.54E-03***</td><td>4.49E-03***</td><td>0.04</td></tr><tr><td></td><td>NbGrainS</td><td>138</td><td>35.38</td><td>0.71</td><td>11.24</td><td>162.56</td><td>8587.89*</td><td>3586.11</td><td>765.63***</td><td>111.84***</td><td>207.08***</td><td>1813.77</td></tr><tr><td></td><td>GW</td><td>g</td><td>5.37</td><td>1.64</td><td>0.68</td><td>26.24</td><td>0.62</td><td>57.14***</td><td>0.36</td><td>1.14***</td><td>0.25***</td><td>0.34***</td><td>2.97</td></tr><tr><td></td><td>TKW</td><td>g</td><td>38.9</td><td>6.04</td><td>0.69</td><td>17.54</td><td>27.79</td><td>676.46***</td><td>96.72*</td><td>10.9***</td><td>2.93***</td><td>3.21***</td><td>20.80</td></tr><tr><td></td><td>GrainArea</td><td>mm<sup>2</sup></td><td>14.0</td><td>1.43</td><td>0.78</td><td>10.68</td><td>4.25*</td><td>30.97***</td><td>6.59*</td><td>0.85***</td><td>0.15***</td><td>0.13***</td><td>1.10</td></tr><tr><td></td><td>GrainWidth</td><td>mm</td><td>3.37</td><td>0.20</td><td>0.75</td><td>6.65</td><td>0.03</td><td>0.68***</td><td>0.23**</td><td>0.01***</td><td>3.36E-03***</td><td>2.24E-03***</td><td>0.02</td></tr><tr><td></td><td>GrainLength</td><td>mm</td><td>5.71</td><td>0.29</td><td>0.86</td><td>4.16</td><td>0.19*</td><td>0.94***</td><td>0.07</td><td>0.05***</td><td>3.93E-03***</td><td>5.48E-03***</td><td>0.03</td></tr></tbody></table>

Abbreviations mean, mean for fitted values; Sd, standard deviation for fitted values; H², broad sense heritability; SI, stress index; varT, variance explained by the treatment (T); varG, variance explained by the genotype (G); varG×T, variance explained by the G×T interaction; varG×Y, variance explained by the G×Year interaction; varRes, residual variance. Wald and LTR test: **P<0.001; *P<0.01; †P<0.05; †P<0.1, non-significant P>0.05**

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 2** Boxplot of TKW for 199 wheat cultivars grown over two years (2016 and 2017) with two post-flowering temperature treatments (NS no stress, S stress). gBLUP values calculated from Eq. 2 were used (see "Materials and methods"). Quartiles and medians are used to construct the boxes. Whiskers extend to 1.5 times the interquartile range from the box. The means of each sample are significantly different from each other (α=0.001) based on Tukey HSD test</figcaption></figure>

(Fig. 3a). By contrast, for the gBLUP of the senescence parameters (Fig. 3b), the CI narrowed at the end of senescence. In fact, varieties are declared as a random effect in the mixed model, so if there are no significant differences among genotypes the shrinkage factor will tend toward infinity and breeding value predictions will tend to the mean (low heritability). Figure 3 shows that while there are differences between genotypes at the beginning of senescence, all the genotypes converged to values which are not discriminating. Furthermore, as the figure summarizes the distribution of senescence kinetics for all genotypes in both treatments, we can see that the heat treatment only had a slight effect on the senescence curves (the SI is at 5.39 for START and 3.75 for XPI).

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Senescence kinetics over thermal time after anthesis for 199 wheat varieties grown with two post-flowering temperature treatments. The mean (curves) ± standard deviation (shaded areas) of the fitted logistic senescence curve calculated from a raw data for both</figcaption></figure>

## Genetics correlation

Relationships among traits and genotypes for each heat treatment were investigated using Pearson's correlation and biplot analysis with gBLUP estimated per treatment from Eq. 7. All correlation data are presented in Online Resource 5. Figure 4 shows the PCA of a subset of traits related to the characteristics of grain from the main shoot, senescence traits, flowering time and number of fertile spikes per plant where vector length shows the extent of variation explained by each trait. The first two axes explained up to 62% of total variability in NS (Fig. 4a) and 64% in S (Fig. 4b) treatments, respectively. The relationships between traits is similar for both treatments. Two main clusters of traits are discernible: (1) traits related to the grain such as TKW, and grain area, length and width, and (2) senescence traits. No significant correlation was detected between D.Z65 and TKW in the NS treatment (r=0.08, p>0.1), but a significant weak correlation was detected between these two traits in the S treatment (r=-0.24, p<0.001). We also detected significant weak correlations (p<0.05) between TKW and the respective senescence traits of YPI (r=0.19), fSTART (r=0.20), END (r=-0.16) and VIT (r=0.19) in the NS treatment. In the S treatment, we still detected correlations between TKW and, respectively, YPI (r=0.26), fSTART (r=0.27), and VIT (r=-0.18) but not between TKW and END, while a positive correlation was detected with AUC (r=0.19).

## Genome-wide association study

A single-locus genome-wide association study, involving a total of 164,198 SNP, was performed using raw data from both years. The main effect and the SNP×treatment interaction were tested using Eq. 9. All phenotypic traits were tested, except the senescence parameters *a* and *D*, and fEND which had not shown significant genotypic effects. Four genomic

<figure><img src="image_9.png" /></figure>

years or *b* from gBLUP parameters are shown for all genotypes in stress (red) and no stress (blue) treatments. The red double-headed arrow represents the period during which the stressed plants were exposed to heat treatment
<figure><img src="image_10.png" /><figcaption>**Fig. 4** Biplots of the principal components of a subset of traits (lines) measured on 199 wheat varieties (dots) grown in **a** no stress and **b** stress post-flowering temperature treatments. The varieties with the</figcaption></figure>

<figure><img src="image_11.png" /></figure>

five highest and five lowest main shoot TKW values are, respectively, plotted as green and red dots with their names given. For details of trait names see "Materials and methods"

<figure><img src="image_12.png" /><figcaption>**Fig. 5** Manhattan plot of the GWAS of the SNP main effect and the SNP × treatment interaction for flowering date (D.Z65) and TKW measured on 199 wheat varieties grown in stress or non-stress post-flowering temperature treatments. **a** SNP main effect for D.Z65; **b**</figcaption></figure>

<figure><img src="image_13.png" /></figure>

SNP by treatment interaction for D.Z65; **d** SNP main effect for TKW; **c** SNP by treatment interaction for TKW. The red line represents the 5.6 significant threshold calculated as described by Gao et al. (2015)

<figure><img src="image_14.png" /></figure>
<table><caption>**Table 2** List of associated SNP following a GWAS conducted on 199 wheat varieties grown for two years with two post-flowering temperature treatments: SNP with no interaction effect</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Trait</th><th rowspan="2">Units</th><th rowspan="2">Chr</th><th rowspan="2">Position (bp)</th><th rowspan="2">SNP</th><th rowspan="2">Major Allele</th><th rowspan="2">Freq</th><th colspan="2">SNP main effect</th><th colspan="2">QTL boundary</th></tr><tr><th>LOD</th><th>Major allele effect</th><th>Effect SD</th><th>QTL</th><th>Nb. SNP</th><th>Min. (bp)</th><th>Max. (bp)</th></tr></thead><tbody><tr><td rowspan="3">Flowering Date</td><td>D.Z65</td><td>d</td><td>4A</td><td>576,159,104</td><td>AX-89667446</td><td>A</td><td>0.64</td><td>5.62</td><td>3.98</td><td>2.35</td><td>QTL.05</td><td>13</td><td>575,485,603</td><td>581,337,826</td></tr><tr><td>D.Z65</td><td>d</td><td>5A</td><td>584,428,207</td><td>AX-89390208</td><td>C</td><td>0.85</td><td>6.93</td><td>6.16</td><td>2.42</td><td>QTL.07</td><td>2</td><td>583,039,479</td><td>585,816,937</td></tr><tr><td>D.Z65</td><td>d</td><td>6A</td><td>13,010,065</td><td>AX-89584501</td><td>G</td><td>0.55</td><td>5.39</td><td>3.63</td><td>2.34</td><td>QTL.08</td><td>1</td><td>10,570,285</td><td>15,449,845</td></tr><tr><td rowspan="3">Senescence</td><td>D.Z65</td><td>d</td><td>2D</td><td>34,995,194</td><td>AX-89382235</td><td>T</td><td>0.52</td><td>10.62</td><td>-5.71</td><td>2.22</td><td>QTL.04</td><td>3</td><td>32,510,457</td><td>36,442,547</td></tr><tr><td>XPI</td><td>°Cd</td><td>4D</td><td>17,083,871</td><td>AX-89398511</td><td>A</td><td>0.67</td><td>5.77</td><td>-65.29</td><td>29.72</td><td>QTL.06</td><td>4</td><td>16,992,006</td><td>18,830,860</td></tr><tr><td>XPI</td><td>°Cd</td><td>2D</td><td>34,995,194</td><td>AX-89382235</td><td>T</td><td>0.52</td><td>6.11</td><td>50.98</td><td>28.82</td><td>QTL.04</td><td>1</td><td>33,532,314</td><td>36,458,074</td></tr><tr><td rowspan="3">VIT</td><td>XPI</td><td>°Cd</td><td>2A</td><td>45,486,914</td><td>AX-89761338</td><td>G</td><td>0.84</td><td>5.54</td><td>0.04</td><td>0.02</td><td>QTL.01</td><td>1</td><td>44,809,106</td><td>46,164,722</td></tr><tr><td>XPI</td><td>°Cd</td><td>2D</td><td>34,995,194</td><td>AX-89382235</td><td>T</td><td>0.52</td><td>6.24</td><td>56.48</td><td>26.04</td><td>QTL.04</td><td>2</td><td>33,217,357</td><td>36,453,954</td></tr><tr><td>XPI</td><td>°Cd</td><td>4D</td><td>17,083,871</td><td>AX-89398511</td><td>A</td><td>0.67</td><td>7.06</td><td>-83.51</td><td>34.52</td><td>QTL.06</td><td>16</td><td>16,826,225</td><td>23,189,173</td></tr><tr><td rowspan="3">START</td><td>XPI</td><td>°Cd</td><td>4D</td><td>17,083,871</td><td>AX-89398511</td><td>A</td><td>0.67</td><td>7.06</td><td>-83.51</td><td>34.52</td><td>QTL.06</td><td>16</td><td>16,826,225</td><td>23,189,173</td></tr><tr><td>XPI</td><td>°Cd</td><td>2B</td><td>25,953,023</td><td>AX-89528410</td><td>A</td><td>0.82</td><td>6.56</td><td>1.00</td><td>0.41</td><td>QTL.02</td><td>3</td><td>24,831,122</td><td>26,711,709</td></tr><tr><td>XPI</td><td>°Cd</td><td>2B</td><td>619,396,366</td><td>AX-89729195</td><td>T</td><td>0.91</td><td>6.47</td><td>1.82</td><td>0.44</td><td>QTL.03</td><td>170</td><td>605,382,213</td><td>652,111,584</td></tr><tr><td rowspan="3">Number of Spikes</td><td>NbGrain</td><td>g</td><td>4D</td><td>19,694,106</td><td>AX-89363486</td><td>G</td><td>0.64</td><td>6.88</td><td>7.44</td><td>3.28</td><td>QTL.06</td><td>21</td><td>16,840,726</td><td>23,189,154</td></tr><tr><td>NbGrain</td><td>g</td><td>4D</td><td>17,083,871</td><td>AX-89398511</td><td>A</td><td>0.67</td><td>9.37</td><td>0.38</td><td>0.14</td><td>QTL.06</td><td>32</td><td>16,664,691</td><td>25,447,819</td></tr><tr><td>NbGrain</td><td>g</td><td>7A</td><td>16,603,832</td><td>AX-89332971</td><td>G</td><td>0.68</td><td>5.87</td><td>3.60</td><td>1.84</td><td>QTL.09</td><td>3</td><td>15,907,587</td><td>17,300,104</td></tr><tr><td rowspan="3">Shoot</td><td>NbGrain</td><td>g</td><td>7A</td><td>16,603,832</td><td>AX-89332971</td><td>G</td><td>0.68</td><td>5.66</td><td>0.87</td><td>0.46</td><td>QTL.09</td><td>2</td><td>15,907,587</td><td>17,300,104</td></tr><tr><td>NbGrain</td><td>g</td><td>7A</td><td>16,603,832</td><td>AX-89332971</td><td>G</td><td>0.68</td><td>6.13</td><td>0.12</td><td>0.06</td><td>QTL.09</td><td>3</td><td>15,907,587</td><td>17,300,104</td></tr><tr><td>NbGrain</td><td>g</td><td>4D</td><td>19,459,578</td><td>AX-89445201</td><td>G</td><td>0.61</td><td>6.44</td><td>35.82</td><td>15.11</td><td>QTL.06</td><td>9</td><td>16,839,627</td><td>22,372,934</td></tr><tr><td rowspan="3">Tillers</td><td>NbGrain</td><td>g</td><td>4D</td><td>19,459,578</td><td>AX-89445201</td><td>G</td><td>0.61</td><td>6.44</td><td>35.82</td><td>15.11</td><td>QTL.06</td><td>9</td><td>16,839,627</td><td>22,372,934</td></tr><tr><td>NbGrain</td><td>g</td><td>4D</td><td>17,430,704</td><td>AX-89345211</td><td>A</td><td>0.64</td><td>8.99</td><td>1.69</td><td>0.61</td><td>QTL.06</td><td>28</td><td>13,163,753</td><td>24,863,780</td></tr><tr><td>NbGrain</td><td>g</td><td>7A</td><td>25,929,273</td><td>AX-89765935</td><td>C</td><td>0.63</td><td>5.86</td><td>3.13</td><td>1.71</td><td>QTL.10</td><td>21</td><td>25,228,927</td><td>26,941,996</td></tr></tbody></table>

*Trait abbreviations* D.Z65, flowering date; XPI and YPI, x-axis and y-axis of the inflexion point of senescence curve; VIT, rate of senescence, START, beginning of senescence; END, end of senescence; SPP, number of fertile spikes per plant; NbGrain, number of grains; GW, grain weight; TKW, thousand kernel weight; GrainArea, area of grain; GrainWidth, width of grain. Position is relative to the reference genome Chinese Spring refSeq1.0 (IWGSC et al. 2018)
areas were significantly associated with the SNP main effect for flowering date considering the estimated 5.26 threshold for -log (*P*value) (Fig. 5a, Table 2). No SNP×treatment interaction was significantly associated with flowering date (Fig. 5b). TKW from the main shoot was only associated with an SNP main effect on chromosome 7A and to an SNP×treatment interaction effect on chromosome 4B (Fig. 5c, d, Tables 2 and 3).

Allelic effects for the most strongly associated SNP for flowering date (chromosome 2D) and the most strongly associated SNP for TKW (4B) are plotted in Fig. 6. First, we can see that there was no genotype dispersion around the bisecting line for flowering date (Fig. 6a), because there was no genotype × treatment interaction, whereas a large dispersion is observed for TKW (Fig. 6b). Second, the figure shows the difference between SNP alleles. For the SNP significantly associated with the main effect (Fig. 6a), genotypes are separated only along the bisecting line (y=x). However, for the SNP significantly associated with treatment interaction (Fig. 6b), genotypes are separated following the second bisecting line (y=-x).

Of the four QTLs associated with flowering date (Table 2), the strongest one is located on chromosome 2D colocalized with the photoperiod sensitivity gene *Ppd-D1* (Beales et al. 2007). This QTL also colocalized with the start of senescence and XPI traits. The effect of the minor allele had opposite effects, meaning that earlier flowering genotypes had a delayed start of senescence. Two more QTLs were associated with senescence traits, one on chromosome 2A was associated with the maximum senescence rate (VIT), and one on chromosome 4D was associated with the start and end of senescence and colocalized with the *Rht-D1* dwarfing gene (Peng et al. 1999). The latter QTL was also associated with the number of grains per spike in the main shoot and tillers, but with an opposite effect to those of the associated senescence traits in that early senescing plants had more grains and heavier spikes than later senescing ones. Lastly, one QTL localized on chromosome 7A was associated with TKW for both the main shoot and the tillers.

Seven QTLs were significantly associated with the SNP×treatment interaction effect (Table 3). Four were related to senescence traits, one located on chromosome 4B was associated with main shoot TKW and related grain traits (area, length, and width), and two located on chromosomes 3B and 6B, respectively, were associated with grain weight and TKW traits for tillers. All the effects of the major allele of the most significant SNP per QTL are presented in Table 3. For example, varieties with allele G at QTL 4B, which is the major allele with a frequency at 0.64, performed less well in S treatment than varieties with allele A. This effect depicted in Fig. 6a amounts to -3.51 (± 2.15) g. None of the QTLs associated with an SNP by treatment interaction showed a significant SNP main effect.

## Discussion

This study explored the impact of post-anthesis heat stress on a European elite bread wheat panel in controlled conditions. Physiological traits such as flowering time and grain morphology were recorded and senescence was modeled. We performed a GWAS to dissect the genetic determinants of heat tolerance. We identified 10 QTLs associated with at least one trait as a main effect and seven QTLs associated with an interaction with post-anthesis heat stress.

## Experiment validation and stress impact

Before interpreting GWAS results, it is important to correctly describe the stress undergone by the plants. Post-anthesis temperature was increased daily for 10 days by on average + 8.6 °C in the day and + 5.8 °C at night. This heat stress reduced grain weight by 24% for the main shoot compared with the optimum temperature of 21 °C in the day and 17 °C at night. This result is consistent with other reports of grain weight loss under post-anthesis heat stress (Stone and Nicolas 1994; Gibson and Paulsen 1999; Spieritz et al. 2006). The TKW decrease was correlated with a reduction in grain width (SI= 4.5%) rather than in grain length (SI= 10.0%, Online Resource 6). The higher sensitivity of grain width over length has already been reported (Jamil et al. 2017).

We decided to apply the S condition three days after the main shoot flowered, so as not to impact the number of grains on the main shoot. This was indeed the case as for both years the number of grains per main shoot was not significantly different between the two temperature treatments (Table 1). The impact of heat stress was therefore focused on only one of the yield components, TKW. However, it has been shown that high pre-flowering temperatures may cause flower sterility (e.g., Farooq et al. 2011; Barber et al. 2017). Yield components develop one after the other and there are multiple compensation mechanisms at work that are liable to bias the interpretation (e.g., Slafer and Andrade 1993; Gate 1995). In our case, fewer grain per spike could have caused an increase in thousand kernel weight. Furthermore, when plants survive a long period of stress, they may set in motion various acclimation mechanisms with different genetic determinisms (Gaspar et al. 2002). Our experiment was designed solely to identify heat stress QTL specific to the early phase of the main shoot grain filling. For the same reason, we chose to analyze grains from the main shoot and from tillers separately. Heat stress may have a different impact on the main shoot and tillers, due to the naturally sequential flowering dates for tillers, up to nine days later than the main shoot (Jones

<figure><img src="image_15.png" /></figure>
<table><caption>**Table 3** List of associated SNP according to a GWAS of 199 wheat varieties grown for two years with stress or no-stress post-flowering temperature treatments: SNP with a significant interaction effect</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Trait</th><th rowspan="2">Units</th><th rowspan="2">Chr</th><th rowspan="2">Position (bp)</th><th rowspan="2">SNP</th><th rowspan="2">Major Allele</th><th rowspan="2">Freq</th><th colspan="2">SNP by treatment interaction</th><th rowspan="2">QTL</th><th colspan="2">QTL boundary</th></tr><tr><th>LOD</th><th>Effect of the major allele effect in S</th><th>Effect SD</th><th>Nb. SNP</th><th>Min. (bp)</th><th>Max. (bp)</th></tr></thead><tbody><tr><td rowspan="4">Senescence</td><td>XPI</td><td>°Cd</td><td>7B</td><td>693,020,724</td><td>AX-89578121</td><td>C</td><td>0.91</td><td>5.71</td><td>60.40</td><td>34.72</td><td>QTL.17</td><td>2</td><td>692,252,452</td><td>693,926,665</td></tr><tr><td>XPI</td><td>°Cd</td><td>1D</td><td>17,642,951</td><td>AX-89766697</td><td>A</td><td>0.62</td><td>5.38</td><td>-32.45</td><td>31.54</td><td>QTL.12</td><td>1</td><td>13,902,574</td><td>21,383,328</td></tr><tr><td>YPI</td><td>µg cm<sup>-3</sup></td><td>1B</td><td>631,711,683</td><td>AX-89638269</td><td>G</td><td>0.81</td><td>6.69</td><td>1.24</td><td>0.82</td><td>QTL.11</td><td>2</td><td>630,767,028</td><td>632,656,378</td></tr><tr><td>AUCK</td><td></td><td>4A</td><td>328,128,331</td><td>AX-89642074</td><td>A</td><td>0.91</td><td>5.55</td><td>-3119.83</td><td>2355.64</td><td>QTL.14</td><td>36</td><td>38,974,540</td><td>472,729,606</td></tr><tr><td rowspan="5">Main Shoot</td><td>TKW</td><td>g</td><td>4B</td><td>1,841,022</td><td>AX-89364414</td><td>G</td><td>0.64</td><td>5.99</td><td>-3.51</td><td>2.15</td><td>QTL.15</td><td>11</td><td>1,339,055</td><td>2,852,831</td></tr><tr><td>GrainArea</td><td>mm<sup>2</sup></td><td>4B</td><td>1,853,262</td><td>AX-89500218</td><td>A</td><td>0.63</td><td>6.42</td><td>-0.78</td><td>0.53</td><td>QTL.15</td><td>34</td><td>1,384,350</td><td>3,700,062</td></tr><tr><td>GrainWidth</td><td>mm</td><td>4B</td><td>1,853,262</td><td>AX-89500218</td><td>A</td><td>0.63</td><td>5.83</td><td>-0.11</td><td>0.07</td><td>QTL.15</td><td>9</td><td>1,310,533</td><td>2,334,232</td></tr><tr><td>Grain.Length</td><td>mm</td><td>4B</td><td>1,854,364</td><td>AX-89738651</td><td>T</td><td>0.64</td><td>5.49</td><td>-1.14E-01</td><td>1.10E-01</td><td>QTL.15</td><td>25</td><td>1,438,558</td><td>3,697,489</td></tr><tr><td>Tillers</td><td>GW</td><td>g</td><td>3B</td><td>150,697,370</td><td>AX-89468413</td><td>C</td><td>0.91</td><td>5.40</td><td>1.48</td><td>0.74</td><td>QTL.13</td><td>1</td><td>141,244,647</td><td>160,150,093</td></tr><tr><td rowspan="4"></td><td>TKW</td><td>g</td><td>6B</td><td>695,794,696</td><td>AX-89726711</td><td>A</td><td>0.58</td><td>5.76</td><td>2.53</td><td>1.96</td><td>QTL.16</td><td>7</td><td>694,970,068</td><td>696,420,120</td></tr><tr><td>GrainArea</td><td>mm<sup>2</sup></td><td>6B</td><td>695,794,696</td><td>AX-89726711</td><td>A</td><td>0.58</td><td>5.85</td><td>0.58</td><td>0.51</td><td>QTL.16</td><td>8</td><td>694,970,068</td><td>696,420,120</td></tr><tr><td>Grain.Length</td><td>mm</td><td>6B</td><td>695,794,696</td><td>AX-89726711</td><td>A</td><td>0.58</td><td>6.13</td><td>1.01E-01</td><td>1.10E-01</td><td>QTL.16</td><td>9</td><td>694,970,068</td><td>696,420,120</td></tr></tbody></table>

*Trait abbreviations* XPI and YPI, respectively, the x-axis and y-axis of the inflection point of the senescence curve; AUCK, total surface under the senescence curve; GW, grain weight; TKW, thousand kernel weight; GrainArea, area of grain; GrainWidth, width of grain; Grain.Length, length of grain. Position is relative to reference genome Chinese Spring refSeq1.0 (IWGSC et al. 2018). S, heat stress treatment
<figure><img src="image_16.png" /><figcaption>**Fig. 6** Correlation between gBLUP values for 199 wheat varieties grown in two post-flowering temperature treatments. **a** Correlation between flowering dates for stressed and non-stressed plants showing the minor allele (T) effect of QTL 2D (blue) and the effect of the other allele (red). **b** Correlation between TKW for stressed and</figcaption></figure>

<figure><img src="image_17.png" /></figure>

non-stressed plants showing the minor allele (G) effect of QTL 4B (blue) and the effect of the other allele (red). Dotted lines correspond to average values and solid lines to linear regressions for the minor allele (blue) and the other allele (red)

et al. 2017). For example, heat stress may impact the num-
ber of grains on tillers by reducing pollen fertility (Farooq
et al. 2011), while the fertilization has already occurred on
the main shoot. In both our experiments (2016 and 2017),
the heat stress decreased the number of grains on tillers
by more than 5%, while decreasing grain weight for both
the main shoot and tillers, with SI at 23.9% and 17.5%,
respectively. The smaller decrease for tillers could be due
either to compensation for having fewer grains and/or the
shorter period of stress experienced due to later flowering.

For monocarpic species exposed to terminal heat stress,
it is useful to determine whether the stress damages the pho-
tosynthetic source more than the reproductive sink. In this
study, we monitored leaf chlorophyll content during grain
filling and observed only a weak impact of stress on senes-
cence traits. The significant treatment effects for senescence
parameters were for the number of °Cd up to the start and
inflexion point of senescence, which had SI of 5.26% and
3.75%, respectively. These observations confirm the view
that high temperatures, as applied in our conditions, shorten
the plant cycle and limit the amount of light intercepted
(Sofield et al. 1977). Clearly, heat stress accelerates the rate
of grain filling while shortening its duration (Sofield et al.
1977; Dias and Lidon 2009). It has been estimated that for
every °C above the optimal growing temperature, the dura-
tion of grain filling is reduced by 2–3 days (Streck 2005).
Such an impact was visualized in our results by plotting the
decrease of chlorophyll content as a function of Julian days
for both treatments (Online Resource 7), but when senes-
cence dynamics was plotted as a function of degree days the
difference was less marked (Fig. 3, Tables 1).

We applied heat stress just after anthesis so it is likely
that the number of endosperm cells was altered with conse-
quences on grain filling and final grain size. Indeed, wheat
grain usually develops in three phases: the lag phase, filling
and maturation (Gate 1995; Acevedo et al. 2002). During the

lag phase, which lasts about 15 days or 250°Cd, the number
of endosperm cells is fixed. While the supply of assimilates
available to the grain during the first two weeks after anthe-
sis is thought to determine the cell number (Brocklehurst
1977), once fixed, it is the number of cells that regulates the
rate of accumulation of dry matter during the filling phase.

This work was carried out on a large collection of elite
European winter varieties that were not specifically bred
for heat stress tolerance. A large variability was observed
among the grain and senescence traits measured. Specifi-
cally, a large G×T variance was observed for grain weight,
meaning that genotypes responded very differently to heat
stress. This information could already be useful for breeders
in assessing the performance of the varieties and selecting
new genitors for breeding crosses. For example, varieties
such as Artdeco and Saint-Ex only lost about 10% of their
TKW while others lost more than 35% (Online Resource 4).
The large variability observed was associated with moderate
to high heritability (Table 2) meaning that this experiment
was suitable for detecting QTL linked with these traits.

QTL associated with measured traits

We detected 10 QTLs associated with at least one trait main effect (Table 2). Three QTLs were associated with senescence traits. QTL 4D is co-located with the *Rht-D1* dwarfing gene. *Rht* genes have many effects on wheat development (e.g., Borrell et al. 1991; Peng et al. 1999). The G allele of molecular marker AX-89398511, the most frequently associated SNP within this QTL, was positively associated with plant height (Touzy et al. 2019), meaning that genotypes with the G allele carry the wild-type form of *Rht-D1*. In the present study, this allele is associated to a higher number of grains for the main shoot and tillers and with earlier senescence (earlier point of inflexion and end of senescence). The effect of *Rht* genes on the number of grains per spike
has often been reported. To investigate the possible causes, Miralles et al. (1988) compared isogenic lines for *Rht-D1* and *Rht-B1*, and found that the semi-dwarf alleles increased the number of fertile florets per spikelet, possibly because assimilate partitioning to the spike was more favorable than to the stem during the pre-anthesis period. Association between *Rht* genes and senescence has also already been reported (Blake et al. 2009; Camargo et al. 2016). Christopher et al. (2018) reported that earlier senescence (defined in that case as thermal time from anthesis to 10% and 50% loss of maximal greenness) was associated with the *Rht* semi-dwarfing allele in the SeriM82/Hartog population. They could not conclude whether this gene was directly associated with senescence or with linked genes. More recently, Jobson et al. (2019) compared the photosynthesis characteristics of *Rht-B1* isogenic lines in space-planted, irrigated field conditions. They concluded that the semi-dwarf allele reduced flag leaf chlorophyll content at anthesis and photosynthetic rate per unit area at anthesis and 14 days after anthesis. Despite these changes, they did not detect any significant modifications in genes global expression due to the presence of the semi-dwarfing allele that could help explain these differences. In addition, in our conditions, plants were placed at random in the greenhouse, so we cannot ignore the possibility that semi-dwarf plants were in competition for light with taller ones, which could have impacted the number of grains and senescence.

The photoperiod sensitivity gene *Ppd-D1* on chromosome 2D co-localizes with a QTL for flowering date and a QTL for the onset of senescence and the point of inflexion of the logistic curve. Co-localizations between phenology and senescence traits are widely reported (e.g., Verma et al. 2004; Bogard et al. 2011; Pinto et al. 2016; Camargo et al. 2016; Christopher et al. 2018). Verma et al. (2004) identified a QTL located on the short arm of chromosome 2D which was associated with higher yield and delayed senescence, but not with flowering date in the conditions tested. By characterizing a Beaver/Soissons bi-parental population those authors showed that the allele responsible for delayed senescence came from the Soissons cultivar. Soissons is in our panel, so we checked the polymorphism of the most strongly associated molecular marker and found a similar result. This means that earlier varieties, with the mutant version of *Ppd-D1* which is insensitive to photoperiod, have flag leaves that stay green for longer. Based on several studies, Bogard et al. (2011) presented different hypotheses for the relationship between wheat precocity and senescence. Three hypotheses involved the direct or indirect effect of comparing varieties differing in precocity in the same field trial. Essentially, in field experiments, terminal stresses, such as high temperature and water and nitrogen deficiencies, often occur at about the same date for all the genotypes. As senescence is highly dependent on these environmental conditions, late flowering

plants would be expected to have shorter stay-green phenotypes. However, as our experiments were carried out in controlled conditions where water and nitrogen were non-limiting, those assumptions do not fit. Another hypothesis, suggested by Wingler et al. (2010) for *A. thaliana*, is that the relationship between precocity and leaf senescence is linked to a change in metabolism soon after anthesis, when the dismantling of structural components of leaves releases sugars after anthesis, which could trigger senescence.

Only one QTL on chromosome 7A was associated with TKW, but it did not co-localize with the QTL identified in our previous work (Touzy et al. 2019). However, the two studies were conducted in very different conditions. In the present study, having only one plant per pot (12 cm wide) resulted in a plant density of about 70 plants m⁻² which is very low compared to typical autumn-sown fields (~200–300 plants m⁻²). The study of a complex trait such as heat tolerance is easier in controlled conditions and we assumed that sowing density and potential differences in tillering and TKW would not strongly affect the assessment of genotype behavior. For agronomical traits, such as yield components, controlled-condition studies for a breeding or agricultural context are not as representative as field experiments. As stressed by Tardieu (2012), both kinds of studies are relevant as they enable different behaviors to be accessed.

## QTL associated with the response to heat stress

Seven QTLs associated with the response to heat treatment were identified. None of them are associated with flowering time, which validated our protocol in allowing us to identify QTL linked to heat stress tolerance, without confounding effects due to escape.

Four QTLs are associated with senescence traits. On chromosome 4A, the QTL associated with the area under the curve (AUCK) is centromeric. As linkage disequilibrium is high in this region, the probability of identifying co-localization with other QTLs or genes from literature is therefore too high. For the remaining QTLs, it is relevant to check co-localization with our previous field trials of the same panel of varieties (Touzy et al. 2019). Trials from a large multi-environmental network were grouped according to their water stress scenarios: optimal, no stress, terminal stress, stress during anthesis, and stress from booting to harvest. Even though drought stress is often associated with heat stress, the heat stress component of the environments was not characterized. Three QTLs associated with the response of grain weight to heat stress were identified in the present study. Two of them co-localize with others found in the multi-environmental network study (Touzy et al. 2019). The SNP most associated with the chromosome 3B QTL was associated with grain yield in environments with water stress
during the anthesis period. The chromosome 4B QTL associated with the response of main shoot TKW co-localized with a QTL associated with TKW in all the water stress scenarios and in optimal conditions. More than 14% of the variability in plant stress responses could be explained by the QTL we have detected here. We have checked reports that place QTLs for TKW and grain yield in hot conditions on chromosome 4B, many of which are reviewed by Zhang et al. (2010), Bouffier (2014), Tricker et al. (2018) and Guan et al. (2020). Pinto et al. (2010) also reported two QTLs on chromosome 4B for TKW and plant height in hot conditions. Noting that several studies had identified the latter QTL the authors suggested that this region may co-localize with the *Rht-B1* gene. However, the bi-parental Seri M82/ Babax population used to identify this QTL did not segregate for any known *Rht* genes. We checked the positions of our QTL and *Rht-B1* (Online Resource 8) and found that both the physical (~29 Mb) and genetic (~40 cM) distances between them are large. Therefore, it seems likely that one or more genes linked to thermal stress tolerance are present in this area, independently of *Rht-B1*.

Even though these QTLs did not explain much of the variability in the stress response, co-localizing QTL between controlled condition and field experiments allowed us to validate our approach.

## Conclusion

We have investigated terminal heat stress tolerance using a European panel of elite winter wheat varieties. An optimal condition and a post-anthesis heat treatment applied to plants grown in the greenhouse showed how widely varieties differ in their ability to tolerate this stress. Grain weight was strongly impacted by this stress, while senescence traits were only slightly influenced. The genetic determinant of heat tolerance was determined by a GWAS. We identified a significant SNP × treatment interaction for TKW on the telomeric region of the short arm of chromosome 4B. Focusing on a well-defined terminal heat stress, we integrate our findings with other stress scenarios to help identify the genomic regions needed to develop heat stress-tolerant crops.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00122-021-04008-5.

**Acknowledgements** Part of the data was obtained thanks to the support of the PIA (Investment for the Future Program) BreedWheat (ANR-10-BTBR-03) project funded by the National Research Agency (ANR), FranceAgriMer, and the French Plant Breeding Support Funds (FSOV-2012D). The authors are also grateful to the ANRT (Association Nationale de la Recherche et de la Technologie) and ARVALIS Institut du végétal for supporting the Ph.D. thesis (2015/0686). The authors would like to thank Pierre Dubreuil (Biogemma), Gilles Charmet and Christine Girousse (INRAE, Université Clermont Auvergne, UMR

GDEC) for their aid in methodology and the interpretation of results and the Biogemma greenhouse team Emmanuel Brignat, Isabelle Delaune, Guillaume Gondal, Valerie Legrand, Sylvain Levadoux. Thanks to the Accent Europe company (http://www.accentEurope.fr) and to Rachel Carol (Bioscience Editing, France) for editing the English in the manuscript.

**Author's contributions statement** The statistical analyses were conducted by GT during his Ph.D. thesis supervised by SP and JLG. GT, JLG, SL and ER wrote the manuscript. VL and XD carried out the greenhouse experiments. All authors contributed to methods and interpretations, and in reviewing the manuscript.

**Funding** Agence nationale de la recherche, ANR-10-BTBR-03, Association nationale de la recherche et de la technologie, 2015/0686, Gaëtan TOUZY, Fonds de soutien à l'obtention végétale, 2012D.

## Declarations

**Conflicts of interests** The authors have no conflicts of interest to declare that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Acevedo E, Silva P, Silva H (2002) Growth and wheat physiology. In: Curtis B, Rajaram S, Gomez-Macpherson H (eds) Bread Wheat, Food and Agriculture Organization of the United Nations, Rome

Barber HM, Lukac M, Simmonds J, Semenov MA, Gooding MJ (2017) Temporally and genetically discrete periods of wheat sensitivity to high temperature. Front Plant Sci. https://doi.org/10.3389/fpls.2017.00051

Barnabás B, Jäger K, Fehér A (2007) The effect of drought and heat stress on reproductive processes in cereals. Plant Cell Environ 31:11–38

Bates DM, Chambers JM (1992) Nonlinear models. In: Chambers SJM, Hastie TJ (eds) Chapter 10 of Statistical Models in S. Wadsworth & Brooks/Cole

Beales J, Turner A, Griffiths S, Snape J, Laurie D (2007) A Pseudo-Response Regulator is misexpressed in the photoperiod insensitive Ppd-D1a mutant of wheat (*Triticum aestivum* L.). Theor Appl Genet 115:721–733

Blake NK, Lanning SP, Martin JM, Doyle M, Sherman JD, Naruoka Y, Talbert LE (2009) Effect of variation for major growth habit genes on maturity and yield in five spring wheat populations. Crop Sci 49:1211–1220

Bogard M, Jourdan M, Allard V, Martre P, Perretant MR, Ravel C, Heumez E, Orford S, Snape J, Griffiths S, Gaju O, Foulkes J, Le Gouis J (2011) Anthesis date mainly explained correlations
between post-anthesis leaf senescence, grain yield, and grain protein concentration in a winter wheat population segregating for flowering time QTLs. J Exp Bot 62:3621–3636

Borrell AK, Incoll LD, Dalling MJ (1991) The influence of the Rht1 and Rht2 alleles on the growth of wheat stems and ears. Ann Bot 67:103–110

Borrell A, Hammer G, Douglas A (2000) Does maintaining green leaf area in sorghum improve yield under drought? I. Leaf Grow Senesce Crop Sci 40:1016–1037

Bouffier B, Derory J, Murigneux A, Reynolds M, Le Gouis J (2015) Clustering of environmental parameters discriminates drought and heat stress bread wheat (Triticum aestivum L.) trials. Agron J 107:1489–1503

Bouffier B (2014) Genetic and ecophysiological dissection of tolerance to drought and heat stress in bread wheat: from environmental characterization to QTL detection. Université Blaise Pascal-Clermont-Ferrand II

Breseghello F, Sorrells ME (2006) Association mapping of kernel size and milling quality in wheat (Triticum aestivum L.) cultivars. Genetics 172:1165–1177

Brisson N, Gate P, Gouache D, Charmet G, Oury F-X, Huard F (2010) Why are wheat yields stagnating in Europe? A comprehensive data analysis for France. Field Crop Res 119:201–212

Brocklehurst PA (1977) Factors controlling grain weight in wheat. Nature 266:348–349

Browning BL, Browning SR (2016) Genotype imputation with millions of reference samples. Am J Hum Genet 98:116–126

Butler D, Cullis BR, Gilmour AR, Gogel BJ (2009) ASReml-R reference manual. Queensland Department of Primary Industries, Australia

Calderini DF, Slafer GA (1998) Changes in yield and yield stability in heat during the 20th century. Field Crop Res 57:335–347

Calderini DF, Abeledo LG, Savin R, Slafer GA (1999) Effect of temperature and carpel size during pre-anthesis on potential grain weight in wheat. J Agric Sci 132:453–459

Camargo AV, Mott R, Gardner KA, Mackay IJ, Corke F, Doonan JH, Kim JT, Bentley AR (2016) Determining phenological patterns associated with the onset of senescence in a wheat MAGIC mapping population. Front Plant Sci. https://doi.org/10.3389/fpls.2016.01540

Christopher JT, Manschadi AM, Hammer GL, Borrell AK (2008) Developmental and physiological traits associated with high yield and stay-green phenotype in wheat. Aust J Agric Res 59:354–364

Christopher JT, Christopher MJ, Borrell AK, Fletcher S, Chenu K (2016) Stay-green traits to improve wheat adaptation in well-watered and water-limited environments. J Exp Bot 67:5159–5172

Christopher M, Chenu K, Jennings R, Fletcher S, Butler D, Borrell A, Christopher J (2018) QTL for stay-green traits in wheat in well-watered and water-limited environments. Field Crop Res 217:32–44

Commuri PD, Jones RJ (2001) High temperatures during endosperm cell division in maize: a genotypic comparison under in vitro and field conditions. Crop Sci 41:1122–1130

Cormier F, Le Gouis J, Dubreuil P, Lafarge S, Praud S (2014) A genome-wide identification of chromosomal regions determining nitrogen use efficiency components in wheat (Triticum aestivum L.). Theor Appl Genet 127:2679–2693

Cossani CM, Reynolds MP (2012) Physiological traits for improving heat tolerance in wheat. Plant Physiol 160:1710–1718

Dias AS, Lidon FC (2009) Evaluation of grain filling rate and duration in bread and durum wheat, under heat stress after anthesis. J Agron Crop Sci 195:137–147

Efron B (1978) Regression and ANOVA with zero-one data: measures of residual variation. J Am Stat Assoc 73:113–121

Farooq M, Bramley H, Palta JA, Siddique KHM (2011) Heat stress in wheat during reproductive and grain-filling phases. Crit Rev Plant Sci 30:491–507

Farooq M, Hussain M, Siddique KHM (2014) Drought stress in wheat during flowering and grain-filling periods. Crit Rev Plant Sci 33:331–349

Gao X, Becker LC, Becker DM, Starmer JD, Province MA (2009) Avoiding the high Bonferroni penalty in genome-wide association studies. Genet Epidemiol 34:100–105

Gao F, Wen W, Liu J, Rasheed A, Yin G, Xia X, Wu X, He Z (2015) Genome-wide linkage mapping of QTL for yield components, plant height and yield-related physiological traits in the Chinese wheat cross Zhou 8425B/Chinese Spring. Front Plant Sci. https://doi.org/10.3389/fpls.2015.01099

Gaspar T, Franck T, Bisbis B, Kevers C, Laurent J, Hausman JF, Dommes J (2002) Concepts in plant stress physiology. Application to plant tissue cultures. Plant Growth Regul 37:263–285

Gate P (2007) Evolution du climat et production du blé. Perspectives Agricoles 336:24–31

Gate P, Blondlot A, Gouache D, Deudon O, Vignier L (2008) Impacts du changement climatique sur la croissance et le développement du blé en France. Quelles solutions et quelles actions à développer ? Corps Gras Lipides 15:332–336

Gate P, Brisson N, Gouache D (2010) Les causes du plafonnement du rendement du blé en France : d'abord une origine climatique. Comptes-Rendus De L'académie D'agriculture De France 96:1–17

Gate P (1995) Ecophysiologie du blé : de la plante à la culture. Lavoisier, France

Gibson LR, Paulsen GM (1999) Yield components of wheat grown under high temperature stress during reproductive growth. Crop Sci 39:1841–1846

Guan PF, Shen XY, Mu Q, Wang YF, Wang XB, Chen YM, Zhao Y, Chen XY, Zhao AJ, Mao WW, Guo YW, Xin MM, Hu ZR, Yao YY, Ni ZF, Sun QX, Peng HR (2020) Dissection and validation of a QTL cluster linked to Rht-B1 locus controlling grain weight in common wheat (Triticum aestivum L.) using near-isogenic lines. Theor Appl Genet 133:2639–2653

Hunt LA, van der Poorten G, Pararajasingham S (1991) Postanthesis temperature effects on duration and rate of grain filling in some winter and spring wheats. Can J Plant Sci 71:609–617

IWGSC (2018) Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science 361:661+

Jamil M, Ali A, Ghafoor A, Akbar KF, Aziz A, Naveed NH, Yasin NA, Gul A, Mujeeb-Kazi A (2017) Digital image analysis of seed shape influenced by heat stress in diverse bread wheat germplasm. Pak J Bot 49:1279–1284

Jobson EM, Johnston RE, Oiestad AJ, Martin JM, Giroux MJ (2019) The impact of the wheat Rht-B1b semi-dwarfing allele on photosynthesis and seed development under field conditions. Front Plant Sci. https://doi.org/10.3389/fpls.2019.00051

Johnson IR, Thornley JHM (1985) Temperature dependence of plant and crop processes. Ann Bot 55:1–24

Jones HE, Lukac M, Brak B, Martinez-Eixarch M, Alhomedhi A, Gooding MJ, Wingen LU, Griffiths S (2017) Photoperiod sensitivity affects flowering duration in wheat. J Agric Sci 155:32–43

Kendall MG, Stuart A (1979) The advanced theory of statistics, 4th edn. Griffin, London

Lim PO, Kim HJ, Gil Nam H (2007) Leaf senescence. Annu Rev Plant Biol 58:115–136

Lobell DB, Schlenker W, Costa-Roberts J (2011) Climate trends and global crop production since 1980. Science 333:616–620

Van Der Loo MPJ (2010) Distribution based outlier detection for univariate data. Statistics Netherland, The Hague/Heerlen

Lopes MS, Reynolds MP (2012) Stay-green in spring wheat can be determined by spectral reflectance measurements (normalized
difference vegetation index) independently from phenology. J Exp Bot 63:3789–3798

Miralles D, Katz S, Colloca A, Slafer G (1988) Floret development in near isogenic wheat lines differing in plant height. Field Crop Res 59:21–30

Nicolas M, Gleadow R, Dalling M (1984) Effects of drought and high temperature on grain growth in wheat. Funct Plant Biol 11:553

Olivares-Villegas JJ, Reynolds MP, McDonald GK (2007) Drought-adaptive attributes in the Seri/Babax hexaploid wheat population. Funct Plant Biol 34:189–203

Pebesma EJ (2004) Multivariable geostatistics in S: the gstat package. Comput Geosci 7:683–691

Peltonen-Sainio P, Jauhiainen L, Laurila IP (2009) Cereal yield trends in northern European conditions: changes in yield potential and its realisation. Field Crop Res 110:85–90

Peng J, Richards DE, Hartley NM, Murphy GP, Devos KM, Flintham JE, Beales J, Fish LJ, Worland AJ, Pelica F, Sudhakar D, Christou P, Snape JW, Gale MD, Harberd NP (1999) ‘Green revolution’ genes encode mutant gibberellin response modulators. Nature 400:256–261

Pinto RS, Reynolds MP, Mathews KL, McIntyre CL, Olivares-Villegas J-J, Chapman SC (2010) Heat and drought adaptive QTL in a wheat population designed to minimize confounding agronomic effects. Theor Appl Genet 121:1001–1021

Pinto RS, Lopes MS, Collins NC, Reynolds MP (2016) Modelling and genetic dissection of staygreen under heat stress. Theor Appl Genet 129:2055–2074

Porter JR, Gawith M (1999) Temperatures and the growth and development of wheat: a review. Eur J Agron 10:23–36

R Development Core Team (2011) A Language and Environment for Statistical Computing. Vienna, Austria: the R Foundation for Statistical Computing. ISBN: 3-900051-07-0. Available online at http://www.R-project.org/

Rife T, Graybosch R, Poland J (2019) A field-based analysis of genetic improvement for grain yield in winter wheat cultivars developed in the US Central Plains from 1992 to 2014. Crop Sci 59:905–910

Rimbert H, Darrier B, Navarro J, Kitt J, Choulet F, Leveugle M, Duarte J, Rivière N, Eversole K, Consortium TIWGS, Le Gouis J, Consortium TB, Davassi A, Balfourier F, Le Paslier M-C, Berard A, Brunel D, Feuillet C, Poncet C, Sourdille P, Paux E (2018) High throughput SNP discovery and genotyping in hexaploid wheat. PLOS One 13:e0186329

Rincent R, Moreau L, Monod H, Kuhn E, Melchinger AE, Malvar RA, Moreno-Gonzalez J, Nicolas S, Madur D, Combes V, Dumas F, Altmann T, Brunel D, Ouzunova M, Flament P, Dubreuil P, Charcosset A, Mary-Huard T (2014) Recovering power in association mapping panels with variable levels of linkage disequilibrium. Genetics 197:375–387

Semenov MA, Shewry PR (2011) Modelling predicts that heat stress, not drought, will increase vulnerability of wheat in Europe. Sci Rep 1:66

Shirdelmoghanloo H, Taylor JD, Lohraseb I, Rabie H, Brien C, Timmins A, Martin P, Mather DE, Emebiri L, Collins NC (2016) A QTL on the short arm of wheat (*Triticum aestivum* L.) chromosome 3B affects the stability of grain weight in plants exposed to a brief heat shock early in grain filling. BMC Plant Biology 16

Slafer GA, Andrade FH (1993) Physiological attributes related to the generation of grain yield in bread wheat cultivars released at different eras. Field Crop Res 31:351–367

Sofield I, Evans L, Cook M, Wardlaw I (1977) Factors influencing the rate and duration of grain filling in wheat. Aust J Plant Physiol 4:785–797

Spiertz HJH, Hamer RJ, Xu H, Primo-Martin C, Don C, van der Putten PEL (2006) Heat stress in wheat (*Triticum aestivum* L.): Effects on grain growth and quality traits. Eur J Agron 25:89–95

Stocker TF, Qin D., Plattner, G.-K., et al. (2013) Climate Change 2013: The Physical Science Basis. In: IPCC CUP, Cambridge, United Kingdom and New York, NY, USA (ed)

Stone P, Nicolas M (1994) Wheat cultivars vary widely in their responses of grain yield and quality to short periods of post-anthesis heat stress. Aust J Plant Physiol 21:887–900

Streck NA (2005) Climate change and agroecosystems: the effect of elevated atmospheric CO2 and temperature on crop growth, development, and yield. Ciênda Rural 35:730–740

Taiz L, Zeiger E (2002) Plant physiology, 3rd edn. W. H. Freeman, New York

Tardieu F (2012) Any trait or trait-related allele can confer drought tolerance: just design the right drought scenario. J Exp Bot 63:25–31

Tashiro T, Wardlaw IF (1990) The effect of high temperature at different stages of ripening on grain set, grain weight and grain dimensions in the semi-dwarf wheat “Banks.” Ann Bot 65:51–61

Telfer P, Edwards J, Norman A, Bennett D, Smith A, Able JA, Kuchel H (2021) Genetic analysis of wheat (*Triticum aestivum*) adaptation to heat stress. Theor Appl Genet 134:1387–1407

Thomas H, Howarth CJ (2000) Five ways to stay green. J Exp Bot 51:329–337

Touzy G, Rincent R, Bogard M, Lafarge S, Dubreuil P, Mini A, Deswarte J-C, Beauchêne K, Le Gouis J, Praud S (2019) Using environmental clustering to identify specific drought tolerance QTLs in bread wheat (T. aestivum L.). Theor Appl Genet 132:2859–2880

Tricker PJ, ElHabti A, Schmidt J, Fleury D (2018) The physiological and genetic basis of combined drought and heat tolerance in wheat. J Exp Bot 69:3195–3210

VanRaden PM (2008) Efficient methods to compute genomic predictions. J Dairy Sci 91:4414–4423

Verma V, Foulkes MJ, Worland AJ, Sylvester-Bradley R, Caligari PDS, Snape JW (2004) Mapping quantitative trait loci for flag leaf senescence as a yield determinant in winter wheat under optimal and drought-stressed environments. Euphytica 135:255–263

Wheeler R, Hong TD, Ellis RH, Batts RR, Morison JIL, Hadley P (1996) The duration and rate of grain growth, and harvest index, of wheat (*Triticum aestivum* L.) in response to temperature and CO₂. J Exp Bot 47:623–630

Wingler A, Purdy SJ, Edwards S-A, Chardon F, Masclaux-Daubresse C (2010) QTL analysis for sugar-regulated leaf senescence supports flowering-dependent and -independent senescence pathways. New Phytol 185:420–433

Xu Q, Paulsen AQ, Guikema JA, Paulsen GM (1995) Functional and ultrastructural injury to photosynthesis in wheat by high temperature during maturation. Environ Exp Bot 35:43–54

Zadoks JC, Chang TT, Konzak CF (1974) A decimal code for the growth stages of cereals. Weed Res 14:415–421

Zhang L-Y, Liu D-C, Guo X-L, Yang W-L, Sun J-Z, Wang D-W, Zhang A (2010) Genomic distribution of quantitative trait loci for yield and yield-related traits in common wheat: Meta-QTL analysis for its related traits in wheat. J Integr Plant Biol 52:996–1007

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.