# Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions

Justin Blancon, Clément Buet, Pierre Dubreuil, Marie-Hélène Tixier, Frédéric Baret, Sébastien Praud

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

## ► To cite this version:

Justin Blancon, Clément Buet, Pierre Dubreuil, Marie-Hélène Tixier, Frédéric Baret, et al.. Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions. TAG Theoretical and Applied Genetics, 2024, 137 (3), pp.1-20. 10.1007/s00122-024-04572-6 . hal-04495117

HAL Id: hal-04495117

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

Submitted on 8 Mar 2024

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
# Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions

Justin Blancon, Clément Buet, Pierre Dubreuil, Marie-Hélène Tixier, Frédéric Baret, Sébastien Praud

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

## ► To cite this version:

Justin Blancon, Clément Buet, Pierre Dubreuil, Marie-Hélène Tixier, Frédéric Baret, et al.. Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions. TAG Theoretical and Applied Genetics, 2024, 137 (3), pp.1-20. 10.1007/s00122-024-04572-6 . hal-04495117

HAL Id: hal-04495117

https://hal.inrae.fr/hal-04495117v1

Submitted on 8 Mar 2024

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>
# Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions

Justin Blancocon<fnref n="1" />,<fnref n="2" /> ID · Clément Buet<fnref n="2" /> · Pierre Dubreuil<fnref n="2" /> · Marie-Hélène Tixier<fnref n="2" /> · Frédéric Baret<fnref n="3" /> · Sébastien Praud<fnref n="2" />

Received: 7 December 2022 / Accepted: 3 February 2024

© The Author(s) 2024

## Abstract

**Key message** Green Leaf Area Index dynamics is a promising secondary trait for grain yield and drought tolerance. Multivariate GWAS is particularly well suited to identify the genetic determinants of the green leaf area index dynamics.

**Abstract** Improvement of maize grain yield is impeded by important genotype-environment interactions, especially under drought conditions. The use of secondary traits, that are correlated with yield, more heritable and less prone to genotype-environment interactions, can increase breeding efficiency. Here, we studied the genetic basis of a new secondary trait: the green leaf area index (GLAI) dynamics over the maize life cycle. For this, we used an unmanned aerial vehicle to characterize the GLAI dynamics of a diverse panel in well-watered and water-deficient trials in two years. From the dynamics, we derived 24 traits (slopes, durations, areas under the curve), and showed that six of them were heritable traits representative of the panel diversity. To identify the genetic determinants of GLAI, we compared two genome-wide association approaches: a univariate (single-trait) method and a multivariate (multi-trait) method combining GLAI traits, grain yield, and precocity. The explicit modeling of correlation structure between secondary traits and grain yield in the multivariate mixed model led to 2.5 times more associations detected. A total of 475 quantitative trait loci (QTLs) were detected. The genetic architecture of GLAI traits appears less complex than that of yield with stronger-effect QTLs that are more stable between environments. We also showed that a subset of GLAI QTLs explains nearly one fifth of yield variability across a larger environmental network of 11 water-deficient trials. GLAI dynamics is a promising grain yield secondary trait in optimal and drought conditions, and the detected QTLs could help to increase breeding efficiency through a marker-assisted approach.

## Introduction

Maize (Zea mays subsp. mays) is the most produced crop in the world, and a major staple food in several cultures. However maize yields are stagnating or decreasing in around 30% of its growing areas despite a continuous genetic progress (Duvick 2005; Tollenaar and Lee 2006; Ray et al. 2012; Laidig et al. 2014). Indeed, in some regions, the genetic

progress is partially countervailed by climate change (Lobell et al. 2011; Ray et al. 2019). Among climatic limiting factors, drought reduces maize yield of 5–6% globally (Lesk et al. 2016) and climate forecasts indicate that drought risk will increase in several maize growing regions around the world by the end of the century (Harrison et al. 2014; Hoegh-Guldberg et al. 2018; Pörtner et al. 2022) leading to even more serious losses (Webber et al. 2018; Leng and Hall 2019). Drought tolerance has therefore become a crucial goal for maize breeding.

Improving grain yield under drought stress is however a difficult task. Grain yield is a highly complex trait resulting from the interaction of the environment and the plant structural characteristics, physiological processes and regulatory pathways along the whole growth cycle. It thus presents large G×E (genotype×environment) interactions and a low heritability, especially under drought stress, that limit breeding efficiency. Conversely, physiological components of grain yield are simpler traits, which may be more heritable

Communicated by Thomas Lubberstedt.

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

Justin Blancocon
justin.blancocon@inrae.fr

<footnote n="1">¹ UMR GDEC, INRAE, Université Clermont Auvergne, 63000 Clermont-Ferrand, France</footnote>

<footnote n="2">² Biogemma, Centre de Recherche de Chappes, 63720 Chappes, France</footnote>

<footnote n="3">³ UMR EMMAH, UMT CAPTE, INRAE, 84914 Avignon, France</footnote>
# Maize green leaf area index dynamics: genetic basis of a new secondary trait for grain yield in optimal and drought conditions

Justin Blancoń<fnref n="1" />,<fnref n="2" /> ID · Clément Buet<fnref n="2" /> · Pierre Dubreuil<fnref n="2" /> · Marie-Hélène Tixier<fnref n="2" /> · Frédéric Baret<fnref n="3" /> · Sébastien Praud<fnref n="2" />

Received: 7 December 2022 / Accepted: 3 February 2024

© The Author(s) 2024

## Abstract

**Key message** Green Leaf Area Index dynamics is a promising secondary trait for grain yield and drought tolerance. Multivariate GWAS is particularly well suited to identify the genetic determinants of the green leaf area index dynamics.

**Abstract** Improvement of maize grain yield is impeded by important genotype-environment interactions, especially under drought conditions. The use of secondary traits, that are correlated with yield, more heritable and less prone to genotype-environment interactions, can increase breeding efficiency. Here, we studied the genetic basis of a new secondary trait: the green leaf area index (GLAI) dynamics over the maize life cycle. For this, we used an unmanned aerial vehicle to characterize the GLAI dynamics of a diverse panel in well-watered and water-deficient trials in two years. From the dynamics, we derived 24 traits (slopes, durations, areas under the curve), and showed that six of them were heritable traits representative of the panel diversity. To identify the genetic determinants of GLAI, we compared two genome-wide association approaches: a univariate (single-trait) method and a multivariate (multi-trait) method combining GLAI traits, grain yield, and precocity. The explicit modeling of correlation structure between secondary traits and grain yield in the multivariate mixed model led to 2.5 times more associations detected. A total of 475 quantitative trait loci (QTLs) were detected. The genetic architecture of GLAI traits appears less complex than that of yield with stronger-effect QTLs that are more stable between environments. We also showed that a subset of GLAI QTLs explains nearly one fifth of yield variability across a larger environmental network of 11 water-deficient trials. GLAI dynamics is a promising grain yield secondary trait in optimal and drought conditions, and the detected QTLs could help to increase breeding efficiency through a marker-assisted approach.

## Introduction

Maize (Zea mays subsp. mays) is the most produced crop in the world, and a major staple food in several cultures. However maize yields are stagnating or decreasing in around 30% of its growing areas despite a continuous genetic progress (Duvick 2005; Tollenaar and Lee 2006; Ray et al. 2012; Laidig et al. 2014). Indeed, in some regions, the genetic

progress is partially countervailed by climate change (Lobell et al. 2011; Ray et al. 2019). Among climatic limiting factors, drought reduces maize yield of 5–6% globally (Lesk et al. 2016) and climate forecasts indicate that drought risk will increase in several maize growing regions around the world by the end of the century (Harrison et al. 2014; Hoegh-Guldberg et al. 2018; Pörtner et al. 2022) leading to even more serious losses (Webber et al. 2018; Leng and Hall 2019). Drought tolerance has therefore become a crucial goal for maize breeding.

Improving grain yield under drought stress is however a difficult task. Grain yield is a highly complex trait resulting from the interaction of the environment and the plant structural characteristics, physiological processes and regulatory pathways along the whole growth cycle. It thus presents large G×E (genotype×environment) interactions and a low heritability, especially under drought stress, that limit breeding efficiency. Conversely, physiological components of grain yield are simpler traits, which may be more heritable

Communicated by Thomas Lubberstedt.

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

Justin Blancoń
justin.blancoń@inrae.fr

<footnote n="1">¹ UMR GDEC, INRAE, Université Clermont Auvergne, 63000 Clermont-Ferrand, France</footnote>

<footnote n="2">² Biogemma, Centre de Recherche de Chappes, 63720 Chappes, France</footnote>

<footnote n="3">³ UMR EMMAH, UMT CAPTE, INRAE, 84914 Avignon, France</footnote>
and less prone to G×E interaction. Using such secondary traits in breeding programs, based on a well-identified ideotype, could enhance selection efficiency when drought stress occurs and allow a better understanding of yield establishment in limiting conditions. To this aim, secondary traits must be correlated with yield under drought stress, be more heritable and less prone to G×E interaction than yield, and must be measurable at low cost in small plots in a nondestructive and rapid manner (Bolaños and Edmeades 1996; Bänziger and Lafitte 1997; Araus et al. 2012). Breeding programs incorporating secondary traits in addition to yield, through selection indices, have been shown to be particularly effective in stressful environments (Monneveux et al. 2006), and notably more effective than selection on yield alone (Fischer et al. 1989; Bänziger and Lafitte 1997). Accordingly, genomic selection studies have shown that using secondary traits in trait-assisted approaches could improve predictive ability compared to classical model (Rutkoski et al. 2016; Sun et al. 2017; Sandhu et al. 2021). Ribaut and Ragot (2007) and Beyene et al. (2016) also achieved substantial yield gains under drought conditions by using secondary traits QTLs in marker-assisted selection programs.

The green leaf area index (GLAI), the cumulative area of one leaf face per square meter of soil, plays a major role in light interception, transpiration and CO₂ exchanges. Its dynamics have thus a predominant impact on grain yield establishment (Monteith 1977; Passioura 1977). It corresponds to the integration, at the canopy level, of the radiation and water use efficiency over the plant cycle, which makes it a promising secondary trait to study drought tolerance. However, its measurement in the field have long been too tedious to use it as a secondary trait. Recent advances now make possible to phenotype hundreds of genotypes in the field for the whole GLAI dynamics (Blancon et al. 2019) and pave the way for its use as a secondary trait to identify and better understand the genetic determinants of drought tolerance.

The benefits of such an approach, both in terms of biological understanding and statistical power, depend on the efforts made to model the temporal evolution of the phenotype (Wu and Lin 2006). The most straightforward approach for detecting longitudinal quantitative trait loci (QTLs) is the discrete approach. It considers the different observations made during the growth cycle as independent traits that are each analyzed through a standard QTL model (Wu and Lin 2006; Hurtado et al. 2012; Würschum et al. 2014; Condorelli et al. 2018). However, the biological interpretation of such analysis is not always obvious because it is difficult to give sense to QTL affecting an arbitrarily chosen point across the curve (Malosetti et al. 2006; Hurtado et al. 2012). Furthermore, measurement dates are often determined by practical rather than biological considerations, and the comparison between different studies is difficult because it is unlikely

that measurements were made at identical time points. This method also assumes the independence of the measurements made at different times, which can lead to a loss of power for the QTL detection (Wu and Lin 2006; Malosetti et al. 2006). To account for the observations' covariance over time, it is possible to integrate them within a multivariate (i.e. multi-trait) model (Galesloot et al. 2014). However, this method, like the previous one, does not account for the temporal dimension and simply represents the studied phenotype as a set of correlated point measurements. It is also difficult to apply in practice because it requires the adjustment of a large number of parameters, and implies that each genotype is characterized at the same date, which is not always possible (Wu and Lin 2006; Malosetti et al. 2006).

A second family of approaches relies on a two-step analysis (Reymond et al. 2003; Hurtado et al. 2012; Crispim et al. 2015; Campbell et al. 2017). The first step consists in estimating a small number of genotypic parameters that summarize the trait dynamics. These parameters can be computed directly from the raw dynamics (slopes, durations, areas, etc.), or extracted from a temporal model most often using exponential, logistic, or Gombertz functions (van Eeuwijk et al. 2019). These genotypic parameters thus incorporate a temporal dimension and present the advantage to be easily interpretable at the biological level. The second step relies on the use of these parameters as response variables in a classical GWAS model. Nevertheless, because these genotypic parameters are often correlated, a joint analysis in a multivariate QTL detection model could better distinguish the effects of markers on each trait, and increase testing power (Campbell et al. 2017). Although this method has been shown to be effective, the quality of the results obtained remains dependent on the accuracy of the estimation of genotypic parameters (Verbeke and Molenberghs 2009). A loss of information between the two steps is also inevitable and may result in a loss of statistical power (Campbell et al. 2019).

Methods that combine temporal modeling and QTL detection in one step, grouped under the functional mapping terminology, have also been proposed, such as the use of mixture models or nonlinear mixed models (Ma et al. 2002; Wu and Lin 2006; Malosetti et al. 2006; Li et al. 2010). However, the use of parametric models may be too restrictive to model the temporal evolution of some traits, including leaf area and senescence (van Eeuwijk et al. 2019). Random regression, using covariance functions (spline or polynomial functions), then offers a nonparametric alternative to model the dynamics of the trait over time, while considering the temporal covariance between each measurement (Kirkpatrick et al. 1990; Wu and Lin 2006; Ning et al. 2017; Campbell et al. 2019). This approach has been widely used in animal genetics (Huisman et al. 2002; Kranis et al. 2007; Howard et al. 2015; Ning et al. 2017, 2018) but is still rare

<figure><img src="image_4.png" /></figure>
and less prone to G×E interaction. Using such secondary traits in breeding programs, based on a well-identified ideotype, could enhance selection efficiency when drought stress occurs and allow a better understanding of yield establishment in limiting conditions. To this aim, secondary traits must be correlated with yield under drought stress, be more heritable and less prone to G×E interaction than yield, and must be measurable at low cost in small plots in a nondestructive and rapid manner (Bolaños and Edmeades 1996; Bänziger and Lafitte 1997; Araus et al. 2012). Breeding programs incorporating secondary traits in addition to yield, through selection indices, have been shown to be particularly effective in stressful environments (Monneveux et al. 2006), and notably more effective than selection on yield alone (Fischer et al. 1989; Bänziger and Lafitte 1997). Accordingly, genomic selection studies have shown that using secondary traits in trait-assisted approaches could improve predictive ability compared to classical model (Rutkoski et al. 2016; Sun et al. 2017; Sandhu et al. 2021). Ribaut and Ragot (2007) and Beyene et al. (2016) also achieved substantial yield gains under drought conditions by using secondary traits QTLs in marker-assisted selection programs.

The green leaf area index (GLAI), the cumulative area of one leaf face per square meter of soil, plays a major role in light interception, transpiration and CO₂ exchanges. Its dynamics have thus a predominant impact on grain yield establishment (Monteith 1977; Passioura 1977). It corresponds to the integration, at the canopy level, of the radiation and water use efficiency over the plant cycle, which makes it a promising secondary trait to study drought tolerance. However, its measurement in the field have long been too tedious to use it as a secondary trait. Recent advances now make possible to phenotype hundreds of genotypes in the field for the whole GLAI dynamics (Blancon et al. 2019) and pave the way for its use as a secondary trait to identify and better understand the genetic determinants of drought tolerance.

The benefits of such an approach, both in terms of biological understanding and statistical power, depend on the efforts made to model the temporal evolution of the phenotype (Wu and Lin 2006). The most straightforward approach for detecting longitudinal quantitative trait loci (QTLs) is the discrete approach. It considers the different observations made during the growth cycle as independent traits that are each analyzed through a standard QTL model (Wu and Lin 2006; Hurtado et al. 2012; Würschum et al. 2014; Condorelli et al. 2018). However, the biological interpretation of such analysis is not always obvious because it is difficult to give sense to QTL affecting an arbitrarily chosen point across the curve (Malosetti et al. 2006; Hurtado et al. 2012). Furthermore, measurement dates are often determined by practical rather than biological considerations, and the comparison between different studies is difficult because it is unlikely

that measurements were made at identical time points. This method also assumes the independence of the measurements made at different times, which can lead to a loss of power for the QTL detection (Wu and Lin 2006; Malosetti et al. 2006). To account for the observations' covariance over time, it is possible to integrate them within a multivariate (i.e. multivariate) model (Galesloot et al. 2014). However, this method, like the previous one, does not account for the temporal dimension and simply represents the studied phenotype as a set of correlated point measurements. It is also difficult to apply in practice because it requires the adjustment of a large number of parameters, and implies that each genotype is characterized at the same date, which is not always possible (Wu and Lin 2006; Malosetti et al. 2006).

A second family of approaches relies on a two-step analysis (Reymond et al. 2003; Hurtado et al. 2012; Crispim et al. 2015; Campbell et al. 2017). The first step consists in estimating a small number of genotypic parameters that summarize the trait dynamics. These parameters can be computed directly from the raw dynamics (slopes, durations, areas, etc.), or extracted from a temporal model most often using exponential, logistic, or Gombertz functions (van Eeuwijk et al. 2019). These genotypic parameters thus incorporate a temporal dimension and present the advantage to be easily interpretable at the biological level. The second step relies on the use of these parameters as response variables in a classical GWAS model. Nevertheless, because these genotypic parameters are often correlated, a joint analysis in a multivariate QTL detection model could better distinguish the effects of markers on each trait, and increase testing power (Campbell et al. 2017). Although this method has been shown to be effective, the quality of the results obtained remains dependent on the accuracy of the estimation of genotypic parameters (Verbeke and Molenberghs 2009). A loss of information between the two steps is also inevitable and may result in a loss of statistical power (Campbell et al. 2019).

Methods that combine temporal modeling and QTL detection in one step, grouped under the functional mapping terminology, have also been proposed, such as the use of mixture models or nonlinear mixed models (Ma et al. 2002; Wu and Lin 2006; Malosetti et al. 2006; Li et al. 2010). However, the use of parametric models may be too restrictive to model the temporal evolution of some traits, including leaf area and senescence (van Eeuwijk et al. 2019). Random regression, using covariance functions (spline or polynomial functions), then offers a nonparametric alternative to model the dynamics of the trait over time, while considering the temporal covariance between each measurement (Kirkpatrick et al. 1990; Wu and Lin 2006; Ning et al. 2017; Campbell et al. 2019). This approach has been widely used in animal genetics (Huisman et al. 2002; Kranis et al. 2007; Howard et al. 2015; Ning et al. 2017, 2018) but is still rare

<figure><img src="image_4.png" /></figure>
in plant genetics. In plants, it is mainly applied in genomic selection (Sun et al. 2017; Campbell et al. 2018), with the exception of Campbell et al. (2019) that published the first GWAS using a random regression model performed on a major crop (rice). This approach is however highly computationally demanding (Ning et al. 2017; Moreira et al. 2020).

The main objective of the present study was to explore the genetic determinism of maize GLAI dynamics in drought and optimal conditions. We also aimed to better understand its link with grain yield, and to evaluate the potential of GLAI dynamics as a secondary trait in selection. For this, we evaluated a multi-parent advanced generation inter-cross (MAGIC) panel in four environments with contrasted hydric regimes. We also phenotyped the GLAI dynamics over the whole maize life cycle using a multispectral camera mounted on a UAV (unmanned aerial vehicle). A simple physiological GLAI model was then used to derive the parameters of the dynamics. To appraise the potential of GLAI dynamics as a secondary trait in selection, its sensitivity to G×E interaction was compared to that of yield and its components. We then applied a univariate and a multivariate two-step GWAS approach based on these parameters to identify and characterize the GLAI dynamics genetic determinants.

# Materials and methods

## Plant material, genotyping and relatedness

We used a broad diversity panel of doubled haploid (DH) maize lines derived from a MAGIC population. This population originated from a funnel cross with 16 parents chosen among historical lines that are representative of the genetic diversity of temperate material (Buet et al. 2013). To limit confounding effects due to precocity differences, certain DH lines were discarded to narrow the range of flowering times, with 95% of the panel flowering within an eight-day period when evaluated in a hybrid context.

The 16 parents of the panel were sequenced with the Illumina HiSeq 2000 at a 12–16 X depth of coverage. The 324 DH lines were genotyped with the Illumina 50 K array (Ganal et al. 2011) and the 600 K Affymetrix array (Unterseer et al. 2014). Using the genotyping data, parental haplotypes were reconstructed with a modified version of the *R* package *qtl* (Broman et al. 2003; R Core Team 2017), then the parental origin was inferred for each allele in every DH line. Based on parental mosaic and sequencing data for each of the 16 parents, the panel was genotyped in silico, resulting in around 17 million single nucleotide polymorphisms (SNPs), that were narrowed down to 7,806,995 SNPs after filtering out any with missing data in parents or a low minor allele frequency (3.5% threshold).

The genomic relatedness between individuals was estimated from the Illumina 50 K array genotyping data. Removal of SNPs with missing data or at biased SYNGENTA loci (Ganal et al. 2011) resulted in 32,466 SNPs evenly distributed on the 10 chromosomes, which were then used to compute a global identity-by-state (IBS) kinship matrix and 10 leave-one-chromosome-out (LOCO) IBS kinship matrices (Rincent et al. 2014) in the *R* package *GAPIT* (Lipka et al. 2012).

## Agronomic evaluation of maize hybrids

DH lines were crossed with the tester line MBS847 and 360 of the resultant hybrids were evaluated in the field in 2016 and 347 were evaluated in 2017, 330 being common to both years. Each year, a well-watered (WW, irrigated) trial and a water-deficient (WD, rain-fed with monitored irrigation) trial were carried out in Saint-Paul-lès-Romans, close to Romans-sur-Isère, France (45°04'06"N, 5°08'01"E). The trials were laid out as alpha-lattice designs with two replicate blocks in each treatment. The microplots were sown at a density of 9 plants/m² with two rows per microplot, a 0.8 m row spacing and a row length of 6 m. The combination of the two years and two water treatments were considered as four environments 16STPAULWD, 16STPAULWW, 17STPAULWD, and 17STPAULWW. More details about climatic conditions and agronomic practices are given in Blancon et al. (2019).

In the 16STPAULWD environment, as evidenced by water balance analysis (Fig. S1), drought stress occurred for a short period during flowering and more severely during grain filling leading to 38% grain yield loss compared to 16STPAULWW. In the 17STPAULWD environment, stress only occurred during flowering (Fig. S1) leading to 20% grain yield loss compared to 17STPAULWW. In each trial, we measured grain yield (GY, in q.ha⁻¹) and thousand kernel weight (TKW, in g), both adjusted to 15% humidity. We computed kernel number per square meter (KN, in grain.m⁻²), and scored the female flowering date (FF, in °C.j) and the anthesis to silking interval (ASI, in °C.j).

Between 2014 and 2017, the performance of the maize hybrid panel was evaluated in 11 field trials under WD conditions in five locations: Blois, Saint-Paul-lès-Romans and Nérac in France, Szeged in Hungary and Graneros in Chile. More details on trial conditions can be found in Table S1. These 11 environments include 16STPAULWD and 17STPAULWD, described above. GY data was obtained in each of these 11 experiments (GY₁₁).
in plant genetics. In plants, it is mainly applied in genomic selection (Sun et al. 2017; Campbell et al. 2018), with the exception of Campbell et al. (2019) that published the first GWAS using a random regression model performed on a major crop (rice). This approach is however highly computationally demanding (Ning et al. 2017; Moreira et al. 2020).

The main objective of the present study was to explore the genetic determinism of maize GLAI dynamics in drought and optimal conditions. We also aimed to better understand its link with grain yield, and to evaluate the potential of GLAI dynamics as a secondary trait in selection. For this, we evaluated a multi-parent advanced generation inter-cross (MAGIC) panel in four environments with contrasted hydric regimes. We also phenotyped the GLAI dynamics over the whole maize life cycle using a multispectral camera mounted on a UAV (unmanned aerial vehicle). A simple physiological GLAI model was then used to derive the parameters of the dynamics. To appraise the potential of GLAI dynamics as a secondary trait in selection, its sensitivity to G×E interaction was compared to that of yield and its components. We then applied a univariate and a multivariate two-step GWAS approach based on these parameters to identify and characterize the GLAI dynamics genetic determinants.

# Materials and methods

## Plant material, genotyping and relatedness

We used a broad diversity panel of doubled haploid (DH) maize lines derived from a MAGIC population. This population originated from a funnel cross with 16 parents chosen among historical lines that are representative of the genetic diversity of temperate material (Buet et al. 2013). To limit confounding effects due to precocity differences, certain DH lines were discarded to narrow the range of flowering times, with 95% of the panel flowering within an eight-day period when evaluated in a hybrid context.

The 16 parents of the panel were sequenced with the Illumina HiSeq 2000 at a 12–16 X depth of coverage. The 324 DH lines were genotyped with the Illumina 50 K array (Ganal et al. 2011) and the 600 K Affymetrix array (Unterseer et al. 2014). Using the genotyping data, parental haplotypes were reconstructed with a modified version of the *R* package *qtl* (Broman et al. 2003; R Core Team 2017), then the parental origin was inferred for each allele in every DH line. Based on parental mosaic and sequencing data for each of the 16 parents, the panel was genotyped in silico, resulting in around 17 million single nucleotide polymorphisms (SNPs), that were narrowed down to 7,806,995 SNPs after filtering out any with missing data in parents or a low minor allele frequency (3.5% threshold).

The genomic relatedness between individuals was estimated from the Illumina 50 K array genotyping data. Removal of SNPs with missing data or at biased SYNGENTA loci (Ganal et al. 2011) resulted in 32,466 SNPs evenly distributed on the 10 chromosomes, which were then used to compute a global identity-by-state (IBS) kinship matrix and 10 leave-one-chromosome-out (LOCO) IBS kinship matrices (Rincent et al. 2014) in the *R* package *GAPIT* (Lipka et al. 2012).

## Agronomic evaluation of maize hybrids

DH lines were crossed with the tester line MBS847 and 360 of the resultant hybrids were evaluated in the field in 2016 and 347 were evaluated in 2017, 330 being common to both years. Each year, a well-watered (WW, irrigated) trial and a water-deficient (WD, rain-fed with monitored irrigation) trial were carried out in Saint-Paul-lès-Romans, close to Romans-sur-Isère, France (45°04'06"N, 5°08'01"E). The trials were laid out as alpha-lattice designs with two replicate blocks in each treatment. The microplots were sown at a density of 9 plants/m² with two rows per microplot, a 0.8 m row spacing and a row length of 6 m. The combination of the two years and two water treatments were considered as four environments 16STPAULWD, 16STPAULWW, 17STPAULWD, and 17STPAULWW. More details about climatic conditions and agronomic practices are given in Blancon et al. (2019).

In the 16STPAULWD environment, as evidenced by water balance analysis (Fig. S1), drought stress occurred for a short period during flowering and more severely during grain filling leading to 38% grain yield loss compared to 16STPAULWW. In the 17STPAULWD environment, stress only occurred during flowering (Fig. S1) leading to 20% grain yield loss compared to 17STPAULWW. In each trial, we measured grain yield (GY, in q.ha⁻¹) and thousand kernel weight (TKW, in g), both adjusted to 15% humidity. We computed kernel number per square meter (KN, in grain.m⁻²), and scored the female flowering date (FF, in °C.j) and the anthesis to silking interval (ASI, in °C.j).

Between 2014 and 2017, the performance of the maize hybrid panel was evaluated in 11 field trials under WD conditions in five locations: Blois, Saint-Paul-lès-Romans and Nérac in France, Szeged in Hungary and Graneros in Chile. More details on trial conditions can be found in Table S1. These 11 environments include 16STPAULWD and 17STPAULWD, described above. GY data was obtained in each of these 11 experiments (GY₁₁).
# Characterization of GLAI dynamics and derived traits

## UAV model-assisted phenotyping of GLAI dynamics

We flew a hexacopter UAV mounted with a six-band multispectral camera nine times over the 2016 trial plots and eleven times over the 2017 trial plots. The multispectral camera is composed of six cameras with 10-nm spectral resolution bands centered, respectively, at 450, 532, 568, 675, 730, and 850 nm. Images were recorded continuously during the flight at a frequency of 1 Hz, and the integration time was adjusted automatically to minimize saturation and maximize the range of variation. The reflectances recorded in the first five bands were normalized by dividing the values by the reflectance measured in the near-infrared to limit the effect of illumination variations during the flight. More details about the processing of multispectral images are given in Fig. S2 and in Blancon et al. (2019).

For each flight, we calibrated empirical relationships linking ground level measurements of GLAI achieved over a small sample of microplots to normalized reflectances.

We then estimated GLAI value at each time point from these transfer functions for each microplot of the trials.

Finally, we used a simple, physiologically-based dynamics model to fit the GLAI estimates and predict the continuous GLAI dynamics for each microplot over the whole season in steps of one growing degree day with a 6 °C base (GDD6). More details can be found in Blancon et al. (2019).

## Estimation of traits describing GLAI dynamics

To describe the GLAI dynamics, we derived integrative traits over five physiological distinct stages: an early vegetative (EV), late vegetative (LV), flowering (F), slow senescence (SS) and rapid senescence phase (RS).

We delimited the first two phases with a bent-cable regression (R package *SiZer*, Sonderegger (2022)) fitted to the GLAI dynamics between shoot emergence (when GLAI becomes positive) and the time when GLAI achieves its maximum (Fig. 1, $S_{EV}$ and $S_{LV}$). The abscissa of the intersection between the two linear parts of the regression defines the limit between the EV and LV phases. The flowering phase extends from the time when GLAI reaches its maximum until the onset of senescence, here defined

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Estimation of 24 parameters derived from maize GLAI dynamics. A Two bent-cable regressions were fitted, one during the vegetative phase and the other during the senescence phase, to estimate four slopes ($S_{EV}$, $S_{LV}$, $S_{SS}$, $S_{RS}$) and to define five distinct phases during the development cycle. B The phases previously defined were used to derive eight areas under the curve ($AUC_{EV}$, $AUC_{LV}$, $AUC_{V}$,</figcaption></figure>

$AUC_{F}$, $AUC_{SS}$, $AUC_{RS}$, $AUC_{S}$, $AUC_{C}$) and the eight corresponding durations ($D_{EV}$, $D_{LV}$, $D_{V}$, $D_{F}$, $D_{SS}$, $D_{RS}$, $D_{S}$, $D_{C}$). Finally, the maximum GLAI reached during the cycle (GLAI$_M$) was extracted from the model and then used to derive three additional durations (C: $D_{75}$, $D_{50}$, $D_{25}$) (color figure online)
# Characterization of GLAI dynamics and derived traits

## UAV model-assisted phenotyping of GLAI dynamics

We flew a hexacopter UAV mounted with a six-band multispectral camera nine times over the 2016 trial plots and eleven times over the 2017 trial plots. The multispectral camera is composed of six cameras with 10-nm spectral resolution bands centered, respectively, at 450, 532, 568, 675, 730, and 850 nm. Images were recorded continuously during the flight at a frequency of 1 Hz, and the integration time was adjusted automatically to minimize saturation and maximize the range of variation. The reflectances recorded in the first five bands were normalized by dividing the values by the reflectance measured in the near-infrared to limit the effect of illumination variations during the flight. More details about the processing of multispectral images are given in Fig. S2 and in Blancon et al. (2019).

For each flight, we calibrated empirical relationships linking ground level measurements of GLAI achieved over a small sample of microplots to normalized reflectances.

We then estimated GLAI value at each time point from these transfer functions for each microplot of the trials.

Finally, we used a simple, physiologically-based dynamics model to fit the GLAI estimates and predict the continuous GLAI dynamics for each microplot over the whole season in steps of one growing degree day with a 6 °C base (GDD6). More details can be found in Blancon et al. (2019).

## Estimation of traits describing GLAI dynamics

To describe the GLAI dynamics, we derived integrative traits over five physiological distinct stages: an early vegetative (EV), late vegetative (LV), flowering (F), slow senescence (SS) and rapid senescence phase (RS).

We delimited the first two phases with a bent-cable regression (R package *SiZer*, Sonderegger (2022)) fitted to the GLAI dynamics between shoot emergence (when GLAI becomes positive) and the time when GLAI achieves its maximum (Fig. 1, $S_{EV}$ and $S_{LV}$). The abscissa of the intersection between the two linear parts of the regression defines the limit between the EV and LV phases. The flowering phase extends from the time when GLAI reaches its maximum until the onset of senescence, here defined

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Estimation of 24 parameters derived from maize GLAI dynamics. A Two bent-cable regressions were fitted, one during the vegetative phase and the other during the senescence phase, to estimate four slopes ($S_{EV}$, $S_{LV}$, $S_{SS}$, $S_{RS}$) and to define five distinct phases during the development cycle. B The phases previously defined were used to derive eight areas under the curve ($AUC_{EV}$, $AUC_{LV}$, $AUC_{V}$,</figcaption></figure>

$AUC_{F}$, $AUC_{SS}$, $AUC_{RS}$, $AUC_{S}$, $AUC_{C}$) and the eight corresponding durations ($D_{EV}$, $D_{LV}$, $D_{V}$, $D_{F}$, $D_{SS}$, $D_{RS}$, $D_{S}$, $D_{C}$). Finally, the maximum GLAI reached during the cycle ($GLAI_{M}$) was extracted from the model and then used to derive three additional durations (C: $D_{75}$, $D_{50}$, $D_{25}$) (color figure online)
as the time when the GLAI becomes less than 95% of the maximum GLAI. We delimited the last two phases by fitting a bent-cable regression to the GLAI dynamics between the start of senescence and complete senescence (when the GLAI becomes zero) (Fig. 1, $S_{SS}$ and $S_{RS}$). The abscissa of the intersection between the two linear parts of this regression delimits the SS from the RS phase. A set of 24 traits counting slopes, durations, and areas under the curve were finally computed from these five phases, alone or combined (Fig. 1). These traits were extracted for each microplot of the four trials.

## Adjusted means, heritability, and variance decomposition

For each combination of trait, year and treatment, we used a linear mixed model with block effect, genotype effect and autocorrelated errors (Model S1) to estimate the adjusted means. The same model, with random genotype effect, was used to compute generalized heritability (Eq. S1, Cullis et al. 2006).

To analyze variance components at the network scale, we used a multi-environmental model with random genotype and $G \times E$ interaction effects (Model S2). Each random effect was tested with a restricted likelihood-ratio test comparing the full model and the null model without the considered effect. The variance components were extracted to be compared, and the standard deviation was expressed as a percentage of the mean $\mu$ ($CV_g$). All linear models were fitted using the *R* package *asreml-R* v3.2 (Butler et al. 2009).

## Trait correlation analysis

To get an overall picture of correlations between traits at the network scale, we performed a principal component analysis (PCA) with the *R* package *FactoMiner* (Lê et al. 2008), based on standardized adjusted means of GLAI and agronomic traits from the four environments. From the PCA results, six GLAI traits were retained for further analysis because they were representative of the panel diversity and characterize complementary biological processes: $D_V$, $D_S$, $GLAI_M$, $AUC_V$, $AUC_F$ and $AUC_S$, the duration of the vegetative phase, the duration of the senescence phase, the maximum of GLAI, the area under the curve during the vegetative, the flowering and the senescence phase respectively.

To better understand the link between GLAI dynamics and GY we evaluated the phenotypic correlation ($r_p$) by computing the Pearson correlation coefficient between these six GLAI traits and yield components, and tested it with a Student test in each of the four environments.

## Genome wide association studies

To identify the genetic determinants of GLAI dynamics, we used two GWAS methods: a univariate approach ($M_{UV}$) and a multivariate approach ($M_{MV}$). For both approaches, the LOCO method was used to avoid proximal contamination (Rincent et al. 2014).

### Univariate approach

A first association study was done trait by trait for each environment of the network, using an univariate model. The six chosen GLAI traits and the five agronomic traits were analyzed with this univariate approach. Model S3 was fitted with *FaST-LMM* v2.07 (Lippert et al. 2011), and the null hypothesis $H_0: \beta = 0$ was tested with a likelihood-ratio test.

### Multivariate approach

In the second association study, we analyzed several GLAI and agronomic traits jointly in the same multivariate model for each environment. Model S4 was fitted using *GEMMA* software (Zhou and Stephens 2014) for $t=8$ traits namely $AUC_V$, $AUC_F$, $AUC_S$, $GLAI_M$, $D_V$, $D_S$, GY and FF. Because of convergence issues, it was not possible to include other agronomic traits in the model. The null hypothesis $H_0: \beta_1 = \cdots = \beta_d = 0$ was tested with a Wald test. *GEMMA* was also used to estimate genetic correlations ($r_g$) between these eight traits.

### SNP clustering

To correct for multiple testing and account for the false discovery rate (Storey 2002), we computed a *qvalue* for each test with the *R* package *qvalue* (Storey et al. 2021). SNPs with *q* value $\le 0.05$ were considered as significantly associated.

For the multivariate approach, rejecting the null hypothesis ($H_0: \beta_1 = \cdots = \beta_d = 0$) does not give information about the trait(s) associated with the significant SNP. Significant SNPs effect for each trait were estimated with their standard deviation. The significance of these effects was tested a posteriori with a Wald test ($\alpha = 0.01$) for each trait, to assign each significant SNP to its trait(s).

For an environment-chromosome-trait combination, associated SNP were clustered with the method proposed by Cormier et al. (2014) with a linkage disequilibrium threshold of $R^2 \ge 0.6$. For a given combination, overlapping clusters were aggregated. Three different types of clusters were computed: $M_{UV}$ clusters built from SNP detected with $M_{UV}$, $M_{MV}$ clusters built from SNP detected with $M_{MV}$; and $M_{ALL}$ clusters built from SNP detected with $M_{UV}$ or $M_{MV}$. Clusters longer than 30% of the chromosome length were discarded as artifacts.
as the time when the GLAI becomes less than 95% of the maximum GLAI. We delimited the last two phases by fitting a bent-cable regression to the GLAI dynamics between the start of senescence and complete senescence (when the GLAI becomes zero) (Fig. 1, $S_{SS}$ and $S_{RS}$). The abscissa of the intersection between the two linear parts of this regression delimits the SS from the RS phase. A set of 24 traits counting slopes, durations, and areas under the curve were finally computed from these five phases, alone or combined (Fig. 1). These traits were extracted for each microplot of the four trials.

## Adjusted means, heritability, and variance decomposition

For each combination of trait, year and treatment, we used a linear mixed model with block effect, genotype effect and autocorrelated errors (Model S1) to estimate the adjusted means. The same model, with random genotype effect, was used to compute generalized heritability (Eq. S1, Cullis et al. 2006).

To analyze variance components at the network scale, we used a multi-environmental model with random genotype and $G \times E$ interaction effects (Model S2). Each random effect was tested with a restricted likelihood-ratio test comparing the full model and the null model without the considered effect. The variance components were extracted to be compared, and the standard deviation was expressed as a percentage of the mean $\mu$ ($CV_g$). All linear models were fitted using the *R* package *asreml-R* v3.2 (Butler et al. 2009).

## Trait correlation analysis

To get an overall picture of correlations between traits at the network scale, we performed a principal component analysis (PCA) with the *R* package *FactoMiner* (Lê et al. 2008), based on standardized adjusted means of GLAI and agronomic traits from the four environments. From the PCA results, six GLAI traits were retained for further analysis because they were representative of the panel diversity and characterize complementary biological processes: $D_V$, $D_S$, $GLAI_M$, $AUC_V$, $AUC_F$ and $AUC_S$, the duration of the vegetative phase, the duration of the senescence phase, the maximum of GLAI, the area under the curve during the vegetative, the flowering and the senescence phase respectively.

To better understand the link between GLAI dynamics and GY we evaluated the phenotypic correlation ($r_p$) by computing the Pearson correlation coefficient between these six GLAI traits and yield components, and tested it with a Student test in each of the four environments.

## Genome wide association studies

To identify the genetic determinants of GLAI dynamics, we used two GWAS methods: a univariate approach ($M_{UV}$) and a multivariate approach ($M_{MV}$). For both approaches, the LOCO method was used to avoid proximal contamination (Rincent et al. 2014).

### Univariate approach

A first association study was done trait by trait for each environment of the network, using an univariate model. The six chosen GLAI traits and the five agronomic traits were analyzed with this univariate approach. Model S3 was fitted with *FaST-LMM* v2.07 (Lippert et al. 2011), and the null hypothesis $H_0: \beta = 0$ was tested with a likelihood-ratio test.

### Multivariate approach

In the second association study, we analyzed several GLAI and agronomic traits jointly in the same multivariate model for each environment. Model S4 was fitted using *GEMMA* software (Zhou and Stephens 2014) for $t=8$ traits namely $AUC_V$, $AUC_F$, $AUC_S$, $GLAI_M$, $D_V$, $D_S$, GY and FF. Because of convergence issues, it was not possible to include other agronomic traits in the model. The null hypothesis $H_0: \beta_1 = \cdots = \beta_d = 0$ was tested with a Wald test. *GEMMA* was also used to estimate genetic correlations ($r_g$) between these eight traits.

### SNP clustering

To correct for multiple testing and account for the false discovery rate (Storey 2002), we computed a *qvalue* for each test with the *R* package *qvalue* (Storey et al. 2021). SNPs with *q* value $\le 0.05$ were considered as significantly associated.

For the multivariate approach, rejecting the null hypothesis ($H_0: \beta_1 = \cdots = \beta_d = 0$) does not give information about the trait(s) associated with the significant SNP. Significant SNPs effect for each trait were estimated with their standard deviation. The significance of these effects was tested a posteriori with a Wald test ($\alpha = 0.01$) for each trait, to assign each significant SNP to its trait(s).

For an environment-chromosome-trait combination, associated SNP were clustered with the method proposed by Cormier et al. (2014) with a linkage disequilibrium threshold of $R^2 \ge 0.6$. For a given combination, overlapping clusters were aggregated. Three different types of clusters were computed: $M_{UV}$ clusters built from SNP detected with $M_{UV}$, $M_{MV}$ clusters built from SNP detected with $M_{MV}$; and $M_{ALL}$ clusters built from SNP detected with $M_{UV}$ or $M_{MV}$. Clusters longer than 30% of the chromosome length were discarded as artifacts.
Each cluster was finally described by its most strongly associ-
ated SNP (peak SNP with the smallest *pvalue*) and the interval
covered by its SNP components.

QTL identification

We used a multi-environment multilocus backward selection
model to identify SNP clusters playing a major role in the
genetic determinism of each trait in the trial network (Model
S5). Initially, the effects of all candidate clusters were incor-
porated in Model S5. The effect of each cluster was tested with
a Wald test when it was added last to the model. The cluster
with the least significant effect was iteratively removed, until
all remaining clusters had been declared significant (α = 0.01).
After the backward elimination, the remaining clusters were
called QTL. QTL effects were estimated in each environ-
ment with Model S5, and tested a posteriori with a Wald test
(α = 0.01) to identify the environment(s) in which each QTL
significantly affects the considered traits. Finally, we defined
the QTL limit as the initial cluster interval extended by the
local linkage disequilibrium extent at each extremity. QTL
effects were defined as B73 allele effect, a parental line of
the panel.

QTL colocalization

To decipher the genetic link between the traits, we studied
QTL colocalization for every combination of trait and environ-
ment. Colocalization was defined as, at least, two overlapping
QTL.

Grain yield prediction from GLAI QTL

A multi-environment model was used to confirm the utility of
GLAI as a yield secondary trait under drought stress condi-
tions using GY₁₁ data from the eleven-environment WD trial
("Agronomic evaluation of maize hybrids" section). We quan-
tified the proportion of GY₁₁ variance explained by the GLAI
QTLs with a backward selection approach. This led to a list of
GLAI QTLs that significantly impact GY at the GY₁₁ network
scale. As this approach was applied to the GLAI QTLs result-
ing from an initial backward elimination ("QTL identification"
section), and not to all the clusters detected, a less conservative
threshold was applied (α=0.05). A Wald test (α = 0.05) was
used a posteriori to identify which QTL were significant in
each environment.

Results

Six parameters of GLAI dynamics are good
descriptors of maize diversity and environmental
constraints

To get a global picture of GLAI and agronomic traits cor-
relation in the trial network, we ran a principal component
analysis of the 24 GLAI related and the 5 agronomic traits
in the four environments together (Fig. 2). The first two PCA
axes explained respectively 40% and 28% of the panel diver-
sity. The first axis differentiated genotypes with a long veg-
etative phase (D<sub>V</sub>) from those with a long senescence phase
(D<sub>S</sub>, S<sub>SS</sub>). The second axis was more related to the GLAI
maximum (GLAI<sub>M</sub>) and the flowering phase (AUC<sub>F</sub>, AUC<sub>C</sub>).
In this plane, the GLAI traits are better represented than the
agronomic traits. Six traits—D<sub>V</sub>, D<sub>S</sub>, GLAI<sub>M</sub>, AUC<sub>V</sub>, AUC
<sub>F</sub> and AUC<sub>S</sub>—are particularly well correlated with the axes,
each showing a Pearson r<sup>2</sup> > 0.57 with one of the two axes.
These traits also captured most of the information carried
by the other GLAI traits at the network scale, for example,
D<sub>S</sub> for D<sub>SS</sub>, D<sub>RS</sub>, D<sub>F</sub>, D<sub>25</sub>, D<sub>50</sub>, and D<sub>75</sub>, or AUC<sub>S</sub> for AUC
<sub>C</sub>, AUC<sub>SS</sub>, AUC<sub>RS</sub> and S<sub>LV</sub> (Fig. 2). The aforementioned six
traits were therefore chosen as good descriptors of the GLAI

<figure><img src="image_6.png" /><figcaption>**Fig. 2** Correlation circle from principal component analysis of 29 maize traits. A normalized PCA was done on the adjusted means for five agronomic traits (green arrows) and 24 GLAI traits (black and red arrows) in 324 maize hybrids in four environments. The first two axes (Dim) are plotted together with their proportion of explained variance. For trait abbreviations see the “Material and methods” section and Fig. 1. Six GLAI traits were chosen as representing the panel diversity (red arrows) (color figure online)</figcaption></figure>
Each cluster was finally described by its most strongly associ-
ated SNP (peak SNP with the smallest *pvalue*) and the interval
covered by its SNP components.

QTL identification

We used a multi-environment multilocus backward selection
model to identify SNP clusters playing a major role in the
genetic determinism of each trait in the trial network (Model
S5). Initially, the effects of all candidate clusters were incor-
porated in Model S5. The effect of each cluster was tested with
a Wald test when it was added last to the model. The cluster
with the least significant effect was iteratively removed, until
all remaining clusters had been declared significant (α = 0.01).
After the backward elimination, the remaining clusters were
called QTL. QTL effects were estimated in each environ-
ment with Model S5, and tested a posteriori with a Wald test
(α = 0.01) to identify the environment(s) in which each QTL
significantly affects the considered traits. Finally, we defined
the QTL limit as the initial cluster interval extended by the
local linkage disequilibrium extent at each extremity. QTL
effects were defined as B73 allele effect, a parental line of
the panel.

QTL colocalization

To decipher the genetic link between the traits, we studied
QTL colocalization for every combination of trait and environ-
ment. Colocalization was defined as, at least, two overlapping
QTL.

Grain yield prediction from GLAI QTL

A multi-environment model was used to confirm the utility of
GLAI as a yield secondary trait under drought stress condi-
tions using GY₁₁ data from the eleven-environment WD trial
("Agronomic evaluation of maize hybrids" section). We quan-
tified the proportion of GY₁₁ variance explained by the GLAI
QTLs with a backward selection approach. This led to a list of
GLAI QTLs that significantly impact GY at the GY₁₁ network
scale. As this approach was applied to the GLAI QTLs result-
ing from an initial backward elimination ("QTL identification"
section), and not to all the clusters detected, a less conservative
threshold was applied (α=0.05). A Wald test (α = 0.05) was
used a posteriori to identify which QTL were significant in
each environment.

Results

Six parameters of GLAI dynamics are good
descriptors of maize diversity and environmental
constraints

To get a global picture of GLAI and agronomic traits cor-
relation in the trial network, we ran a principal component
analysis of the 24 GLAI related and the 5 agronomic traits
in the four environments together (Fig. 2). The first two PCA
axes explained respectively 40% and 28% of the panel diver-
sity. The first axis differentiated genotypes with a long veg-
etative phase (D<sub>V</sub>) from those with a long senescence phase
(D<sub>S</sub>, S<sub>SS</sub>). The second axis was more related to the GLAI
maximum (GLAI<sub>M</sub>) and the flowering phase (AUC<sub>F</sub>, AUC<sub>C</sub>).
In this plane, the GLAI traits are better represented than the
agronomic traits. Six traits—D<sub>V</sub>, D<sub>S</sub>, GLAI<sub>M</sub>, AUC<sub>V</sub>, AUC
<sub>F</sub> and AUC<sub>S</sub>—are particularly well correlated with the axes,
each showing a Pearson r<sup>2</sup> > 0.57 with one of the two axes.
These traits also captured most of the information carried
by the other GLAI traits at the network scale, for example,
D<sub>S</sub> for D<sub>SS</sub>, D<sub>RS</sub>, D<sub>F</sub>, D<sub>25</sub>, D<sub>50</sub>, and D<sub>75</sub>, or AUC<sub>S</sub> for AUC
<sub>C</sub>, AUC<sub>SS</sub>, AUC<sub>RS</sub> and S<sub>LV</sub> (Fig. 2). The aforementioned six
traits were therefore chosen as good descriptors of the GLAI

<figure><img src="image_6.png" /><figcaption>**Fig. 2** Correlation circle from principal component analysis of 29 maize traits. A normalized PCA was done on the adjusted means for five agronomic traits (green arrows) and 24 GLAI traits (black and red arrows) in 324 maize hybrids in four environments. The first two axes (Dim) are plotted together with their proportion of explained variance. For trait abbreviations see the “Material and methods” section and Fig. 1. Six GLAI traits were chosen as representing the panel diversity (red arrows) (color figure online)</figcaption></figure>
diversity in the panel. Of the agronomic traits, GY is the best represented, mainly by the second component ($r^2=0.39$). GY is thus more correlated to the $GLAI_M$ and the AUC traits than to the duration of the vegetative or senescence phase. Conversely, grain yield components TKW and KN were more weakly linked with the GLAI traits, while flowering traits FF and ASI, are poorly related to GLAI dynamics and GY on the whole network.

Environmental conditions affected differently GLAI and agronomics traits in 2016 and 2017 (Fig. 3). Drought stress affected agronomic traits in the same direction in 2016 and 2017, with decreases in GY, KN and TKW, and a delay in FF, leading to a slight increase in ASI. TKW was barely impacted in 2017 by the drought stress as it

occurred during the flowering phase. GLAI traits were more affected by the different timing of the drought stress in 2016 and 2017, which reveal contrasting physiological constraints between the two years. GLAI traits related to the vegetative phase were less affected by drought stress than those related to the magnitude or senescence phase, perhaps because drought had no impact on the duration of the vegetative phase ($D_V$ Fig. 3, $D_{EV}$ and $D_{LV}$ Fig. S4). Indeed, drought stress began a few days before flowering in 2016 and 2017, so had only a marginal impact on vegetative processes. However, drought stress was early enough to decrease the dynamics amplitude ($GLAI_M$) in 2017, mainly by reducing leaf area (Blancon et al. 2019), and thus $AUC_V$ and $AUC_S$. The lower $GLAI_M$ in

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Drought stress impact on five agronomic traits and six modeled GLAI traits representative of maize diversity. The agronomic traits are GY (q.ha⁻¹), KN (kernels.m⁻²), TKW (g), FF (GDD6) and ASI (GDD6). The GLAI traits are $D_V$ (GDD6), $D_S$ (GDD6), $GLAI_M$ (m².m⁻²), $AUC_V$ (m².m⁻².GDD6), $AUC_F$ (m².m⁻².GDD6), and $AUC_S$ (m².m⁻².GDD6). Boxplots were constructed from the adjusted means</figcaption></figure>

of traits measured in 324 maize hybrids grown in the field in well-watered (WW) and water-deficient (WD) conditions. Boxes represent the inter-quartile range and whiskers extend to 1.5 times the inter-quartile range. Medians (horizontal lines inside the box), means (diamond symbols), and outliers (circle symbols) are indicated for each trait in each environment (color figure online)
diversity in the panel. Of the agronomic traits, GY is the best represented, mainly by the second component ($r^2=0.39$). GY is thus more correlated to the $GLAI_M$ and the AUC traits than to the duration of the vegetative or senescence phase. Conversely, grain yield components TKW and KN were more weakly linked with the GLAI traits, while flowering traits FF and ASI, are poorly related to GLAI dynamics and GY on the whole network.

Environmental conditions affected differently GLAI and agronomics traits in 2016 and 2017 (Fig. 3). Drought stress affected agronomic traits in the same direction in 2016 and 2017, with decreases in GY, KN and TKW, and a delay in FF, leading to a slight increase in ASI. TKW was barely impacted in 2017 by the drought stress as it

occurred during the flowering phase. GLAI traits were more affected by the different timing of the drought stress in 2016 and 2017, which reveal contrasting physiological constraints between the two years. GLAI traits related to the vegetative phase were less affected by drought stress than those related to the magnitude or senescence phase, perhaps because drought had no impact on the duration of the vegetative phase ($D_V$ Fig. 3, $D_{EV}$ and $D_{LV}$ Fig. S4). Indeed, drought stress began a few days before flowering in 2016 and 2017, so had only a marginal impact on vegetative processes. However, drought stress was early enough to decrease the dynamics amplitude ($GLAI_M$) in 2017, mainly by reducing leaf area (Blancon et al. 2019), and thus $AUC_V$ and $AUC_S$. The lower $GLAI_M$ in

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Drought stress impact on five agronomic traits and six modeled GLAI traits representative of maize diversity. The agronomic traits are GY (q.ha⁻¹), KN (kernels.m⁻²), TKW (g), FF (GDD6) and ASI (GDD6). The GLAI traits are $D_V$ (GDD6), $D_S$ (GDD6), $GLAI_M$ (m².m⁻²), $AUC_V$ (m².m⁻².GDD6), $AUC_F$ (m².m⁻².GDD6), and $AUC_S$ (m².m⁻².GDD6). Boxplots were constructed from the adjusted means</figcaption></figure>

of traits measured in 324 maize hybrids grown in the field in well-watered (WW) and water-deficient (WD) conditions. Boxes represent the inter-quartile range and whiskers extend to 1.5 times the inter-quartile range. Medians (horizontal lines inside the box), means (diamond symbols), and outliers (circle symbols) are indicated for each trait in each environment (color figure online)
17STPAULWD seems to have led to a delayed senescence ($D_S$) that could be explained by better light distribution in the canopy strata (Borrás et al. 2003; Huang et al. 2017; Yang et al. 2019). The drought effect on GLAI<sub>M</sub> and the consequences on $D_S$ compensated each other, so the decrease in AUC<sub>S</sub> in 2017 was only slight. Conversely, in 2016 the drought stress mainly occurred during grain filling, which curtailed the stay-green phase ($D_F$ and $D_S$, Fig. 3, Fig. 1), so AUC<sub>F</sub> and AUC<sub>S</sub> were smaller.

The relationships between the 11 traits are likely to vary under different conditions, so Pearson correlation coefficients ($r$) were calculated between pairs of traits analyzed within each environment (Fig. S3). The most strongly correlated traits were $D_V$ and $D_S$ ($r$ values range from 0.23 to -0.83), $D_V$ and AUC<sub>V</sub> ($r$ from 0.77 to 0.92), $D_S$ and AUC<sub>V</sub> ($r$ from -0.22 to -0.73), GLAI<sub>M</sub> and AUC<sub>V</sub> ($r$ from 0.84 to 0.90), GLAI<sub>M</sub> and AUC<sub>S</sub> ($r$ from 0.89 to 0.95), and GY and KN ($r$ from 0.64 to 0.86). There was therefore a strong compensation between the duration of the vegetative phase and the duration of the senescence phase, as described by the first axis of the PCA. Moreover, the AUC traits seemed more dependent on the amplitude of the dynamics (GLAI<sub>M</sub>) than on the duration of the phases. Correlations between yield and GLAI traits were strongest for AUC<sub>S</sub>, AUC<sub>F</sub> and GLAI<sub>M</sub> ($r$ from 0.23 to 0.50). Globally the latter three traits were more closely correlated to yield than were the flowering traits FF and ASI. TKW was never significantly correlated with GLAI traits, except in the 16STPAULWD environment where it was correlated with AUC<sub>F</sub> ($r$=0.16) and with the duration of senescence ($r$=0.26). The correlation between TKW and GY was also the strongest ($r$=0.33) in 16STPAULWD where the drought stress that occurred during grain filling had an impact on the duration of senescence.

<table><caption>**Table 1** Variance decomposition of six modeled GLAI traits and five agronomic traits in the maize trial network</caption><thead><tr><th>Trait</th><th>Unit</th><th>Min-Max</th><th>Mean</th><th>σ<sup>2</sup><sub>g</sub></th><th>CV<sub>g</sub></th><th>σ<sup>2</sup><sub>ge</sub></th><th>σ<sup>2</sup><sub>ge</sub>/σ<sup>2</sup><sub>g</sub>(%)</th><th>H<sup>2</sup></th></tr></thead><tbody><tr><td>D<sub>V</sub></td><td>GDD6</td><td>654-1005</td><td>829</td><td>4215</td><td>*** 7.83</td><td>3.91 10<sup>-8</sup></td><td>n.s</td><td>0.91</td></tr><tr><td>D<sub>S</sub></td><td>GDD6</td><td>620-892</td><td>790</td><td>475</td><td>*** 2.76</td><td>246</td><td>*** 51.86</td><td>0.85</td></tr><tr><td>GLAI<sub>M</sub></td><td>m<sup>2</sup>.m<sup>-2</sup></td><td>1.44-5.66</td><td>3.20</td><td>0.10</td><td>*** 9.64</td><td>2.4110<sup>-3</sup></td><td>*** 2.53</td><td>0.88</td></tr><tr><td>AUC<sub>V</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>362-1502</td><td>844</td><td>15,296</td><td>*** 14.65</td><td>188</td><td>*** 1.23</td><td>0.89</td></tr><tr><td>AUC<sub>F</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>512-2220</td><td>1287</td><td>23,764</td><td>*** 11.98</td><td>2545</td><td>*** 10.71</td><td>0.88</td></tr><tr><td>AUC<sub>S</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>681-2632</td><td>1589</td><td>17,779</td><td>*** 8.39</td><td>917</td><td>*** 5.16</td><td>0.87</td></tr><tr><td>GY</td><td>q.ha<sup>-1</sup></td><td>21-134</td><td>91</td><td>43.6</td><td>*** 7.25</td><td>22.5</td><td>*** 51.67</td><td>0.83</td></tr><tr><td>KN</td><td>kernels</td><td>954-4964</td><td>3446</td><td>69,258</td><td>*** 7.64</td><td>21,847</td><td>*** 31.54</td><td>0.82</td></tr><tr><td>TKW</td><td>g</td><td>140-367</td><td>265</td><td>240</td><td>*** 5.85</td><td>56.2</td><td>*** 23.40</td><td>0.85</td></tr><tr><td>FF</td><td>GDD6</td><td>823-1140</td><td>1007</td><td>514</td><td>*** 2.25</td><td>63.3</td><td>*** 12.32</td><td>0.85</td></tr><tr><td>ASI</td><td>GDD6</td><td>-87-187</td><td>31</td><td>51.2</td><td>*** 22.77</td><td>27.6</td><td>** 53.86</td><td>0.50</td></tr></tbody></table>

Min, Max and Mean were computed from the adjusted means of the traits of 324 hybrids over four environments. $\sigma_g^2$, genetic variance; $CV_g$, genetic coefficient of variation; $\sigma_{ge}^2$, G×E interaction variance; $\sigma_{ge}^2/\sigma_g^2$, ratio of G×E interaction variance to genetic variance (%); $H^2$, generalized heritability

## GLAI traits exhibit limited G×E interaction

The 11 traits analyzed all showed significant genetic effects (Table 1). The genetic variance in GY is noteworthy with $CV_g$=7%. All GLAI traits, except $D_S$, showed more genetic variance (higher $CV_g$) than GY, while ASI was the only agronomic traits showing more variance than GY. The two traits with the least genetic variance were $D_S$ and FF with $CV_g$ values around 2.5%. On the contrary, ASI and AUC<sub>V</sub> had the strongest genetic variances, with a $CV_g$ of 23% and 15% respectively. All agronomic traits had a high generalized heritability ($H^2$>0.8) except ASI, but GLAI traits were even more heritable. The G×E interaction was significant for all traits, except $D_V$. Apart from $D_S$, G×E interaction is smaller in GLAI traits (1–11% of the genetic variance $\sigma_g^2$) than in agronomic traits (12–54% of the genetic variance). The G×E interactions of $D_S$ and GY were however similar (52% of the genetic variance).

## Different associations between polymorphisms and traits in environments found with univariate and multivariate approaches

In a GWAS, a total of 7,806,995 SNPs were tested at 5% false discovery rate threshold for each environment-trait combination with two methods. With the univariate approach ($M_{UV}$) 49,459 significant associations (combinations between a trait, a SNP and an environment) were detected for all 11 analyzed traits. For all environments combined, FF represented 55% of the associations detected, GLAI<sub>M</sub> 12%, and AUC<sub>F</sub>, AUC<sub>S</sub>, GY and KN between 5 and 10% each (Fig. 4). No association was found for $D_V$, $D_S$ or ASI. Most other traits had at least

<footnote>***pvalue ≤ 0.001; **pvalue ≤ 0.01; *pvalue ≤ 0.05; n.s., pvalue > 0.05.</footnote>
17STPAULWD seems to have led to a delayed senescence (D<sub>S</sub>) that could be explained by better light distribution in the canopy strata (Borrás et al. 2003; Huang et al. 2017; Yang et al. 2019). The drought effect on GLAI<sub>M</sub> and the consequences on D<sub>S</sub> compensated each other, so the decrease in AUC<sub>S</sub> in 2017 was only slight. Conversely, in 2016 the drought stress mainly occurred during grain filling, which curtailed the stay-green phase (D<sub>F</sub> and D<sub>S</sub>, Fig. 3, Fig. 1), so AUC<sub>F</sub> and AUC<sub>S</sub> were smaller.

The relationships between the 11 traits are likely to vary under different conditions, so Pearson correlation coefficients (r) were calculated between pairs of traits analyzed within each environment (Fig. S3). The most strongly correlated traits were D<sub>V</sub> and D<sub>S</sub> (r values range from 0.23 to -0.83), D<sub>V</sub> and AUC<sub>V</sub> (r from 0.77 to 0.92), D<sub>S</sub> and AUC<sub>V</sub> (r from -0.22 to -0.73), GLAI<sub>M</sub> and AUC<sub>V</sub> (r from 0.84 to 0.90), GLAI<sub>M</sub> and AUC<sub>S</sub> (r from 0.89 to 0.95), and GY and KN (r from 0.64 to 0.86). There was therefore a strong compensation between the duration of the vegetative phase and the duration of the senescence phase, as described by the first axis of the PCA. Moreover, the AUC traits seemed more dependent on the amplitude of the dynamics (GLAI<sub>M</sub>) than on the duration of the phases. Correlations between yield and GLAI traits were strongest for AUC<sub>S</sub>, AUC<sub>F</sub> and GLAI<sub>M</sub> (r from 0.23 to 0.50). Globally the latter three traits were more closely correlated to yield than were the flowering traits FF and ASI. TKW was never significantly correlated with GLAI traits, except in the 16STPAULWD environment where it was correlated with AUC<sub>F</sub> (r = 0.16) and with the duration of senescence (r = 0.26). The correlation between TKW and GY was also the strongest (r = 0.33) in 16STPAULWD where the drought stress that occurred during grain filling had an impact on the duration of senescence.

<table><caption>**Table 1** Variance decomposition of six modeled GLAI traits and five agronomic traits in the maize trial network</caption><thead><tr><th>Trait</th><th>Unit</th><th>Min-Max</th><th>Mean</th><th>σ<sup>2</sup><sub>g</sub></th><th>CV<sub>g</sub></th><th>σ<sup>2</sup><sub>ge</sub></th><th>σ<sup>2</sup><sub>ge</sub>/σ<sup>2</sup><sub>g</sub>(%)</th><th>H<sup>2</sup></th></tr></thead><tbody><tr><td>D<sub>V</sub></td><td>GDD6</td><td>654-1005</td><td>829</td><td>4215</td><td>*** 7.83</td><td>3.91 10<sup>-8</sup></td><td>n.s</td><td>0.91</td></tr><tr><td>D<sub>S</sub></td><td>GDD6</td><td>620-892</td><td>790</td><td>475</td><td>*** 2.76</td><td>246</td><td>*** 51.86</td><td>0.85</td></tr><tr><td>GLAI<sub>M</sub></td><td>m<sup>2</sup>.m<sup>-2</sup></td><td>1.44-5.66</td><td>3.20</td><td>0.10</td><td>*** 9.64</td><td>2.4110<sup>-3</sup></td><td>*** 2.53</td><td>0.88</td></tr><tr><td>AUC<sub>V</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>362-1502</td><td>844</td><td>15,296</td><td>*** 14.65</td><td>188</td><td>*** 1.23</td><td>0.89</td></tr><tr><td>AUC<sub>F</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>512-2220</td><td>1287</td><td>23,764</td><td>*** 11.98</td><td>2545</td><td>*** 10.71</td><td>0.88</td></tr><tr><td>AUC<sub>S</sub></td><td>m<sup>2</sup>.m<sup>-2</sup>.GDD6</td><td>681-2632</td><td>1589</td><td>17,779</td><td>*** 8.39</td><td>917</td><td>*** 5.16</td><td>0.87</td></tr><tr><td>GY</td><td>q.ha<sup>-1</sup></td><td>21-134</td><td>91</td><td>43.6</td><td>*** 7.25</td><td>22.5</td><td>*** 51.67</td><td>0.83</td></tr><tr><td>KN</td><td>kernels</td><td>954-4964</td><td>3446</td><td>69,258</td><td>*** 7.64</td><td>21,847</td><td>*** 31.54</td><td>0.82</td></tr><tr><td>TKW</td><td>g</td><td>140-367</td><td>265</td><td>240</td><td>*** 5.85</td><td>56.2</td><td>*** 23.40</td><td>0.85</td></tr><tr><td>FF</td><td>GDD6</td><td>823-1140</td><td>1007</td><td>514</td><td>*** 2.25</td><td>63.3</td><td>*** 12.32</td><td>0.85</td></tr><tr><td>ASI</td><td>GDD6</td><td>-87-187</td><td>31</td><td>51.2</td><td>*** 22.77</td><td>27.6</td><td>** 53.86</td><td>0.50</td></tr></tbody></table>

Min, Max and Mean were computed from the adjusted means of the traits of 324 hybrids over four environments. σ<sup>2</sup><sub>g</sub>, genetic variance; CV<sub>g</sub>, genetic coefficient of variation; σ<sup>2</sup><sub>ge</sub>, G×E interaction variance; σ<sup>2</sup><sub>ge</sub>/σ<sup>2</sup><sub>g</sub>, ratio of G×E interaction variance to genetic variance (%); H<sup>2</sup>, generalized heritability

## GLAI traits exhibit limited G×E interaction

The 11 traits analyzed all showed significant genetic effects (Table 1). The genetic variance in GY is noteworthy with CV<sub>g</sub> = 7%. All GLAI traits, except D<sub>S</sub>, showed more genetic variance (higher CV<sub>g</sub>) than GY, while ASI was the only agronomic traits showing more variance than GY. The two traits with the least genetic variance were D<sub>S</sub> and FF with CV<sub>g</sub> values around 2.5%. On the contrary, ASI and AUC<sub>V</sub> had the strongest genetic variances, with a CV<sub>g</sub> of 23% and 15% respectively. All agronomic traits had a high generalized heritability (H<sup>2</sup> > 0.8) except ASI, but GLAI traits were even more heritable. The G×E interaction was significant for all traits, except D<sub>V</sub>. Apart from D<sub>S</sub>, G×E interaction is smaller in GLAI traits (1–11% of the genetic variance σ<sup>2</sup><sub>g</sub>) than in agronomic traits (12–54% of the genetic variance). The G×E interactions of D<sub>S</sub> and GY were however similar (52% of the genetic variance).

## Different associations between polymorphisms and traits in environments found with univariate and multivariate approaches

In a GWAS, a total of 7,806,995 SNPs were tested at 5% false discovery rate threshold for each environment-trait combination with two methods. With the univariate approach (M<sub>UV</sub>) 49,459 significant associations (combinations between a trait, a SNP and an environment) were detected for all 11 analyzed traits. For all environments combined, FF represented 55% of the associations detected, GLAI<sub>M</sub> 12%, and AUC<sub>F</sub>, AUC<sub>S</sub>, GY and KN between 5 and 10% each (Fig. 4). No association was found for D<sub>V</sub>, D<sub>S</sub> or ASI. Most other traits had at least

<footnote>***pvalue ≤ 0.001; **pvalue ≤ 0.01; *pvalue ≤ 0.05; n.s., pvalue > 0.05.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Number of associations detected in univariate ($M_{UV}$) and multivariate ($M_{MV}$) GWAS for 11 traits measured on 324 maize hybrids in four environments. The diameter of each circle is proportional to</figcaption></figure>

the total number of associations detected for each approach in that environment. Associations with KN, TKW and ASI were not tested in the $M_{MV}$ approach (color figure online)

one association in two different environments, but not GY, TKW and $AUC_V$. Flowering was the only trait with associations in all four environments.

With the multivariate approach ($M_{MV}$) 33,043 associated SNPs were revealed, but 9% of them could not be assigned to a trait with the a posteriori Wald test ($\alpha=0.01$). Thus, 113,164 significant associations were detected for eight traits in $M_{MV}$, which compares to 44,572 associations for the same eight traits for $M_{UV}$. All eight traits exhibited significant associations in each environment in $M_{MV}$, which contrasts sharply with the univariate approach (Fig. 4). The distribution of associations between each trait and environment was also much more even. For instance, FF, $GLAI_M$, $AUC_V$ and $AUC_S$ each represented 15–20% of associations.

It must be noted that in the univariate approach, most of the associated SNPs (86%) were associated with only one trait in a given environment, 11% with two, 3% with three and 0.4% with four environments. In contrast, the multivariate approach detected many more SNPs associated with several traits in a given environment. Specifically, 79% of SNPs were associated with 2–8 traits, and 23%, 22% and 17% of SNPs were associated with 4, 5 and 6 traits, respectively.

Most associations were specific to the GWAS model used, with 54% specific associations in $M_{UV}$ and 85% in $M_{MV}$. Only 13% of GLAI trait associations were detected by both models (Fig. S5), 9% of GY associations and 19% FF associations.

## GLAI traits have a simpler genetic architecture than yield

Significant SNPs were clustered for each trait-environment combination based on linkage disequilibrium. Within mixed clusters (i.e., defined as a mixture of SNPs detected by $M_{UV}$ and $M_{MV}$), most associated SNPs in each metacluster were from the multivariate approach (81% of cases).

Among the 1444 clusters identified, 19 were longer than 10% of the corresponding chromosome. The lengths of the clusters were quite short compared to the large number of associated SNPs, with 75% of clusters shorter than 0.25 Mbp (90% < 3.5 Mbp).

As expected, for a given trait, the number of clusters identified depends strongly on the number of associated SNPs (Table S2). However, the 25,591 SNPs associated with FF in 17STPAULWD were grouped into just 66 clusters. On the contrary, the 2,079 SNPs associated with GY in the same environment were grouped into as many as 96 clusters, the maximum detected in this study for a single trait in a given environment. This is consistent with the highly polygenic architecture of yield.

A backward elimination approach was applied to the large number of genomic regions detected. This had the twofold advantage of limiting the risk of false positives and reducing the number of regions to analyze. Thus, only a third of the clusters were kept with 1 to 20 QTLs depending on the trait and environment considered, thus accounting for 475 QTL in total (Tables S2 and S3). Moreover, the multi-environment
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Number of associations detected in univariate ($M_{UV}$) and multivariate ($M_{MV}$) GWAS for 11 traits measured on 324 maize hybrids in four environments. The diameter of each circle is proportional to</figcaption></figure>

the total number of associations detected for each approach in that environment. Associations with KN, TKW and ASI were not tested in the $M_{MV}$ approach (color figure online)

one association in two different environments, but not GY, TKW and $AUC_V$. Flowering was the only trait with associations in all four environments.

With the multivariate approach ($M_{MV}$) 33,043 associated SNPs were revealed, but 9% of them could not be assigned to a trait with the a posteriori Wald test ($\alpha=0.01$). Thus, 113,164 significant associations were detected for eight traits in $M_{MV}$, which compares to 44,572 associations for the same eight traits for $M_{UV}$. All eight traits exhibited significant associations in each environment in $M_{MV}$, which contrasts sharply with the univariate approach (Fig. 4). The distribution of associations between each trait and environment was also much more even. For instance, FF, $GLAI_M$, $AUC_V$ and $AUC_S$ each represented 15–20% of associations.

It must be noted that in the univariate approach, most of the associated SNPs (86%) were associated with only one trait in a given environment, 11% with two, 3% with three and 0.4% with four environments. In contrast, the multivariate approach detected many more SNPs associated with several traits in a given environment. Specifically, 79% of SNPs were associated with 2–8 traits, and 23%, 22% and 17% of SNPs were associated with 4, 5 and 6 traits, respectively.

Most associations were specific to the GWAS model used, with 54% specific associations in $M_{UV}$ and 85% in $M_{MV}$. Only 13% of GLAI trait associations were detected by both models (Fig. S5), 9% of GY associations and 19% FF associations.

## GLAI traits have a simpler genetic architecture than yield

Significant SNPs were clustered for each trait-environment combination based on linkage disequilibrium. Within mixed clusters (i.e., defined as a mixture of SNPs detected by $M_{UV}$ and $M_{MV}$), most associated SNPs in each metacluster were from the multivariate approach (81% of cases).

Among the 1444 clusters identified, 19 were longer than 10% of the corresponding chromosome. The lengths of the clusters were quite short compared to the large number of associated SNPs, with 75% of clusters shorter than 0.25 Mbp (90% < 3.5 Mbp).

As expected, for a given trait, the number of clusters identified depends strongly on the number of associated SNPs (Table S2). However, the 25,591 SNPs associated with FF in 17STPAULWD were grouped into just 66 clusters. On the contrary, the 2,079 SNPs associated with GY in the same environment were grouped into as many as 96 clusters, the maximum detected in this study for a single trait in a given environment. This is consistent with the highly polygenic architecture of yield.

A backward elimination approach was applied to the large number of genomic regions detected. This had the twofold advantage of limiting the risk of false positives and reducing the number of regions to analyze. Thus, only a third of the clusters were kept with 1 to 20 QTLs depending on the trait and environment considered, thus accounting for 475 QTL in total (Tables S2 and S3). Moreover, the multi-environment
backward selection model can also be used to identify QTLs in environments where they had not been initially detected (Table S2). While most QTLs are concentrated on chromosomes 1, 3, 4, 5 and 8, there are at least 15 QTLs on each chromosome.

In each environment, the average proportion of phenotypic variance explained by a QTL is 4%, with a maximum of 11%, and no major differences between GLAI traits and agronomic traits. Nevertheless, among the 16 QTLs explaining less than 1% of the phenotypic variance, eight were yield QTLs. Some differences are also apparent at the network scale: the mean and maximum proportion of variance explained by QTL were, respectively, 3.5% and 9.5% for GLAI traits, 2.3% and 5.2% for yield and its components, and 2.8% and 6.1% for flowering. This confirms that GLAI traits have a weaker G×E interaction than yield and flowering.

Overall, the QTLs detected captured 13–61% of the phenotypic variance depending on the trait considered (Fig. 5, $M_{ALL}$). To evaluate the respective contributions of $M_{UV}$ and $M_{MV}$ approaches to explain the phenotypic variance, two new sets of QTLs were selected by backward elimination, from the SNPs detected independently by each approach. It appeared that the QTLs detected by $M_{UV}$ explained less of the variance in all the traits analyzed, except FF. Indeed, with the $M_{MV}$ approach $r^2$ values for the six GLAI traits and GY were 4–43% higher. The combined use of $M_{UV}$ and $M_{MV}$

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Proportion of phenotypic variance explained by all the QTLs detected on 324 maize hybrids evaluated in four environments. For each trait, we compared the variance explained by the QTLs from a univariate GWAS approach ($M_{UV}$), the QTLs from a multivariate approach ($M_{MV}$), and the QTLs from $M_{ALL}$ built from the combination of SNPs detected by $M_{UV}$ and $M_{MV}$. Associations with KN, TKW were not tested with $M_{MV}$. No association was found for ASI</figcaption></figure>

brought little gain compared to $M_{MV}$ alone, at most +11% for $AUC_S$.

## GLAI QTLs show many colocalizations

According to the genomic sequence, many QTLs are colocalizing suggesting significant pleiotropy or genetic linkage between the analyzed traits, in particular the GLAI traits. Indeed, 53 to 89% of GLAI QTLs colocalized with QTLs for at least one other trait in one or more environments, compared to only 27–43% of agronomic QTLs. Moreover, QTLs of the same trait in different environments often colocalized, implying that trait- and environment-specific QTLs are rare (Fig. S6). It is also noteworthy that most of the QTLs had similar effects in the four environments and it was rare that the direction of the effect reversed (Figs. S4 and S7).

Overall, GLAI traits shared a significant proportion of QTLs with same-direction effects, especially $AUC_V$, $AUC_F$ and $AUC_S$, but also $AUC_V$ and $D_V$ (Fig. S6). The colocalizations between GLAI traits and GY or FF, distributed over one to seven genomic regions, were rarer. Despite this marked trend, there were some QTL colocalizations with opposite effects for almost all trait pairs.

Yield QTLs detected under drought stress showed fewer colocalizations than those detected under optimal conditions, especially for 16STPAULWD where the stress lasted longer (Fig. S6). However, interestingly, in both stressed environments, some yield QTLs colocalized with GLAI QTLs detected under optimal conditions. Additionally, yield and flowering QTLs rarely colocalized, with only three common regions, each affecting both traits in the same direction. In particular, there was a single colocalization between yield (17STPAULWD) and flowering (16STPAULWD) under stressed conditions.

Some traits showed only colocalizations with opposite effects, such as $D_V$ and $D_S$, $AUC_V$ and $D_S$, and $GLAI_M$ and $D_S$ (Fig. S6), which is consistent with the mainly negative phenotypic (Fig. S3) and genetic (Fig. S8) correlations between these pairs of traits. Notably, three colocalizations with opposite effects were detected for $AUC_S$ and GY versus only one with effects in the same direction. This trend is surprising because it contradicts the marked positive correlations between $AUC_S$ and GY at the phenotypic and genetic level and probably reflects the weaknesses of the QTL colocalization method, or the backward elimination approach for the selection of the final QTL list.

## GLAI QTLs can explain yield variability under drought stress conditions

In the $GY_{11}$ network, maize GY measured over 11 WD trials spread over five locations in three countries and four years varied on average from 53 to 89 q.ha⁻¹ with
backward selection model can also be used to identify QTLs in environments where they had not been initially detected (Table S2). While most QTLs are concentrated on chromosomes 1, 3, 4, 5 and 8, there are at least 15 QTLs on each chromosome.

In each environment, the average proportion of phenotypic variance explained by a QTL is 4%, with a maximum of 11%, and no major differences between GLAI traits and agronomic traits. Nevertheless, among the 16 QTLs explaining less than 1% of the phenotypic variance, eight were yield QTLs. Some differences are also apparent at the network scale: the mean and maximum proportion of variance explained by QTL were, respectively, 3.5% and 9.5% for GLAI traits, 2.3% and 5.2% for yield and its components, and 2.8% and 6.1% for flowering. This confirms that GLAI traits have a weaker G×E interaction than yield and flowering.

Overall, the QTLs detected captured 13–61% of the phenotypic variance depending on the trait considered (Fig. 5, $M_{ALL}$). To evaluate the respective contributions of $M_{UV}$ and $M_{MV}$ approaches to explain the phenotypic variance, two new sets of QTLs were selected by backward elimination, from the SNPs detected independently by each approach. It appeared that the QTLs detected by $M_{UV}$ explained less of the variance in all the traits analyzed, except FF. Indeed, with the $M_{MV}$ approach $r^2$ values for the six GLAI traits and GY were 4–43% higher. The combined use of $M_{UV}$ and $M_{MV}$

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Proportion of phenotypic variance explained by all the QTLs detected on 324 maize hybrids evaluated in four environments. For each trait, we compared the variance explained by the QTLs from a univariate GWAS approach ($M_{UV}$), the QTLs from a multivariate approach ($M_{MV}$), and the QTLs from $M_{ALL}$ built from the combination of SNPs detected by $M_{UV}$ and $M_{MV}$. Associations with KN, TKW were not tested with $M_{MV}$. No association was found for ASI</figcaption></figure>

brought little gain compared to $M_{MV}$ alone, at most +11% for $AUC_S$.

## GLAI QTLs show many colocalizations

According to the genomic sequence, many QTLs are colocalizing suggesting significant pleiotropy or genetic linkage between the analyzed traits, in particular the GLAI traits. Indeed, 53 to 89% of GLAI QTLs colocalized with QTLs for at least one other trait in one or more environments, compared to only 27–43% of agronomic QTLs. Moreover, QTLs of the same trait in different environments often colocalized, implying that trait- and environment-specific QTLs are rare (Fig. S6). It is also noteworthy that most of the QTLs had similar effects in the four environments and it was rare that the direction of the effect reversed (Figs. S4 and S7).

Overall, GLAI traits shared a significant proportion of QTLs with same-direction effects, especially $AUC_V$, $AUC_F$ and $AUC_S$, but also $AUC_V$ and $D_V$ (Fig. S6). The colocalizations between GLAI traits and GY or FF, distributed over one to seven genomic regions, were rarer. Despite this marked trend, there were some QTL colocalizations with opposite effects for almost all trait pairs.

Yield QTLs detected under drought stress showed fewer colocalizations than those detected under optimal conditions, especially for 16STPAULWD where the stress lasted longer (Fig. S6). However, interestingly, in both stressed environments, some yield QTLs colocalized with GLAI QTLs detected under optimal conditions. Additionally, yield and flowering QTLs rarely colocalized, with only three common regions, each affecting both traits in the same direction. In particular, there was a single colocalization between yield (17STPAULWD) and flowering (16STPAULWD) under stressed conditions.

Some traits showed only colocalizations with opposite effects, such as $D_V$ and $D_S$, $AUC_V$ and $D_S$, and $GLAI_M$ and $D_S$ (Fig. S6), which is consistent with the mainly negative phenotypic (Fig. S3) and genetic (Fig. S8) correlations between these pairs of traits. Notably, three colocalizations with opposite effects were detected for $AUC_S$ and GY versus only one with effects in the same direction. This trend is surprising because it contradicts the marked positive correlations between $AUC_S$ and GY at the phenotypic and genetic level and probably reflects the weaknesses of the QTL colocalization method, or the backward elimination approach for the selection of the final QTL list.

## GLAI QTLs can explain yield variability under drought stress conditions

In the $GY_{11}$ network, maize GY measured over 11 WD trials spread over five locations in three countries and four years varied on average from 53 to 89 q.ha⁻¹ with
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Grain yield variability explained by GLAI QTLs in the GY<sub>11</sub> drought stressed network. **A** Effect of a wide range of drought stress on yield of 324 hybrids measured in eleven environments. The black boxplot corresponds to the whole GY<sub>11</sub> network, which is composed of the eleven stressed environments detailed below. The gray boxplots correspond to the two WD trials used to identify the GLAI QTLs then used in the backward selection approach to explain grain yield variability in the whole network. Boxplots were constructed from the adjusted means. The vertical line in the boxplots corresponds to the</figcaption></figure>

median, while the diamond corresponds to the mean. The box length represents the inter-quartile range, and the whiskers extend to 1.5 times the inter-quartile range. Generalized heritability is given at the left of each boxplot. **B** Proportion of grain yield variance explained by GLAI QTLs. Proportion of explained variance, computed as r², is given in percent. Color use is the same as panel A. The number of GLAI QTLs with a significant effect on grain yield in the network, or each environment, is given at the right of each bar

a generalized heritability between 0.63 and 0.80, which underlines the diversity of drought stress that occurred during these trials (Fig. 6A). The use of a multi-environment backward selection model from all GLAI QTLs (120 unique regions for D<sub>V</sub>, D<sub>S</sub>, AUC<sub>V</sub>, AUC<sub>F</sub>, AUC<sub>S</sub>, GLAI<sub>M</sub>) led to the identification of 16 GLAI QTLs distributed over seven chromosomes having a significant effect on GY in this drought stressed network (Fig. 6B, Fig. S9). Taken together, the 16 QTLs explained almost one-fifth (r<sup>2</sup>=0.18, *p*value < 10<sup>-3</sup>, Fig. 6B) of the phenotypic variability observed over the 11 environments, and from 7 to 18% of the phenotypic variability observed in each environment, with 3–8 significant QTLs. Interestingly, each GLAI trait was associated with at least one of the 16 QTLs selected by the backward selection model, confirming the relevance of the impact of the six GLAI traits on yield.

## Discussion

To quantify the changes in GLAI in diverse maize genotypes, a panel of 324 DH lines from a MAGIC population was phenotyped by UAV spectral imaging in four trials under optimal conditions and drought stress. Six traits derived from the GLAI dynamics were shown to capture a significant proportion of the phenotypic diversity of the panel. To identify the genetic basis for these secondary traits and their influence on yield and drought tolerance, two GWAS approaches were compared: a univariate model and a multivariate model.

## Traits governing GLAI dynamics are less complex than yield

The traits evaluated within the trial network showed a high heritability, except for ASI. The low heritability of ASI (H<sup>2</sup>=0.5) is consistent with the results of Cairns et al. (2013), Trachsel et al. (2016) and Bouchet et al. (2017) but much lower than the heritability reported for American and Chinese nested association mapping populations (Li et al. 2016c). The heritability of FF was also similar to previous estimates (Li et al. 2016c; Bouchet et al. 2017). Yield and its components show strong heritability compared to previous studies (Trachsel et al. 2016; Li et al. 2016a; Bouchet et al. 2017; Cerrudo et al. 2018), as do the AUC (Christopher et al. 2014; Trachsel et al. 2016; Cerrudo et al. 2018) and senescence-related traits (Messmer et al. 2011; Ziyomo and Bernardo 2013; Almeida et al. 2014; Yang et al. 2017). The new traits describing vegetative growth of the canopy (D<sub>V</sub>, GLAI<sub>M</sub>, AUC<sub>V</sub>) showed similar heritability to traits such as leaf number and leaf size (Bouchet et al. 2017; Pan et al. 2017; Wang et al. 2017; Zhao et al. 2019).

GLAI traits were more heritable than agronomic traits, including yield and its components, and often had higher genetic variance. Moreover, although heritability gap between GY and GLAI traits is modest when computed at the network level, it can be noted that differences increase when heritability is evaluated at the trial level with a mean (maximum) gain of heritability for GLAI traits compared to GY of 0.10 (0.15) and 0.21 (0.3) for 16STPAULWW and 16STPAULWD, respectively, and 0 (0.12), 0.17 (0.28) for 17STPAULWW and 17STPAULWD, respectively. In
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Grain yield variability explained by GLAI QTLs in the GY<sub>11</sub> drought stressed network. **A** Effect of a wide range of drought stress on yield of 324 hybrids measured in eleven environments. The black boxplot corresponds to the whole GY<sub>11</sub> network, which is composed of the eleven stressed environments detailed below. The gray boxplots correspond to the two WD trials used to identify the GLAI QTLs then used in the backward selection approach to explain grain yield variability in the whole network. Boxplots were constructed from the adjusted means. The vertical line in the boxplots corresponds to the</figcaption></figure>

median, while the diamond corresponds to the mean. The box length represents the inter-quartile range, and the whiskers extend to 1.5 times the inter-quartile range. Generalized heritability is given at the left of each boxplot. **B** Proportion of grain yield variance explained by GLAI QTLs. Proportion of explained variance, computed as r², is given in percent. Color use is the same as panel A. The number of GLAI QTLs with a significant effect on grain yield in the network, or each environment, is given at the right of each bar

a generalized heritability between 0.63 and 0.80, which underlines the diversity of drought stress that occurred during these trials (Fig. 6A). The use of a multi-environment backward selection model from all GLAI QTLs (120 unique regions for D<sub>V</sub>, D<sub>S</sub>, AUC<sub>V</sub>, AUC<sub>F</sub>, AUC<sub>S</sub>, GLAI<sub>M</sub>) led to the identification of 16 GLAI QTLs distributed over seven chromosomes having a significant effect on GY in this drought stressed network (Fig. 6B, Fig. S9). Taken together, the 16 QTLs explained almost one-fifth (r<sup>2</sup>=0.18, *p*value < 10<sup>-3</sup>, Fig. 6B) of the phenotypic variability observed over the 11 environments, and from 7 to 18% of the phenotypic variability observed in each environment, with 3–8 significant QTLs. Interestingly, each GLAI trait was associated with at least one of the 16 QTLs selected by the backward selection model, confirming the relevance of the impact of the six GLAI traits on yield.

## Discussion

To quantify the changes in GLAI in diverse maize genotypes, a panel of 324 DH lines from a MAGIC population was phenotyped by UAV spectral imaging in four trials under optimal conditions and drought stress. Six traits derived from the GLAI dynamics were shown to capture a significant proportion of the phenotypic diversity of the panel. To identify the genetic basis for these secondary traits and their influence on yield and drought tolerance, two GWAS approaches were compared: a univariate model and a multivariate model.

## Traits governing GLAI dynamics are less complex than yield

The traits evaluated within the trial network showed a high heritability, except for ASI. The low heritability of ASI (H<sup>2</sup>=0.5) is consistent with the results of Cairns et al. (2013), Trachsel et al. (2016) and Bouchet et al. (2017) but much lower than the heritability reported for American and Chinese nested association mapping populations (Li et al. 2016c). The heritability of FF was also similar to previous estimates (Li et al. 2016c; Bouchet et al. 2017). Yield and its components show strong heritability compared to previous studies (Trachsel et al. 2016; Li et al. 2016a; Bouchet et al. 2017; Cerrudo et al. 2018), as do the AUC (Christopher et al. 2014; Trachsel et al. 2016; Cerrudo et al. 2018) and senescence-related traits (Messmer et al. 2011; Ziyomo and Bernardo 2013; Almeida et al. 2014; Yang et al. 2017). The new traits describing vegetative growth of the canopy (D<sub>V</sub>, GLAI<sub>M</sub>, AUC<sub>V</sub>) showed similar heritability to traits such as leaf number and leaf size (Bouchet et al. 2017; Pan et al. 2017; Wang et al. 2017; Zhao et al. 2019).

GLAI traits were more heritable than agronomic traits, including yield and its components, and often had higher genetic variance. Moreover, although heritability gap between GY and GLAI traits is modest when computed at the network level, it can be noted that differences increase when heritability is evaluated at the trial level with a mean (maximum) gain of heritability for GLAI traits compared to GY of 0.10 (0.15) and 0.21 (0.3) for 16STPAULWW and 16STPAULWD, respectively, and 0 (0.12), 0.17 (0.28) for 17STPAULWW and 17STPAULWD, respectively. In
addition, the G×E interaction was lower for most GLAI traits than for agronomic traits (Table 1, Table S4). Bouchet et al. (2017) and Yang et al. (2014) made similar observations regarding architectural canopy traits (number of leaves, size of leaves) and yield components, with comparable values of the ratio $\sigma_{ge}^2/\sigma_g^2$. In our study, the heritability was slightly higher when G×E interaction ratios were lower, whereas this relationship was more pronounced in the results of Yang et al. (2014) but not evidenced in results from Bouchet et al. (2017). These differences confirm that beyond the characteristics specific to the traits analyzed, the conclusions drawn from such genetic analyses depend on the population considered and the experimental network. For this study, we used a panel produced from 16 lines representative of the genetic diversity used for hybrid production in temperate zones, encompassing a very broad diversity of canopy structures (number, size, color and angle of the leaves). While the test-cross was necessary for the results to be pertinent to elite or breeding material as correlations are often weak between lines and hybrids (Cairns et al. 2012; Trachsel et al. 2016), phenotypic evaluation after the test cross clearly reduced the panel diversity. The limited size of this experimental network probably also did not reveal all the existing variability in our panel for the considered traits. Our results suggest that the traits extracted from GLAI dynamics are less complex than yield-related traits.

## Larger green leaf area increases grain yield under optimal conditions and moderate drought

GLAI traits were highly correlated with each other at the phenotypic and genetic levels (Figs. S3 and S8). They were also significantly correlated with yield and flowering. However, yield was more correlated with GLAI than with flowering, which is consistent with the results of Trachsel et al. (2016) and Cerrudo et al. (2018). It is interesting to note that the duration of the vegetative phase ($D_v$) and the duration of the senescence phase ($D_s$) were negatively correlated, and that this opposition was the main factor differentiating the genotypes of the panel (Fig. 2). This relationship shows that it is possible to lengthen the senescence phase by shortening the vegetative phase while maintaining the same total duration. Thus, the duration of grain filling could possibly be increased without altering the maturity date of a genotype. Trachsel et al. (2017) obtained similar results showing under optimal and stressed conditions that, provided vegetative development was vigorous, earlier flowering extended the grain filling duration, thus increasing grain yield. In our study, the traits most consistently and strongly correlated overall with yield were those that reflect vigorous vegetative development and hence a large green leaf area during grain filling: $AUC_v$, $AUC_F$, $GLAI_M$ and $AUC_S$ with $0.14 \le r_p \le 0.50$ and $-0.04 \le r_g \le 0.94$. These moderate

correlations were expected as the proportion of yield variation explained by individual traits is known to be low at intermediate levels of water-deficit and stress, like the ones experienced in this study (Cooper and Messina 2023).

Stay-green, defined as the maintenance of photosynthetic activity during filling, is one of the major determinants of yield under optimal and stressed conditions (Wolfe et al. 1988; Cairns et al. 2012; Almeida et al. 2014; Kante et al. 2016) and has played a major role in the improvement of new hybrids in recent decades (Duvick 2005). Stay-green may be the consequence of increased green surface area before flowering, postponement of senescence, or slower senescence (Thomas and Howarth 2000; Christopher et al. 2014, 2016). $AUC_S$ is a good indicator of stay-green because it incorporates each of these aspects. Indeed, $AUC_S$ was the GLAI trait most correlated with yield in our four environments (Figs. S3 and S8) which confirms the major interest of this trait. While numerous studies stated that stay-green can originate from 'shifting water' from pre- to post-anthesis, whether by optimal transpiration efficiency, root uptake or reduced transpiration area, that would lengthen senescence duration (Borrell et al. 2000, 2014), it is interesting to note that in this study, stay-green was the result of more green surface grown before flowering and an acceleration of senescence during the rapid phase as seen by the strong correlations between $AUC_S$, $GLAI_M$ and $S_{RS}$. The weak correlation of $AUC_S$ and $D_S$ shows that slower senescence was not a contributory factor. Moreover, $D_S$ was only weakly correlated with GY. A well-developed foliar system was therefore central to the development of yield within our experimental network whether in optimal conditions or moderate drought stress with full irrigation at the end of the cycle. Maintaining leaf growth and delaying senescence increase transpiration, but if water reserves are restored before the end of the cycle, gas exchange and photosynthesis can still be maximized (Tardieu 2012). If not, with a severe terminal stress, maintaining leaf growth would be deleterious. In such conditions, typical to the Australian environment type experienced by Borrell et al. (2000, 2014), it would be expected that GY would be negatively correlated with traits that increase leaf transpiration, such as areas under the curve or $GLAI_M$, while being positively correlated with $D_S$.

These specificities highlight that secondary traits are only of interest if associated with well-characterized types of environments, which refers to ideotyping analysis (Hammer et al. 2014; Bustos-Korts et al. 2019b; Cooper and Messina 2023). Indeed, correlation between secondary traits and GY varies between environment type, and even across the crop cycle in a given environment. With a good understanding of environment type, one can expect identifying the most relevant secondary traits that maximize correlation with GY to build one or more associated ideotype, if there is sufficient knowledge about environment variability (Cooper and Messina 2023).
addition, the G×E interaction was lower for most GLAI traits than for agronomic traits (Table 1, Table S4). Bouchet et al. (2017) and Yang et al. (2014) made similar observations regarding architectural canopy traits (number of leaves, size of leaves) and yield components, with comparable values of the ratio $\sigma_{ge}^2/\sigma_g^2$. In our study, the heritability was slightly higher when G×E interaction ratios were lower, whereas this relationship was more pronounced in the results of Yang et al. (2014) but not evidenced in results from Bouchet et al. (2017). These differences confirm that beyond the characteristics specific to the traits analyzed, the conclusions drawn from such genetic analyses depend on the population considered and the experimental network. For this study, we used a panel produced from 16 lines representative of the genetic diversity used for hybrid production in temperate zones, encompassing a very broad diversity of canopy structures (number, size, color and angle of the leaves). While the test-cross was necessary for the results to be pertinent to elite or breeding material as correlations are often weak between lines and hybrids (Cairns et al. 2012; Trachsel et al. 2016), phenotypic evaluation after the test cross clearly reduced the panel diversity. The limited size of this experimental network probably also did not reveal all the existing variability in our panel for the considered traits. Our results suggest that the traits extracted from GLAI dynamics are less complex than yield-related traits.

## Larger green leaf area increases grain yield under optimal conditions and moderate drought

GLAI traits were highly correlated with each other at the phenotypic and genetic levels (Figs. S3 and S8). They were also significantly correlated with yield and flowering. However, yield was more correlated with GLAI than with flowering, which is consistent with the results of Trachsel et al. (2016) and Cerrudo et al. (2018). It is interesting to note that the duration of the vegetative phase ($D_v$) and the duration of the senescence phase ($D_s$) were negatively correlated, and that this opposition was the main factor differentiating the genotypes of the panel (Fig. 2). This relationship shows that it is possible to lengthen the senescence phase by shortening the vegetative phase while maintaining the same total duration. Thus, the duration of grain filling could possibly be increased without altering the maturity date of a genotype. Trachsel et al. (2017) obtained similar results showing under optimal and stressed conditions that, provided vegetative development was vigorous, earlier flowering extended the grain filling duration, thus increasing grain yield. In our study, the traits most consistently and strongly correlated overall with yield were those that reflect vigorous vegetative development and hence a large green leaf area during grain filling: $AUC_v$, $AUC_F$, $GLAI_M$ and $AUC_S$ with $0.14 \le r_p \le 0.50$ and $-0.04 \le r_g \le 0.94$. These moderate

correlations were expected as the proportion of yield variation explained by individual traits is known to be low at intermediate levels of water-deficit and stress, like the ones experienced in this study (Cooper and Messina 2023).

Stay-green, defined as the maintenance of photosynthetic activity during filling, is one of the major determinants of yield under optimal and stressed conditions (Wolfe et al. 1988; Cairns et al. 2012; Almeida et al. 2014; Kante et al. 2016) and has played a major role in the improvement of new hybrids in recent decades (Duvick 2005). Stay-green may be the consequence of increased green surface area before flowering, postponement of senescence, or slower senescence (Thomas and Howarth 2000; Christopher et al. 2014, 2016). $AUC_S$ is a good indicator of stay-green because it incorporates each of these aspects. Indeed, $AUC_S$ was the GLAI trait most correlated with yield in our four environments (Figs. S3 and S8) which confirms the major interest of this trait. While numerous studies stated that stay-green can originate from 'shifting water' from pre- to post-anthesis, whether by optimal transpiration efficiency, root uptake or reduced transpiration area, that would lengthen senescence duration (Borrell et al. 2000, 2014), it is interesting to note that in this study, stay-green was the result of more green surface grown before flowering and an acceleration of senescence during the rapid phase as seen by the strong correlations between $AUC_S$, $GLAI_M$ and $S_{RS}$. The weak correlation of $AUC_S$ and $D_S$ shows that slower senescence was not a contributory factor. Moreover, $D_S$ was only weakly correlated with GY. A well-developed foliar system was therefore central to the development of yield within our experimental network whether in optimal conditions or moderate drought stress with full irrigation at the end of the cycle. Maintaining leaf growth and delaying senescence increase transpiration, but if water reserves are restored before the end of the cycle, gas exchange and photosynthesis can still be maximized (Tardieu 2012). If not, with a severe terminal stress, maintaining leaf growth would be deleterious. In such conditions, typical to the Australian environment type experienced by Borrell et al. (2000, 2014), it would be expected that GY would be negatively correlated with traits that increase leaf transpiration, such as areas under the curve or $GLAI_M$, while being positively correlated with $D_S$.

These specificities highlight that secondary traits are only of interest if associated with well-characterized types of environments, which refers to ideotyping analysis (Hammer et al. 2014; Bustos-Korts et al. 2019b; Cooper and Messina 2023). Indeed, correlation between secondary traits and GY varies between environment type, and even across the crop cycle in a given environment. With a good understanding of environment type, one can expect identifying the most relevant secondary traits that maximize correlation with GY to build one or more associated ideotype, if there is sufficient knowledge about environment variability (Cooper and Messina 2023).
This will thus allow a better breeding efficiency through indi-
rect selection.

Multivariate GWAS is the most powerful approach

With the combination of very dense genotyping data
(7,806,995 SNPs) and two GWAS approaches, many associa-
tions (167,854) were detected for 11 traits measured in four
environments. The multivariate approach (M<sub>MV</sub>) detected
many more associations, including most of the associations
revealed by the univariate approach (M<sub>UV</sub>, Fig. S5). Moreo-
ver, M<sub>MV</sub> revealed many more associations common to several
traits (79% of associations) than M<sub>UV</sub> and explained a greater
part of the phenotypic variance, particularly for yield (Fig. 5).

The power gained in GWAS with the multivariate approach compared to the univariate approach has been demonstrated by several studies (Zhu and Zhang 2009; O'Reilly et al. 2012; Stephens 2013; Galesloot et al. 2014; Porter and O'Reilly 2017). Taking into account the covariances among traits brings additional information that increases the power in detecting genetic variants involved in the control of several traits, but also variants affecting a single trait (Ferreira and Purcell 2009; Stephens 2013; Carlson et al. 2019). However, the use of multivariate models may be limited by the number of phenotypes that can be integrated. For example, in this study, a maximum of eight traits could be integrated into the multivariate model fitted with the GEMMA software. In recent years, many tools have been developed for multi-trait GWAS that differ in terms of the number of analyzable traits, calculation times and assumptions on the distribution of phenotypes (O'Reilly et al. 2012; Sluis et al. 2013; Stephens 2013; Zhou and Stephens 2014). These tools also have a very variable power depending on the phenotypic correlations and the genetic architecture of the traits studied (Porter and O'Reilly 2017).

In this work, we used a two-step approach for its com-
putational efficiency at the cost of a probable loss of infor-
mation between the two steps. The gain in power noted in
several studies using functional mapping, and especially
random regression (Wu and Lin 2006; Li et al. 2010; Ning
et al. 2017, 2018; Campbell et al. 2019; Moreira et al. 2020)
may have better enhanced our understanding of the timing
of the genetic control of GLAI dynamics in response to
drought stress. However, computationally it would not have
been practicable with our dense genotyping data and the
backward-elimination method (Ning et al. 2017; Moreira
et al. 2020).

GLAI and agronomic traits exhibit polygenic
architecture and common genetic determinants

The associations we detected overlap strongly with the
genetic determinisms of GLAI traits, yield and flowering

(Fig. S6). Our results, like those of Cheverud (1988), Bou-
chet et al. (2017), Pan et al. (2017), Sodini et al. (2018)
and Yuan et al. (2019) show that correlations at the phe-
notypic level are good indicators of the sharing of genetic
determinants (pleiotropy or genetic linkage) between two
traits (Fig. S10A). This indicates that it is possible to select
secondary yield traits quite efficiently based on phenotypic
correlations. However, we observed only a moderate link
between phenotypic and genetic correlations (Fig. S10B),
which demonstrates the complementarity of these two pieces
of information. This partial discrepancy is probably due to
permanent environmental effects at the plot level throughout
the cycle, captured by the temporal approach implemented in
this study (Hadfield et al. 2007; Canela-Xandri et al. 2018).

Common genetic determinisms between earliness, plant
architecture and changes in leaf area during the cycle are
well known (Li et al. 2016b; Bouchet et al. 2017; Pan et al.
2017; Condorelli et al. 2018; Spindel et al. 2018). The entire
foliar and reproductive apparatus can be traced back to the
activity of the apical meristem and the floral meristems it
initiates (Kwiatkowska 2008), which is probably the reason
for these common determinisms (Thompson et al. 2015;
Baute et al. 2015). The impact of changes in leaf area on
light interception, gas exchange and water fluxes explains
the large number of QTLs affecting both yield and leaf area
dynamics under optimal or stressed conditions (Trachsel
et al. 2016; Cerrudo et al. 2018; Condorelli et al. 2018).
In our study, however, there is less colocalization between
yield and flowering than between yield and GLAI traits,
which could be the consequence of purposely reducing the
flowering time range of the panel. Having limited the flower-
ing date confounding effect, it can be assumed that a large
proportion of the common QTLs between yield and GLAI
traits that do not colocalize with flowering are due to direct
QTL effects on the GLAI. Moreover, we showed that some
GLAI QTLs detected under optimal conditions colocalize
with yield QTLs detected under stressed conditions. This
means yield under drought stress is partly due to GLAI
QTLs detected in the absence of stress, which may facilitate
the prediction of grain yield in breeding.

Individually, the detected QTLs reveal a complex archi-
tecture for each trait, with 4% of phenotypic variance
explained on average, and a maximum of 11%. A polygenic
architecture with many low-effect QTLs and a few medium-
effect QTLs is fairly consistent with the results of Bouchet
et al. (2017) and Pan et al. (2017) for leaf architecture, yield
and flowering, and those of Trachsel et al. (2016) and Cer-
rudo et al. (2018) for normalized difference vegetation index
dynamics, senescence and yield. In addition, our results
show that individually, GLAI QTLs explain a larger pro-
portion of the phenotypic variance at the network scale than
yield or flowering QTLs, which is in good agreement with
the high stability of the GLAI QTLs between the different
This will thus allow a better breeding efficiency through indi-
rect selection.

Multivariate GWAS is the most powerful approach

With the combination of very dense genotyping data
(7,806,995 SNPs) and two GWAS approaches, many associa-
tions (167,854) were detected for 11 traits measured in four
environments. The multivariate approach (M<sub>MV</sub>) detected
many more associations, including most of the associations
revealed by the univariate approach (M<sub>UV</sub>, Fig. S5). Moreo-
ver, M<sub>MV</sub> revealed many more associations common to several
traits (79% of associations) than M<sub>UV</sub> and explained a greater
part of the phenotypic variance, particularly for yield (Fig. 5).

The power gained in GWAS with the multivariate approach compared to the univariate approach has been demonstrated by several studies (Zhu and Zhang 2009; O'Reilly et al. 2012; Stephens 2013; Galesloot et al. 2014; Porter and O'Reilly 2017). Taking into account the covariances among traits brings additional information that increases the power in detecting genetic variants involved in the control of several traits, but also variants affecting a single trait (Ferreira and Purcell 2009; Stephens 2013; Carlson et al. 2019). However, the use of multivariate models may be limited by the number of phenotypes that can be integrated. For example, in this study, a maximum of eight traits could be integrated into the multivariate model fitted with the GEMMA software. In recent years, many tools have been developed for multi-trait GWAS that differ in terms of the number of analyzable traits, calculation times and assumptions on the distribution of phenotypes (O'Reilly et al. 2012; Sluis et al. 2013; Stephens 2013; Zhou and Stephens 2014). These tools also have a very variable power depending on the phenotypic correlations and the genetic architecture of the traits studied (Porter and O'Reilly 2017).

In this work, we used a two-step approach for its com-
putational efficiency at the cost of a probable loss of infor-
mation between the two steps. The gain in power noted in
several studies using functional mapping, and especially
random regression (Wu and Lin 2006; Li et al. 2010; Ning
et al. 2017, 2018; Campbell et al. 2019; Moreira et al. 2020)
may have better enhanced our understanding of the timing
of the genetic control of GLAI dynamics in response to
drought stress. However, computationally it would not have
been practicable with our dense genotyping data and the
backward-elimination method (Ning et al. 2017; Moreira
et al. 2020).

GLAI and agronomic traits exhibit polygenic
architecture and common genetic determinants

The associations we detected overlap strongly with the
genetic determinisms of GLAI traits, yield and flowering

(Fig. S6). Our results, like those of Cheverud (1988), Bou-
chet et al. (2017), Pan et al. (2017), Sodini et al. (2018)
and Yuan et al. (2019) show that correlations at the phe-
notypic level are good indicators of the sharing of genetic
determinants (pleiotropy or genetic linkage) between two
traits (Fig. S10A). This indicates that it is possible to select
secondary yield traits quite efficiently based on phenotypic
correlations. However, we observed only a moderate link
between phenotypic and genetic correlations (Fig. S10B),
which demonstrates the complementarity of these two pieces
of information. This partial discrepancy is probably due to
permanent environmental effects at the plot level throughout
the cycle, captured by the temporal approach implemented in
this study (Hadfield et al. 2007; Canela-Xandri et al. 2018).

Common genetic determinisms between earliness, plant
architecture and changes in leaf area during the cycle are
well known (Li et al. 2016b; Bouchet et al. 2017; Pan et al.
2017; Condorelli et al. 2018; Spindel et al. 2018). The entire
foliar and reproductive apparatus can be traced back to the
activity of the apical meristem and the floral meristems it
initiates (Kwiatkowska 2008), which is probably the reason
for these common determinisms (Thompson et al. 2015;
Baute et al. 2015). The impact of changes in leaf area on
light interception, gas exchange and water fluxes explains
the large number of QTLs affecting both yield and leaf area
dynamics under optimal or stressed conditions (Trachsel
et al. 2016; Cerrudo et al. 2018; Condorelli et al. 2018).
In our study, however, there is less colocalization between
yield and flowering than between yield and GLAI traits,
which could be the consequence of purposely reducing the
flowering time range of the panel. Having limited the flower-
ing date confounding effect, it can be assumed that a large
proportion of the common QTLs between yield and GLAI
traits that do not colocalize with flowering are due to direct
QTL effects on the GLAI. Moreover, we showed that some
GLAI QTLs detected under optimal conditions colocalize
with yield QTLs detected under stressed conditions. This
means yield under drought stress is partly due to GLAI
QTLs detected in the absence of stress, which may facilitate
the prediction of grain yield in breeding.

Individually, the detected QTLs reveal a complex archi-
tecture for each trait, with 4% of phenotypic variance
explained on average, and a maximum of 11%. A polygenic
architecture with many low-effect QTLs and a few medium-
effect QTLs is fairly consistent with the results of Bouchet
et al. (2017) and Pan et al. (2017) for leaf architecture, yield
and flowering, and those of Trachsel et al. (2016) and Cer-
rudo et al. (2018) for normalized difference vegetation index
dynamics, senescence and yield. In addition, our results
show that individually, GLAI QTLs explain a larger pro-
portion of the phenotypic variance at the network scale than
yield or flowering QTLs, which is in good agreement with
the high stability of the GLAI QTLs between the different
environments. Similarly, the total $r^2$ explained by all QTLs is greater for GLAI (AUC<sub>F</sub>, AUC<sub>S</sub>, GLAI<sub>M</sub>) than for yield (Fig. 5, M<sub>UV</sub>), while more QTLs were detected for yield (+25%). It therefore appears that the genetic architecture of GLAI is simpler than that of yield with fewer genetic factors, more stable between environments, exhibiting somewhat larger individual effects.

## GLAI dynamics is a promising secondary trait for drought tolerance

Sixteen GLAI QTLs were shown to explain GY measured in 11 environments under drought stress conditions (GY<sub>11</sub>). They explained a significant part of the variance of GY<sub>11</sub> (18%), despite the small size of the initial experimental network in which the QTLs were detected (only four environments, two of which under optimal conditions). It is interesting to note that the 16 QTLs had much more variable effects between the 11 environments of GY<sub>11</sub> (Fig. S9) than in the four relatively similar environments where GLAI data was collected (Fig. S7). GLAI QTLs detected in a specific environment are therefore informative of observed yield in a wide range of contrasting environments.

Among the 16 QTLs associated with GY<sub>11</sub>, five have similar effects in all environments, such as C1:57,391,598 whose allele from B73 increases grain yield. This QTL colocalizes with a QTL identified in the reference network whose B73 allele increases AUC<sub>S</sub> and GLAI<sub>M</sub>. It also colocalizes with a stay-green QTL found by Wang et al. (2016) and two leaf-width QTLs detected by Wang et al. (2018). A gain in GY when there is more green leaf area at flowering and during grain filling is consistent with the observations made on the reference network. Another example is C5:10,876,996, whose allele from B73 leads to a yield decrease, colocalizing with a QTL whose B73 allele increases D<sub>S</sub> and decreases GY under drought stress in the reference network. Consistently, Almeida et al. (2014) previously detected three QTLs affecting chlorophyll content during grain filling that also colocalize with C5:10,876,996. However, as observed by Trachsel et al. (2017), one would expect that increasing the duration of filling would lead to a yield increase, at least in environments without severe terminal stress. But, as D<sub>S</sub> is negatively correlated with D<sub>V</sub>, AUC<sub>V</sub>, and GLAI<sub>M</sub>, the decreased vigor of vegetative development might explain the yield loss (Fig. S3). Finally, QTL C3:19,131,545 is typical of QTLs with unstable effects between environments. Indeed, the reference allele of this QTL has the strongest GY effect detected in this study, varying from -6.9 to +3.1 q.ha<sup>-1</sup> (-14.4 to +4.8% of the mean yield) for 17SZEGEDWD and 16BLOISWD, respectively. It colocalizes with QTLs whose reference allele increases FF, AUC<sub>V</sub>, AUC<sub>S</sub>, and/or GLAI<sub>M</sub> under optimal and stressed conditions and GY under optimal conditions in the reference network. This QTL

therefore reveals the variable effect of a developed leaf area over the entire cycle. While both environments experienced drought stress during flowering and grain filling, it appears that 17SZEGEDWD where B73 allele has the strongest del-eterious effect on GY, was much warmer than 16BLOISWD where B73 allele increases GY (23 °C versus 18 °C in average between flowering and harvest). It would therefore seem that thermal stress coupled with drought stress is the reason for the negative effect of a more developed foliar area on GY, possibly cause by increased evapotranspiration induced by the rise in temperature. The interpretations proposed here for these three QTLs remain tentative, because although GLAI traits show little G×E interaction in the reference network, it is possible that the effect of these QTLs on GLAI changes direction within the nine additional environments considered here.

All these 16 QTL colocalize with previously identified QTL or meta-QTL for leaf architecture, plant architecture and/or GY and GY components (Trachsel et al. 2016; Zhang et al. 2016, 2017; Chen et al. 2017; Zhao et al. 2018). Validated on this 11 WD trials network, these GLAI QTL affecting GY seem interesting candidates to be further investigated.

## Perspectives and implications for breeding

In this study, we describe six new traits related to GLAI dynamics that constitutes relevant secondary traits for GY improvement under optimal and drought conditions. The other 18 GLAI traits extracted from the dynamics that were not analyzed in detail in this study show overall similar characteristics to those of the six traits analyzed (Fig. S4, Table S4). These traits could therefore be particularly interesting for new breeding programs, depending on the considered ideotype. Indeed, although the selection efficiency expected by using these traits individually can be close to that of yield alone, due to a limited gain in heritability or a genetic correlation with yield that can be reduced in certain cases, it as to be noted that using an ideotype selection index based on several of these traits can greatly increase the efficiency of indirect selection (Bänziger and Lafitte 1997; Monneveux et al. 2006; Ziyomo and Bernardo 2013). Additionally, it can be valuable to decrease selection accuracy if economic efficiency is increased in breeding programs (van Eeuwijk et al. 2019). Indeed, by selecting secondary traits such as GLAI dynamics, it is possible to quickly phenotype a large number of individuals at low cost early in the cycle, which could at least allow eliminating the individuals that are strongly divergent from the desired ideotype (Araus et al. 2018). Consistently, Bustos-Korts et al. (2019a) restated in a recent simulation study that intermediate secondary traits, such as green canopy dynamics, are particularly well suited
environments. Similarly, the total $r^2$ explained by all QTLs is greater for GLAI (AUC<sub>F</sub>, AUC<sub>S</sub>, GLAI<sub>M</sub>) than for yield (Fig. 5, M<sub>UV</sub>), while more QTLs were detected for yield (+25%). It therefore appears that the genetic architecture of GLAI is simpler than that of yield with fewer genetic factors, more stable between environments, exhibiting somewhat larger individual effects.

## GLAI dynamics is a promising secondary trait for drought tolerance

Sixteen GLAI QTLs were shown to explain GY measured in 11 environments under drought stress conditions (GY<sub>11</sub>). They explained a significant part of the variance of GY<sub>11</sub> (18%), despite the small size of the initial experimental network in which the QTLs were detected (only four environments, two of which under optimal conditions). It is interesting to note that the 16 QTLs had much more variable effects between the 11 environments of GY<sub>11</sub> (Fig. S9) than in the four relatively similar environments where GLAI data was collected (Fig. S7). GLAI QTLs detected in a specific environment are therefore informative of observed yield in a wide range of contrasting environments.

Among the 16 QTLs associated with GY<sub>11</sub>, five have similar effects in all environments, such as C1:57,391,598 whose allele from B73 increases grain yield. This QTL colocalizes with a QTL identified in the reference network whose B73 allele increases AUC<sub>S</sub> and GLAI<sub>M</sub>. It also colocalizes with a stay-green QTL found by Wang et al. (2016) and two leaf-width QTLs detected by Wang et al. (2018). A gain in GY when there is more green leaf area at flowering and during grain filling is consistent with the observations made on the reference network. Another example is C5:10,876,996, whose allele from B73 leads to a yield decrease, colocalizing with a QTL whose B73 allele increases D<sub>S</sub> and decreases GY under drought stress in the reference network. Consistently, Almeida et al. (2014) previously detected three QTLs affecting chlorophyll content during grain filling that also colocalize with C5:10,876,996. However, as observed by Trachsel et al. (2017), one would expect that increasing the duration of filling would lead to a yield increase, at least in environments without severe terminal stress. But, as D<sub>S</sub> is negatively correlated with D<sub>V</sub>, AUC<sub>V</sub>, and GLAI<sub>M</sub>, the decreased vigor of vegetative development might explain the yield loss (Fig. S3). Finally, QTL C3:19,131,545 is typical of QTLs with unstable effects between environments. Indeed, the reference allele of this QTL has the strongest GY effect detected in this study, varying from -6.9 to +3.1 q.ha<sup>-1</sup> (-14.4 to +4.8% of the mean yield) for 17SZEGEDWD and 16BLOISWD, respectively. It colocalizes with QTLs whose reference allele increases FF, AUC<sub>V</sub>, AUC<sub>S</sub>, and/or GLAI<sub>M</sub> under optimal and stressed conditions and GY under optimal conditions in the reference network. This QTL

therefore reveals the variable effect of a developed leaf area over the entire cycle. While both environments experienced drought stress during flowering and grain filling, it appears that 17SZEGEDWD where B73 allele has the strongest del-eterious effect on GY, was much warmer than 16BLOISWD where B73 allele increases GY (23 °C versus 18 °C in average between flowering and harvest). It would therefore seem that thermal stress coupled with drought stress is the reason for the negative effect of a more developed foliar area on GY, possibly cause by increased evapotranspiration induced by the rise in temperature. The interpretations proposed here for these three QTLs remain tentative, because although GLAI traits show little G×E interaction in the reference network, it is possible that the effect of these QTLs on GLAI changes direction within the nine additional environments considered here.

All these 16 QTL colocalize with previously identified QTL or meta-QTL for leaf architecture, plant architecture and/or GY and GY components (Trachsel et al. 2016; Zhang et al. 2016, 2017; Chen et al. 2017; Zhao et al. 2018). Validated on this 11 WD trials network, these GLAI QTL affecting GY seem interesting candidates to be further investigated.

## Perspectives and implications for breeding

In this study, we describe six new traits related to GLAI dynamics that constitutes relevant secondary traits for GY improvement under optimal and drought conditions. The other 18 GLAI traits extracted from the dynamics that were not analyzed in detail in this study show overall similar characteristics to those of the six traits analyzed (Fig. S4, Table S4). These traits could therefore be particularly interesting for new breeding programs, depending on the considered ideotype. Indeed, although the selection efficiency expected by using these traits individually can be close to that of yield alone, due to a limited gain in heritability or a genetic correlation with yield that can be reduced in certain cases, it as to be noted that using an ideotype selection index based on several of these traits can greatly increase the efficiency of indirect selection (Bänziger and Lafitte 1997; Monneveux et al. 2006; Ziyomo and Bernardo 2013). Additionally, it can be valuable to decrease selection accuracy if economic efficiency is increased in breeding programs (van Eeuwijk et al. 2019). Indeed, by selecting secondary traits such as GLAI dynamics, it is possible to quickly phenotype a large number of individuals at low cost early in the cycle, which could at least allow eliminating the individuals that are strongly divergent from the desired ideotype (Araus et al. 2018). Consistently, Bustos-Korts et al. (2019a) restated in a recent simulation study that intermediate secondary traits, such as green canopy dynamics, are particularly well suited
for field evaluation during early breeding cycles to improve selection efficiency.

To improve breeding efficiency, the many QTL identified in this work could be used in a marker-assisted selection program. However, it can be noted that the most heritable traits with the least G×E interaction were not necessarily the best explained by the detected QTLs in the multivariate GWAS approach (for example, $D_v$). This suggests that these traits probably have a polygenic architecture with very weak effects, as is postulated for plant height (Peiffer et al. 2014). In this context, a genomic selection approach would undoubtedly be more effective than an approach based on QTL detection to increase the efficiency of breeding programs (Spindel et al. 2015), as shown by Bouidghaghen et al. (2023) for $D_v$. Indeed, the use of secondary traits in genomic selection to predict yield greatly increases the quality of prediction (Rutkoski et al. 2016; Sun et al. 2017; Crain et al. 2018; Bustos-Korts et al. 2019a). However, genomic selection should not be substitute to QTL detection, which may be more efficient in certain cases (Spindel et al. 2015; Yuan et al. 2019). Cerrudo et al. (2018) proposed an optimal approach of using the QTLs identified by GWAS in the first stages of a breeding program to enrich the material with strong-effect alleles, then using genomic selection to improve grain yield by taking advantage of low-effect loci. Alternatively, the integration of previously detected QTLs into a genomic prediction model has also performed well, provided the QTLs explain more than 10% of the variance (Bernardo 2014; Spindel et al. 2016).

Finally, the best strategies to take advantage of these new secondary traits will only be identified by deepening our knowledge and understanding of how and when they confer adaption to a specific environment type, probably through a mixture of ecophysiological modeling and empirical validation.

## Coupling ecophysiological and genetic models to gain in knowledge, understanding and power

More and more approaches combine ecophysiological modeling and genetic models through a framework that can go from sequential to fully nested. This coupling allows all these approaches to gain knowledge, understanding and/or predictive ability depending on whether the emphasis is placed on ecophysiological or genetic modeling (Messina et al. 2023).

Our two-step approach uses a simple physiological model and high-throughput field phenotyping to extract the components of GLAI dynamics and then characterize their genetic determinism. While the main focus is on genetic modeling, the use of a physiological model has several advantages: the dynamic modeling of GLAI allows (i) to extract descriptive genotypic parameters for different physiological phases,

which is highly desired since the correlation between secondary traits and GY varies during the cycle depending on the type of environment, and (ii) to increase the heritability of the extracted traits (Bustos-Korts et al. 2019a, b). The main contribution of this work is a better knowledge of the genetic determinism of GLAI dynamics and its link with GY made possible by fine field phenotyping and a computationally efficient multi-trait and multi-environment QTL approach for a large number of genotypes and a dense genotyping. Indeed, this study is one of the first to study the genetic determinism of maize GLAI temporal evolution over the whole cycle in the field. As seen in similar works (Millet et al. 2016, 2019; Touzy et al. 2019; Bouidghaghen et al. 2023), this coupled modeling approach allowed us to gain power in GWAS but also to better understand the physiological mechanisms relevant for adaptation to a particular environment through the analysis of correlations between traits and colocalization of QTL (Bustos-Korts et al. 2019b), which constitutes a first step toward ideotype design. This understanding is however limited by the limited range of environments and the relationships between evaluated traits.

Studies that put more effort into ecophysiological modeling are able to greatly expand our understanding of G×E interactions and adaptive traits for each type of environment and therefore propose one or more ideotypes when necessary (Bustos-Korts et al. 2019b; Cooper and Messina 2023; Messina et al. 2023). Today, the most advanced approaches to couple ecophysiological and genetic modeling (Technow et al. 2015; Cooper et al. 2016; Messina et al. 2018; Bustos-Korts et al. 2019b) consist of a one-step approach that integrates biological information from several environments through an ecophysiological model, directly into genomic prediction algorithms to increase the predictive power and understanding of G×E interactions in the breeding network. While they surely constitute the next step to take advantage of a better G×E understanding to improve the breeding efficiency, this kind of approaches remains complementary with approaches like our, based on targeted phenotyping in a well-characterized environment, because they allow to gain empirical knowledge which can then enrich the ecophysiological and genetic model couple (Yin et al. 2004; Bustos-Korts et al. 2019a; Hammer et al. 2019).

## Conclusion

Genetic determinants of maize canopy characteristics were found by measuring GLAI dynamics over the crop cycle in a diverse panel from a MAGIC population genotyped at a very high density. Due to the difficulty of monitoring GLAI in the field for a large number of genotypes, this study is one of the first to study the genetic determinism of its temporal evolution over the whole cycle. From the GLAI dynamics,
for field evaluation during early breeding cycles to improve selection efficiency.

To improve breeding efficiency, the many QTL identified in this work could be used in a marker-assisted selection program. However, it can be noted that the most heritable traits with the least G×E interaction were not necessarily the best explained by the detected QTLs in the multivariate GWAS approach (for example, $D_v$). This suggests that these traits probably have a polygenic architecture with very weak effects, as is postulated for plant height (Peiffer et al. 2014). In this context, a genomic selection approach would undoubtedly be more effective than an approach based on QTL detection to increase the efficiency of breeding programs (Spindel et al. 2015), as shown by Bouidghaghen et al. (2023) for $D_v$. Indeed, the use of secondary traits in genomic selection to predict yield greatly increases the quality of prediction (Rutkoski et al. 2016; Sun et al. 2017; Crain et al. 2018; Bustos-Korts et al. 2019a). However, genomic selection should not be substitute to QTL detection, which may be more efficient in certain cases (Spindel et al. 2015; Yuan et al. 2019). Cerrudo et al. (2018) proposed an optimal approach of using the QTLs identified by GWAS in the first stages of a breeding program to enrich the material with strong-effect alleles, then using genomic selection to improve grain yield by taking advantage of low-effect loci. Alternatively, the integration of previously detected QTLs into a genomic prediction model has also performed well, provided the QTLs explain more than 10% of the variance (Bernardo 2014; Spindel et al. 2016).

Finally, the best strategies to take advantage of these new secondary traits will only be identified by deepening our knowledge and understanding of how and when they confer adaption to a specific environment type, probably through a mixture of ecophysiological modeling and empirical validation.

## Coupling ecophysiological and genetic models to gain in knowledge, understanding and power

More and more approaches combine ecophysiological modeling and genetic models through a framework that can go from sequential to fully nested. This coupling allows all these approaches to gain knowledge, understanding and/or predictive ability depending on whether the emphasis is placed on ecophysiological or genetic modeling (Messina et al. 2023).

Our two-step approach uses a simple physiological model and high-throughput field phenotyping to extract the components of GLAI dynamics and then characterize their genetic determinism. While the main focus is on genetic modeling, the use of a physiological model has several advantages: the dynamic modeling of GLAI allows (i) to extract descriptive genotypic parameters for different physiological phases,

which is highly desired since the correlation between secondary traits and GY varies during the cycle depending on the type of environment, and (ii) to increase the heritability of the extracted traits (Bustos-Korts et al. 2019a, b). The main contribution of this work is a better knowledge of the genetic determinism of GLAI dynamics and its link with GY made possible by fine field phenotyping and a computationally efficient multi-trait and multi-environment QTL approach for a large number of genotypes and a dense genotyping. Indeed, this study is one of the first to study the genetic determinism of maize GLAI temporal evolution over the whole cycle in the field. As seen in similar works (Millet et al. 2016, 2019; Touzy et al. 2019; Bouidghaghen et al. 2023), this coupled modeling approach allowed us to gain power in GWAS but also to better understand the physiological mechanisms relevant for adaptation to a particular environment through the analysis of correlations between traits and colocalization of QTL (Bustos-Korts et al. 2019b), which constitutes a first step toward ideotype design. This understanding is however limited by the limited range of environments and the relationships between evaluated traits.

Studies that put more effort into ecophysiological modeling are able to greatly expand our understanding of G×E interactions and adaptive traits for each type of environment and therefore propose one or more ideotypes when necessary (Bustos-Korts et al. 2019b; Cooper and Messina 2023; Messina et al. 2023). Today, the most advanced approaches to couple ecophysiological and genetic modeling (Technow et al. 2015; Cooper et al. 2016; Messina et al. 2018; Bustos-Korts et al. 2019b) consist of a one-step approach that integrates biological information from several environments through an ecophysiological model, directly into genomic prediction algorithms to increase the predictive power and understanding of G×E interactions in the breeding network. While they surely constitute the next step to take advantage of a better G×E understanding to improve the breeding efficiency, this kind of approaches remains complementary with approaches like our, based on targeted phenotyping in a well-characterized environment, because they allow to gain empirical knowledge which can then enrich the ecophysiological and genetic model couple (Yin et al. 2004; Bustos-Korts et al. 2019a; Hammer et al. 2019).

## Conclusion

Genetic determinants of maize canopy characteristics were found by measuring GLAI dynamics over the crop cycle in a diverse panel from a MAGIC population genotyped at a very high density. Due to the difficulty of monitoring GLAI in the field for a large number of genotypes, this study is one of the first to study the genetic determinism of its temporal evolution over the whole cycle. From the GLAI dynamics,
24 traits were derived, six of which were shown to capture a significant part of the panel diversity. Six GLAI traits cor- related with yield under optimal and drought conditions are more variable, more heritable, and less subject to G×E interaction than yield, making them particularly suitable secondary traits for drought tolerance improvement. The multivariate GWAS approach explained a large proportion (40–60%) of the traits variance. The QTLs identified show a strong overlap between the genetic determinants of GLAI, GY and flowering. A subset of GLAI QTLs was validated as explaining nearly one-fifth of the GY variation measured on 11 environments under drought stress. GLAI traits are promising for the relative simplicity of their genetic archi- tecture and medium genetic effects that are relatively stable between environments.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00122-024-04572-6.

**Acknowledgements** The authors are grateful to the ANRT (Association Nationale de la Recherche et de la Technologie) for the financial support of the PhD thesis (2015/1190). The authors would also like to thank Brigitte Mangin and Gaëtan Touzy for the helpful discussions during data analysis and Rachel Carol (Bioscience Editing, France) for English editing, as well as the reviewers for their precious comments.

**Author contribution statement** JB, SP, PD and FB contributed to the study conception and design. Data analysis was performed by JB, CB and M-HT. The first draft of the manuscript was written by JB, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

**Funding** This work was supported by Association Nationale de la Recherche et de la Technologie, grant N°2015/1190.

**Data availability** Data and plant material are available on reasonable request.

## Declarations

**Conflict of interest** 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

Almeida GD, Nair S, Borém A et al (2014) Molecular mapping across three populations reveals a QTL hotspot region on chromosome 3 for secondary traits associated with drought tolerance in tropical maize. Mol Breed 34:701–715. https://doi.org/10.1007/s11032-014-0068-5

Araus JL, Serret MD, Edmeades G (2012) Phenotyping maize for adaptation to drought. Front Physiol 3:305. https://doi.org/10.3389/fphys.2012.00305

Araus JL, Kefauver SC, Zaman-Allah M et al (2018) Translating high-throughput phenotyping into genetic gain. Trends Plant Sci 23:451–466. https://doi.org/10.1016/j.tplants.2018.02.001

Bänziger M, Lafitte HR (1997) Efficiency of secondary traits for improving maize for low-nitrogen target environments. Crop Sci. https://doi.org/10.2135/cropsci1997.0011183X003700040013x

Baute J, Herman D, Coppens F et al (2015) Correlation analysis of the transcriptome of growing leaves with mature leaf parameters in a maize RIL population. Genome Biol 16:168. https://doi.org/10.1186/s13059-015-0735-9

Bernardo R (2014) Genomewide selection when major genes are known. Crop Sci 54:68–75. https://doi.org/10.2135/cropsci2013.05.0315

Beyene Y, Semagn K, Mugo S et al (2016) Performance and grain yield stability of maize populations developed using marker-assisted recurrent selection and pedigree selection procedures. Euphytica Neth J Plant Breed 208:285–297. https://doi.org/10.1007/s10681-015-1590-1

Blanco J, Dutartre D, Tixier M-H et al (2019) A high-throughput model-assisted method for phenotyping maize green leaf area index dynamics using unmanned aerial vehicle imagery. Front Plant Sci 10:685. https://doi.org/10.3389/fpls.2019.00685

Bolaños J, Edmeades GO (1996) The importance of the anthesis-silking interval in breeding for drought tolerance in tropical maize. Field Crops Res 48:65–80. https://doi.org/10.1016/0378-4290(96)00036-6

Borrás L, Maddonni GA, Otegui ME (2003) Leaf senescence in maize hybrids: plant population, row spacing and kernel set effects. Field Crops Res 82:13–26. https://doi.org/10.1016/S0378-4290(03)00002-9

Borrell AK, Hammer GL, Henzell RG (2000) Does maintaining green leaf area in sorghum improve yield under drought? II. dry matter production and yield. Crop Sci 40:1037–1048. https://doi.org/10.2135/cropsci2000.4041037x

Borrell AK, van Oosterom EJ, Mullet JE et al (2014) Stay-green alleles individually enhance grain yield in sorghum under drought by modifying canopy development and water uptake patterns. New Phytol 203:817–830. https://doi.org/10.1111/nph.12869

Bouchet S, Bertin P, Presterl T et al (2017) Association mapping for phenology and plant architecture in maize shows higher power for developmental traits compared with growth influenced traits. Heredity 118:249–259. https://doi.org/10.1038/hdy.2016.88

Bouidghaghen J, Moreau L, Beauchène K et al (2023) Robotized indoor phenotyping allows genomic prediction of adaptive traits in the field. Nat Commun 14:6603. https://doi.org/10.1038/s41467-023-42298-z

Broman KW, Wu H, Sen S, Churchill GA (2003) R/qlt: QTL mapping in experimental crosses. Bioinformatics 19:889–890. https://doi.org/10.1093 bioinformatics/btg112

Buet C, Dubreuil P, Tixier M-H, et al (2013) The molecular characterization of a MAGIC population reveals high potential for gene discovery. Poster session presented at the 55th Annual Maize Genetics Conference, St. Charles, Illinois

Bustos-Korts D, Boer MP, Malosetti M et al (2019a) Combining crop growth modeling and statistical genetic modeling to evaluate

<figure><img src="image_11.png" /></figure>
24 traits were derived, six of which were shown to capture a significant part of the panel diversity. Six GLAI traits cor- related with yield under optimal and drought conditions are more variable, more heritable, and less subject to G×E interaction than yield, making them particularly suitable secondary traits for drought tolerance improvement. The multivariate GWAS approach explained a large proportion (40–60%) of the traits variance. The QTLs identified show a strong overlap between the genetic determinants of GLAI, GY and flowering. A subset of GLAI QTLs was validated as explaining nearly one-fifth of the GY variation measured on 11 environments under drought stress. GLAI traits are promising for the relative simplicity of their genetic archi- tecture and medium genetic effects that are relatively stable between environments.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00122-024-04572-6.

**Acknowledgements** The authors are grateful to the ANRT (Association Nationale de la Recherche et de la Technologie) for the financial support of the PhD thesis (2015/1190). The authors would also like to thank Brigitte Mangin and Gaëtan Touzy for the helpful discussions during data analysis and Rachel Carol (Bioscience Editing, France) for English editing, as well as the reviewers for their precious comments.

**Author contribution statement** JB, SP, PD and FB contributed to the study conception and design. Data analysis was performed by JB, CB and M-HT. The first draft of the manuscript was written by JB, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

**Funding** This work was supported by Association Nationale de la Recherche et de la Technologie, grant N°2015/1190.

**Data availability** Data and plant material are available on reasonable request.

## Declarations

**Conflict of interest** 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

Almeida GD, Nair S, Borém A et al (2014) Molecular mapping across three populations reveals a QTL hotspot region on chromosome 3 for secondary traits associated with drought tolerance in tropical maize. Mol Breed 34:701–715. https://doi.org/10.1007/s11032-014-0068-5

Araus JL, Serret MD, Edmeades G (2012) Phenotyping maize for adaptation to drought. Front Physiol 3:305. https://doi.org/10.3389/fphys.2012.00305

Araus JL, Kefauver SC, Zaman-Allah M et al (2018) Translating high-throughput phenotyping into genetic gain. Trends Plant Sci 23:451–466. https://doi.org/10.1016/j.tplants.2018.02.001

Bänziger M, Lafitte HR (1997) Efficiency of secondary traits for improving maize for low-nitrogen target environments. Crop Sci. https://doi.org/10.2135/cropsci1997.0011183X003700040013x

Baute J, Herman D, Coppens F et al (2015) Correlation analysis of the transcriptome of growing leaves with mature leaf parameters in a maize RIL population. Genome Biol 16:168. https://doi.org/10.1186/s13059-015-0735-9

Bernardo R (2014) Genomewide selection when major genes are known. Crop Sci 54:68–75. https://doi.org/10.2135/cropsci2013.05.0315

Beyene Y, Semagn K, Mugo S et al (2016) Performance and grain yield stability of maize populations developed using marker-assisted recurrent selection and pedigree selection procedures. Euphytica Neth J Plant Breed 208:285–297. https://doi.org/10.1007/s10681-015-1590-1

Blanco J, Dutartre D, Tixier M-H et al (2019) A high-throughput model-assisted method for phenotyping maize green leaf area index dynamics using unmanned aerial vehicle imagery. Front Plant Sci 10:685. https://doi.org/10.3389/fpls.2019.00685

Bolaños J, Edmeades GO (1996) The importance of the anthesis-silking interval in breeding for drought tolerance in tropical maize. Field Crops Res 48:65–80. https://doi.org/10.1016/0378-4290(96)00036-6

Borrás L, Maddonni GA, Otegui ME (2003) Leaf senescence in maize hybrids: plant population, row spacing and kernel set effects. Field Crops Res 82:13–26. https://doi.org/10.1016/S0378-4290(03)00002-9

Borrell AK, Hammer GL, Henzell RG (2000) Does maintaining green leaf area in sorghum improve yield under drought? II. dry matter production and yield. Crop Sci 40:1037–1048. https://doi.org/10.2135/cropsci2000.4041037x

Borrell AK, van Oosterom EJ, Mullet JE et al (2014) Stay-green alleles individually enhance grain yield in sorghum under drought by modifying canopy development and water uptake patterns. New Phytol 203:817–830. https://doi.org/10.1111/nph.12869

Bouchet S, Bertin P, Presterl T et al (2017) Association mapping for phenology and plant architecture in maize shows higher power for developmental traits compared with growth influenced traits. Heredity 118:249–259. https://doi.org/10.1038/hdy.2016.88

Bouidghaghen J, Moreau L, Beauchène K et al (2023) Robotized indoor phenotyping allows genomic prediction of adaptive traits in the field. Nat Commun 14:6603. https://doi.org/10.1038/s41467-023-42298-z

Broman KW, Wu H, Sen S, Churchill GA (2003) R/qlt: QTL mapping in experimental crosses. Bioinformatics 19:889–890. https://doi.org/10.1093 bioinformatics/btg112

Buet C, Dubreuil P, Tixier M-H, et al (2013) The molecular characterization of a MAGIC population reveals high potential for gene discovery. Poster session presented at the 55th Annual Maize Genetics Conference, St. Charles, Illinois

Bustos-Korts D, Boer MP, Malosetti M et al (2019a) Combining crop growth modeling and statistical genetic modeling to evaluate
phenotyping strategies. *Front Plant Sci* 10:1491. https://doi.org/
10.3389/fpls.2019.01491

Bustos-Korts D, Malosetti M, Chenu K et al (2019b) From QTLs to
adaptation landscapes: using genotype-to-phenotype models to
characterize G×E over time. *Front Plant Sci* 10:1540. https://doi.
org/10.3389/fpls.2019.01540

Butler DG, Cullis BR, Gilmour AR, Gogel BJ (2009) ASReml-R user guide release 3.0. VSN International Ltd. https://www.vsnco.uk, Hemel Hempstead

Cairns JE, Sanchez C, Vargas M et al (2012) Dissecting maize produc-
tivity: ideotypes associated with grain yield under drought stress
and well-watered conditions. J Integr Plant Biol 54:1007–1020.
https://doi.org/10.1111/j.1744-7909.2012.01156.x

Cairns JE, Crossa J, Zaidi PH et al (2013) Identification of drought, heat, and combined drought and heat tolerant donors in Maize. Crop Sci 53:1335–1346. https://doi.org/10.2135/cropsci2012.09.0545

Campbell MT, Du Q, Liu K et al (2017) A comprehensive image-
based phenomic analysis reveals the complex genetic architecture
of shoot growth dynamics in rice (Oryza sativa). Plant Genom
10:2016. https://doi.org/10.3835/plantgenome2016.07.0064

Campbell M, Walia H, Morota G (2018) Utilizing random regression models for genomic prediction of a longitudinal trait derived from high-throughput phenotyping. Plant Direct 2:e00080. https://doi.org/10.1002/pld3.80

Campbell M, Momen M, Walia H, Morota G (2019) Leveraging breeding values obtained from random regression models for genetic inference of longitudinal traits. Plant Genome 12:180075. https://doi.org/10.3835/plantgenome2018.10.0075

Canela-Xandri O, Rawlik K, Tenesa A (2018) An atlas of genetic asso-
ciations in UK Biobank. Nat Genet 50:1593–1599. https://doi.
org/10.1038/s41588-018-0248-z

Carlson MO, Montilla-Bascon G, Hoekenga OA et al (2019) Multivari-
ate genome-wide association analyses reveal the genetic basis of
seed fatty acid composition in oat (Avena sativa L.). G3 Genes
Genom Genetcs 9:2963–2975. https://doi.org/10.1534/g3.119.
400228

Cerrudo D, Cao S, Yuan Y et al (2018) Genomic selection outperforms marker assisted selection for grain yield and physiological traits in a maize doubled haploid population across water treatments. *Front Plant Sci* 9:366. https://doi.org/10.3389/fpls.2018.00366

Chen L, An Y, Li Y et al (2017) Candidate loci for yield-related traits in maize revealed by a combination of MetaQTL analysis and regional association mapping. Front Plant Sci 8:2190. https://doi.org/10.3389/fpls.2017.02190

Cheverud JM (1988) A comparison of genetic and phenotypic corre-
lations. Evolution 42:958-968. https://doi.org/10.1111/j.1558-5646.1988.tb02514.x

Christopher JT, Veyradier M, Borrell AK et al (2014) Phenotyping novel stay-green traits to capture genetic variation in senescence dynamics. Funct Plant Biol 41:1035–1048. https://doi.org/10.1071/FP14052

Christopher JT, Christopher MJ, Borrell AK et al (2016) Stay-green traits to improve wheat adaptation in well-watered and water-limited environments. J Exp Bot 67:5159–5172. https://doi.org/10.1093/jxb/erw276

Condorelli GE, Maccaferri M, Newcomb M et al (2018) Comparative aerial and ground based high throughput phenotyping for the genetic dissection of NDVI as a proxy for drought adaptive traits in durum wheat. Front Plant Sci 9:983. https://doi.org/10.3389/fpls.2018.00893

Cooper M, Messina CD (2023) Breeding crops for drought-affected environments and improved climate resilience. Plant Cell 35:162–186. https://doi.org/10.1093/plcell/koac321

Cooper M, Technow F, Messina C et al (2016) Use of crop growth models with whole-genome prediction: application to a maize

multienvironment trial. Crop Sci 56:2141-2156. https://doi.org/
10.2135/cropsci2015.08.0512

Cormier F, Le Gouis J, Dubreuil P et al (2014) A genome-wide iden-
tification of chromosomal regions determining nitrogen use effi-
ciency components in wheat (Triticum aestivum L.). Theor Appl
Genet 127:2679-2693

Crain J, Mondal S, Rutkoski J et al (2018) Combining high-throughput phenotyping and genomic information to increase prediction and selection accuracy in wheat breeding. Plant Genom 11:170043. https://doi.org/10.3835/plantgenome2017.05.0043

Crispim AC, Kelly MJ, Guimarães SEF et al (2015) Multi-trait GWAS and new candidate genes annotation for growth curve parameters in brahman cattle. PLoS ONE 10:e0139906. https://doi.org/10.1371/journal.pone.0139906

Cullis BR, Smith AB, Coombes NE (2006) On the design of early gen-
eration variety trials with correlated data. J Agric Biol Environ
Stat 11:381. https://doi.org/10.1198/108571106X154443

Duvick DN (2005) Genetic progress in yield of United States maize (Zea mays L.). Maydica 50:193-202

Ferreira MAR, Purcell SM (2009) A multivariate test of associa-
tion. Bioinformatics 25:132-133. https://doi.org/10.1093/bioin
formatics/btn563

Fischer KS, Edmeades GO, Johnson EC (1989) Selection for the improvement of maize yield under moisture-deficits. Field Crops Res 22:227-243. https://doi.org/10.1016/0378-4290(89)90094-4

Galesloot TE, Steen K, Kiemeney LALM et al (2014) A comparison of multivariate genome-wide association methods. PLoS ONE 9:e95923. https://doi.org/10.1371/journal.pone.0095923

Ganal MW, Durstewitz G, Polley A et al (2011) A large maize (Zea mays L.) SNP genotyping array: development and germplasm genotyping, and genetic mapping to compare with the B73 reference genome. PLoS ONE 6:e28334. https://doi.org/10.1371/journal.pone.0028334

Hadfield JD, Nutall A, Osorio D, Owens IPF (2007) Testing the pheno-
typic gambit: phenotypic, genetic and environmental correlations
of colour. J Evol Biol 20:549-557. https://doi.org/10.1111/j.
1420-9101.2006.01262.x

Hammer GL, McLean G, Chapman S et al (2014) Crop design for specific adaptation in variable dryland production environments. Crop Pasture Sci 65:614–626. https://doi.org/10.1071/CP14088

Hammer G, Messina C, Wu A, Cooper M (2019) Biological reality and parsimony in crop models—Why we need both in crop improvement! Silico Plants. https://doi.org/10.1093/insilicoplants/diz010

Harrison MT, Tardieu F, Dong Z et al (2014) Characterizing drought stress and trait influence on maize yield under current and future conditions. Glob Change Biol 20:867–878. https://doi.org/10.1111/gcb.12381

Hoegh-Guldberg O, Jacob D, Bindi M, et al (2018) Impacts of 1.5 C global warming on natural and human systems. Glob Warm 15 C IPCC Spec Rep

Howard JT, Jiao S, Tiezzi F et al (2015) Genome-wide association study on legendre random regression coefficients for the growth and feed intake trajectory on Duroc Boars. BMC Genet. https://doi.org/10.1186/s12863-015-0218-8

Huang S, Gao Y, Li Y et al (2017) Influence of plant architecture on maize physiology and yield in the Heilonggang River valley. Crop J 5:52–62. https://doi.org/10.1016/j.cj.2016.06.018

Huisman AE, Veerkamp RF, Arendonk JAM (2002) Genetic parame-
ters for various random regression models to describe the weight
data of pigs. J Anim Sci 80:575-582. https://doi.org/10.2527/
2002.803575x

Hurtado PX, Schnabel SK, Zaban A et al (2012) Dynamics of senes-
cence-related QTLs in potato. Euphytica 183:289-302. https://
doi.org/10.1007/s10681-011-0464-4
phenotyping strategies. *Front Plant Sci* 10:1491. https://doi.org/
10.3389/fpls.2019.01491

Bustos-Korts D, Malosetti M, Chenu K et al (2019b) From QTLs to
adaptation landscapes: using genotype-to-phenotype models to
characterize G×E over time. *Front Plant Sci* 10:1540. https://doi.
org/10.3389/fpls.2019.01540

Butler DG, Cullis BR, Gilmour AR, Gogel BJ (2009) ASReml-R user guide release 3.0. VSN International Ltd. https://www.vsn.co.uk, Hemel Hempstead

Cairns JE, Sanchez C, Vargas M et al (2012) Dissecting maize produc-
tivity: ideotypes associated with grain yield under drought stress
and well-watered conditions. J Integr Plant Biol 54:1007–1020.
https://doi.org/10.1111/j.1744-7909.2012.01156.x

Cairns JE, Crossa J, Zaidi PH et al (2013) Identification of drought, heat, and combined drought and heat tolerant donors in Maize. Crop Sci 53:1335–1346. https://doi.org/10.2135/cropsci2012.09.0545

Campbell MT, Du Q, Liu K et al (2017) A comprehensive image-
based phenomic analysis reveals the complex genetic architecture
of shoot growth dynamics in rice (Oryza sativa). Plant Genom
10:2016. https://doi.org/10.3835/plantgenome2016.07.0064

Campbell M, Walia H, Morota G (2018) Utilizing random regression models for genomic prediction of a longitudinal trait derived from high-throughput phenotyping. Plant Direct 2:e00080. https://doi.org/10.1002/pld3.80

Campbell M, Momen M, Walia H, Morota G (2019) Leveraging breeding values obtained from random regression models for genetic inference of longitudinal traits. Plant Genome 12:180075. https://doi.org/10.3835/plantgenome2018.10.0075

Canela-Xandri O, Rawlik K, Tenesa A (2018) An atlas of genetic asso-
ciations in UK Biobank. Nat Genet 50:1593–1599. https://doi.
org/10.1038/s41588-018-0248-z

Carlson MO, Montilla-Bascon G, Hoekenga OA et al (2019) Multivariate genome-wide association analyses reveal the genetic basis of seed fatty acid composition in oat (Avena sativa L.). G3 Genes Genom Genetcs 9:2963–2975. https://doi.org/10.1534/g3.119.400228

Cerrudo D, Cao S, Yuan Y et al (2018) Genomic selection outperforms marker assisted selection for grain yield and physiological traits in a maize doubled haploid population across water treatments. *Front Plant Sci* 9:366. https://doi.org/10.3389/fpls.2018.00366

Chen L, An Y, Li Y et al (2017) Candidate loci for yield-related traits in maize revealed by a combination of MetaQTL analysis and regional association mapping. Front Plant Sci 8:2190. https://doi.org/10.3389/fpls.2017.02190

Cheverud JM (1988) A comparison of genetic and phenotypic corre-
lations. Evolution 42:958-968. https://doi.org/10.1111/j.1558-5646.1988.tb02514.x

Christopher JT, Veyradier M, Borrell AK et al (2014) Phenotyping novel stay-green traits to capture genetic variation in senescence dynamics. Funct Plant Biol 41:1035–1048. https://doi.org/10.1071/FP14052

Christopher JT, Christopher MJ, Borrell AK et al (2016) Stay-green traits to improve wheat adaptation in well-watered and water-limited environments. J Exp Bot 67:5159–5172. https://doi.org/10.1093/jxb/erw276

Condorelli GE, Maccaferri M, Newcomb M et al (2018) Comparative aerial and ground based high throughput phenotyping for the genetic dissection of NDVI as a proxy for drought adaptive traits in durum wheat. Front Plant Sci 9:983. https://doi.org/10.3389/fpls.2018.00893

Cooper M, Messina CD (2023) Breeding crops for drought-affected environments and improved climate resilience. Plant Cell 35:162–186. https://doi.org/10.1093/plcell/koac321

Cooper M, Technow F, Messina C et al (2016) Use of crop growth models with whole-genome prediction: application to a maize

multienvironment trial. Crop Sci 56:2141–2156. https://doi.org/
10.2135/cropsci2015.08.0512

Cormier F, Le Gouis J, Dubreuil P et al (2014) A genome-wide ident-
fication of chromosomal regions determining nitrogen use effi-
ciency components in wheat (Triticum aestivum L.). Theor Appl
Genet 127:2679–2693

Crain J, Mondal S, Rutkoski J et al (2018) Combining high-throughput phenotyping and genomic information to increase prediction and selection accuracy in wheat breeding. Plant Genom 11:170043. https://doi.org/10.3835/plantgenome2017.05.0043

Crispim AC, Kelly MJ, Guimarães SEF et al (2015) Multi-trait GWAS and new candidate genes annotation for growth curve parameters in brahman cattle. PLoS ONE 10:e0139906. https://doi.org/10.1371/journal.pone.0139906

Cullis BR, Smith AB, Coombes NE (2006) On the design of early gen-
eration variety trials with correlated data. J Agric Biol Environ
Stat 11:381. https://doi.org/10.1198/108571106X154443

Duvick DN (2005) Genetic progress in yield of United States maize (Zea mays L.). Maydica 50:193-202

Ferreira MAR, Purcell SM (2009) A multivariate test of associa-
tion. Bioinformatics 25:132–133. https://doi.org/10.1093/bioin
formatics/btn563

Fischer KS, Edmeades GO, Johnson EC (1989) Selection for the improvement of maize yield under moisture-deficits. Field Crops Res 22:227–243. https://doi.org/10.1016/0378-4290(89)90094-4

Galesloot TE, Steen K, Kiemeney LALM et al (2014) A comparison of multivariate genome-wide association methods. PLoS ONE 9:e95923. https://doi.org/10.1371/journal.pone.0095923

Ganal MW, Durstewitz G, Polley A et al (2011) A large maize (Zea mays L.) SNP genotyping array: development and germplasm genotyping, and genetic mapping to compare with the B73 reference genome. PLoS ONE 6:e28334. https://doi.org/10.1371/journal.pone.0028334

Hadfield JD, Nutall A, Osorio D, Owens IPF (2007) Testing the pheno-
typic gambit: phenotypic, genetic and environmental correlations
of colour. J Evol Biol 20:549–557. https://doi.org/10.1111/j.
1420-9101.2006.01262.x

Hammer GL, McLean G, Chapman S et al (2014) Crop design for specific adaptation in variable dryland production environments. Crop Pasture Sci 65:614–626. https://doi.org/10.1071/CP14088

Hammer G, Messina C, Wu A, Cooper M (2019) Biological reality and parsimony in crop models—Why we need both in crop improvement! Silico Plants. https://doi.org/10.1093/insilicoplants/diz010

Harrison MT, Tardieu F, Dong Z et al (2014) Characterizing drought stress and trait influence on maize yield under current and future conditions. Glob Change Biol 20:867–878. https://doi.org/10.1111/gcb.12381

Hoegh-Guldberg O, Jacob D, Bindi M, et al (2018) Impacts of 1.5 C global warming on natural and human systems. Glob Warm 15 C IPCC Spec Rep

Howard JT, Jiao S, Tiezzi F et al (2015) Genome-wide association study on legendre random regression coefficients for the growth and feed intake trajectory on Duroc Boars. BMC Genet. https://doi.org/10.1186/s12863-015-0218-8

Huang S, Gao Y, Li Y et al (2017) Influence of plant architecture on maize physiology and yield in the Heilonggang River valley. Crop J 5:52–62. https://doi.org/10.1016/j.cj.2016.06.018

Huisman AE, Veerkamp RF, Arendonk JAM (2002) Genetic parame-
ters for various random regression models to describe the weight
data of pigs. J Anim Sci 80:575–582. https://doi.org/10.2527/
2002.803575x

Hurtado PX, Schnabel SK, Zaban A et al (2012) Dynamics of senes-
cence-related QTLs in potato. Euphytica 183:289–302. https://
doi.org/10.1007/s10681-011-0464-4
Kante M, Revilla P, De La Fuente M et al (2016) Stay-green QTLs in temperate elite maize. Euphytica 207:463–473. https://doi.org/10.1007/s10681-015-1575-0

Kirkpatrick M, Lofsvold D, Bulmer M (1990) Analysis of the inheritance, selection and evolution of growth trajectories. Genetics 124:979–993. https://doi.org/10.1093/genetics/124.4.979

Kranis A, Su G, Sorensen D, Woolliams JA (2007) The application of random regression models in the genetic analysis of monthly egg production in turkeys and a comparison with alternative longitudinal models. Poult Sci 86:470–475. https://doi.org/10.1093/ps/86.3.470

Kwiatkowska D (2008) Flowering and apical meristem growth dynamics. J Exp Bot 59:187–201. https://doi.org/10.1093/jxb/erm290

Laidig F, Piepho H-P, Drobek T, Meyer U (2014) Genetic and non-genetic long-term trends of 12 different crops in German official variety performance trials and on-farm yield trends. Theor Appl Genet 127:2599–2617. https://doi.org/10.1007/s00122-014-2402-z

Lê S, Josse J, Husson F (2008) FactoMineR: a package for multivariate analysis. J Stat Softw 25:1–18. https://doi.org/10.18637/jss.v025.i01

Leng G, Hall J (2019) Crop yield sensitivity of global major agricultural countries to droughts and the projected changes in the future. Sci Total Environ 654:811–821. https://doi.org/10.1016/j.scitotenv.2018.10.434

Lesk C, Rowhani P, Ramankutty N (2016) Influence of extreme weather disasters on global crop production. Nature 529:84–87. https://doi.org/10.1038/nature16467

Li Q, Huang Z, Xu M et al (2010) Functional mapping of genotype-environment interactions for soybean growth by a semiparametric approach. Plant Methods 6:13. https://doi.org/10.1186/1746-4811-6-13

Li C, Sun B, Li Y et al (2016a) Numerous genetic loci identified for drought tolerance in the maize nested association mapping populations. BMC Genom 17:894. https://doi.org/10.1186/s12864-016-3170-8

Li D, Wang X, Zhang X et al (2016b) The genetic architecture of leaf number and its genetic relationship to flowering time in maize. New Phytol 210:256–268. https://doi.org/10.1111/nph.13765

Li Y, Li C, Bradbury PJ et al (2016c) Identification of genetic variants associated with maize flowering time using an extremely large multi-genetic background population. Plant J 86:391–402. https://doi.org/10.1111/tpj.13174

Lipka AE, Tian F, Wang Q et al (2012) GAPIT: genome association and prediction integrated tool. Bioinformatics 28:2397–2399. https://doi.org/10.1093/bioinformatics/bts444

Lippert C, Listgarten J, Liu Y et al (2011) FaST linear mixed models for genome-wide association studies. Nat Methods 8:833–835. https://doi.org/10.1038/nmeth.1681

Lobell DB, Schlenker W, Costa-Roberts J (2011) Climate trends and global crop production since 1980. Science 333:616–620. https://doi.org/10.1126/science.1204531

Ma C-X, Casella G, Wu R (2002) Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics 161:1751–1762. https://doi.org/10.1093/genetics/161.4.1751

Malosetti M, Visser RGF, Celis-Gamboa C, van Eeuwijk FA (2006) QTL methodology for response curves on the basis of non-linear mixed models, with an illustration to senescence in potato. Theor Appl Genet 113:288–300. https://doi.org/10.1007/s00122-006-0294-2

Messina CD, Technow F, Tang T et al (2018) Leveraging biological insight and environmental variation to improve phenotypic prediction: integrating crop growth models (CGM) with whole

genome prediction (WGP). Eur J Agron 100:151–162. https://doi.org/10.1016/j.eja.2018.01.007

Messina CD, Gho C, Hammer GL, Cooper M (2023) Two decades of harnessing standing genetic variation for physiological traits to improve drought tolerance in maize (Zea mays L.). J Exp Bot. https://doi.org/10.1093/jxb/erad231

Messmer R, Fracheboud Y, Bänziger M et al (2011) Drought stress and tropical maize: QTLs for leaf greenness, plant senescence, and root capacitance. Field Crops Res 124:93–103. https://doi.org/10.1016/j.fcr.2011.06.010

Millet EJ, Welcker C, Kruijer W et al (2016) Genome-wide analysis of yield in europe: allelic effects vary with drought and heat scenarios. Plant Physiol 172:749–764. https://doi.org/10.1104/pp.16.00621

Millet EJ, Kruijer W, Coupel-Ledru A et al (2019) Genomic prediction of maize yield across European environmental conditions. Nat Genet 51:952–956. https://doi.org/10.1038/s41588-019-0414-y

Monneveux P, Sánchez C, Beck D, Edmeades GO (2006) Drought tolerance improvement in tropical maize source populations: evidence of progress. Crop Sci 46:180–191. https://doi.org/10.2135/cropsci2005.04-0034

Monteith JL (1977) Climate and the efficiency of crop production in Britain. Philos Trans R Soc Lond B Biol Sci 281:277–294

Moreira FF, Oliveira HR, Volenec JJ et al (2020) Integrating high-throughput phenotyping and statistical genomic methods to genetically improve longitudinal traits in crops. Front Plant Sci 11:681. https://doi.org/10.3389/fpls.2020.00681

Ning C, Kang H, Zhou L et al (2017) Performance gains in genome-wide association studies for longitudinal traits via modeling time-varied effects. Sci Rep 7:590. https://doi.org/10.1038/s41598-017-00638-2

Ning C, Wang D, Zheng X et al (2018) Eigen decomposition expedites longitudinal genome-wide association studies for milk production traits in Chinese Holstein. Genet Sel Evol 50:12. https://doi.org/10.1186/s12711-018-0383-0

O'Reilly PF, Hoggart CJ, Pomyen Y et al (2012) MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS. PLoS ONE 7:e34861. https://doi.org/10.1371/journal.pone.0034861

Pan Q, Xu Y, Li K et al (2017) The genetic basis of plant architecture in 10 maize recombinant inbred line populations. Plant Physiol 175:858–873. https://doi.org/10.1104/pp.17.00709

Passioura J (1977) Grain yield, harvest index, and water use of wheat. J Aust Inst Agric Sci 43:117–120

Peiffer JA, Romay MC, Gore MA et al (2014) The genetic architecture of maize height. Genetics 196:1337–1356. https://doi.org/10.1534/genetics.113.159152

Porter HF, O'Reilly PF (2017) Multivariate simulation framework reveals performance of multi-trait GWAS methods. Sci Rep. https://doi.org/10.1038/srep38837

Pörtner HO, Roberts DC, Adams H et al (2022) Climate change 2022: impacts, adaptation and vulnerability. Cambridge University Press

R Core Team (2017) R: A language and environment for statistical computing. R foundation for statistical computing. https://www.R-project.org, Vienna

Ray DK, Ramankutty N, Mueller ND et al (2012) Recent patterns of crop yield growth and stagnation. Nat Commun 3:1293. https://doi.org/10.1038/ncomms2296

Ray DK, West PC, Clark M et al (2019) Climate change has likely already affected global food production. PLoS ONE 14: https://doi.org/10.1371/journal.pone.0217148

Reymond M, Muller B, Leonardi A et al (2003) Combining quantitative trait loci analysis and an ecophysiological model to analyze the genetic variability of the responses of maize leaf growth
Kante M, Revilla P, De La Fuente M et al (2016) Stay-green QTLs in temperate elite maize. Euphytica 207:463–473. https://doi.org/10.1007/s10681-015-1575-0

Kirkpatrick M, Lofsvold D, Bulmer M (1990) Analysis of the inheritance, selection and evolution of growth trajectories. Genetics 124:979–993. https://doi.org/10.1093/genetics/124.4.979

Kranis A, Su G, Sorensen D, Woolliams JA (2007) The application of random regression models in the genetic analysis of monthly egg production in turkeys and a comparison with alternative longitudinal models. Poult Sci 86:470–475. https://doi.org/10.1093/ps/86.3.470

Kwiatkowska D (2008) Flowering and apical meristem growth dynamics. J Exp Bot 59:187–201. https://doi.org/10.1093/jxb/erm290

Laidig F, Piepho H-P, Drobek T, Meyer U (2014) Genetic and non-genetic long-term trends of 12 different crops in German official variety performance trials and on-farm yield trends. Theor Appl Genet 127:2599–2617. https://doi.org/10.1007/s00122-014-2402-z

Lê S, Josse J, Husson F (2008) FactoMineR: a package for multivariate analysis. J Stat Softw 25:1–18. https://doi.org/10.18637/jss.v025.i01

Leng G, Hall J (2019) Crop yield sensitivity of global major agricultural countries to droughts and the projected changes in the future. Sci Total Environ 654:811–821. https://doi.org/10.1016/j.scitotenv.2018.10.434

Lesk C, Rowhani P, Ramankutty N (2016) Influence of extreme weather disasters on global crop production. Nature 529:84–87. https://doi.org/10.1038/nature16467

Li Q, Huang Z, Xu M et al (2010) Functional mapping of genotype-environment interactions for soybean growth by a semiparametric approach. Plant Methods 6:13. https://doi.org/10.1186/1746-4811-6-13

Li C, Sun B, Li Y et al (2016a) Numerous genetic loci identified for drought tolerance in the maize nested association mapping populations. BMC Genom 17:894. https://doi.org/10.1186/s12864-016-3170-8

Li D, Wang X, Zhang X et al (2016b) The genetic architecture of leaf number and its genetic relationship to flowering time in maize. New Phytol 210:256–268. https://doi.org/10.1111/nph.13765

Li Y, Li C, Bradbury PJ et al (2016c) Identification of genetic variants associated with maize flowering time using an extremely large multi-genetic background population. Plant J 86:391–402. https://doi.org/10.1111/tpj.13174

Lipka AE, Tian F, Wang Q et al (2012) GAPIT: genome association and prediction integrated tool. Bioinformatics 28:2397–2399. https://doi.org/10.1093/bioinformatics/bts444

Lippert C, Listgarten J, Liu Y et al (2011) FaST linear mixed models for genome-wide association studies. Nat Methods 8:833–835. https://doi.org/10.1038/nmeth.1681

Lobell DB, Schlenker W, Costa-Roberts J (2011) Climate trends and global crop production since 1980. Science 333:616–620. https://doi.org/10.1126/science.1204531

Ma C-X, Casella G, Wu R (2002) Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics 161:1751–1762. https://doi.org/10.1093/genetics/161.4.1751

Malosetti M, Visser RGF, Celis-Gamboa C, van Eeuwijk FA (2006) QTL methodology for response curves on the basis of non-linear mixed models, with an illustration to senescence in potato. Theor Appl Genet 113:288–300. https://doi.org/10.1007/s00122-006-0294-2

Messina CD, Technow F, Tang T et al (2018) Leveraging biological insight and environmental variation to improve phenotypic prediction: integrating crop growth models (CGM) with whole

genome prediction (WGP). Eur J Agron 100:151–162. https://doi.org/10.1016/j.eja.2018.01.007

Messina CD, Gho C, Hammer GL, Cooper M (2023) Two decades of harnessing standing genetic variation for physiological traits to improve drought tolerance in maize (Zea mays L.). J Exp Bot. https://doi.org/10.1093/jxb/erad231

Messmer R, Fracheboud Y, Bänziger M et al (2011) Drought stress and tropical maize: QTLs for leaf greenness, plant senescence, and root capacitance. Field Crops Res 124:93–103. https://doi.org/10.1016/j.fcr.2011.06.010

Millet EJ, Welcker C, Kruijer W et al (2016) Genome-wide analysis of yield in europe: allelic effects vary with drought and heat scenarios. Plant Physiol 172:749–764. https://doi.org/10.1104/pp.16.00621

Millet EJ, Kruijer W, Coupel-Ledru A et al (2019) Genomic prediction of maize yield across European environmental conditions. Nat Genet 51:952–956. https://doi.org/10.1038/s41588-019-0414-y

Monneveux P, Sánchez C, Beck D, Edmeades GO (2006) Drought tolerance improvement in tropical maize source populations: evidence of progress. Crop Sci 46:180–191. https://doi.org/10.2135/cropsci2005.04-0034

Monteith JL (1977) Climate and the efficiency of crop production in Britain. Philos Trans R Soc Lond B Biol Sci 281:277–294

Moreira FF, Oliveira HR, Volenec JJ et al (2020) Integrating high-throughput phenotyping and statistical genomic methods to genetically improve longitudinal traits in crops. Front Plant Sci 11:681. https://doi.org/10.3389/fpls.2020.00681

Ning C, Kang H, Zhou L et al (2017) Performance gains in genome-wide association studies for longitudinal traits via modeling time-varied effects. Sci Rep 7:590. https://doi.org/10.1038/s41598-017-00638-2

Ning C, Wang D, Zheng X et al (2018) Eigen decomposition expedites longitudinal genome-wide association studies for milk production traits in Chinese Holstein. Genet Sel Evol 50:12. https://doi.org/10.1186/s12711-018-0383-0

O'Reilly PF, Hoggart CJ, Pomyen Y et al (2012) MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS. PLoS ONE 7:e34861. https://doi.org/10.1371/journal.pone.0034861

Pan Q, Xu Y, Li K et al (2017) The genetic basis of plant architecture in 10 maize recombinant inbred line populations. Plant Physiol 175:858–873. https://doi.org/10.1104/pp.17.00709

Passioura J (1977) Grain yield, harvest index, and water use of wheat. J Aust Inst Agric Sci 43:117–120

Peiffer JA, Romay MC, Gore MA et al (2014) The genetic architecture of maize height. Genetics 196:1337–1356. https://doi.org/10.1534/genetics.113.159152

Porter HF, O'Reilly PF (2017) Multivariate simulation framework reveals performance of multi-trait GWAS methods. Sci Rep. https://doi.org/10.1038/srep38837

Pörtner HO, Roberts DC, Adams H et al (2022) Climate change 2022: impacts, adaptation and vulnerability. Cambridge University Press

R Core Team (2017) R: A language and environment for statistical computing. R foundation for statistical computing. https://www.R-project.org, Vienna

Ray DK, Ramankutty N, Mueller ND et al (2012) Recent patterns of crop yield growth and stagnation. Nat Commun 3:1293. https://doi.org/10.1038/ncomms2296

Ray DK, West PC, Clark M et al (2019) Climate change has likely already affected global food production. PLoS ONE 14: https://doi.org/10.1371/journal.pone.0217148

Reymond M, Muller B, Leonardi A et al (2003) Combining quantitative trait loci analysis and an ecophysiological model to analyze the genetic variability of the responses of maize leaf growth
to temperature and water deficit. Plant Physiol 131:664-675.
https://doi.org/10.1104/pp.013839

Ribaut J-M, Ragot M (2007) Marker-assisted selection to improve drought adaptation in maize: the backcross approach, perspectives, limitations, and alternatives. J Exp Bot 58:351-360. https://doi.org/10.1093/jxb/erl214

Rincent R, Moreau L, Monod H et al (2014) Recovering power in association mapping panels with variable levels of linkage disequilibrium. Genetics 197:375-387. https://doi.org/10.1534/genetics.113.159731

Rutkoski J, Poland J, Mondal S et al (2016) Canopy temperature and vegetation indices from high-throughput phenotyping improve accuracy of pedigree and genomic selection for grain yield in wheat. G3 Genes Genom Genet 6:2799-2808. https://doi.org/10.1534/g3.116.032888

Sandhu KS, Mihalyov PD, Lewien MJ et al (2021) Combining genomic and phenomic information for predicting grain protein content and grain yield in spring wheat. Front Plant Sci 12:170. https://doi.org/10.3389/fpls.2021.613300

Sluis S, Posthuma D, Dolan CV (2013) TATES: Efficient Multivariate Genotype-Phenotype Analysis for Genome-Wide Association Studies. PLOS Genet 9:e1003235. https://doi.org/10.1371/journal.pgen.1003235

Sodini SM, Kemper KE, Wray NR, Trzaskowski M (2018) Comparison of genotypic and phenotypic correlations: cheverud's conjecture in humans. Genetics 209:941-948. https://doi.org/10.1534/genetics.117.300630

Sonderegger D (2022) SiZer: significant zero crossings

Spindel J, Begum H, Akdemir D et al (2015) Genomic selection and association mapping in rice (Oryza sativa): effect of trait genetic architecture, training population composition, marker number and statistical model on accuracy of rice genomic selection in elite. Trop Rice Breed Lines PLOS Genet 11:e1004982. https://doi.org/10.1371/journal.pgen.1004982

Spindel JE, Begum H, Akdemir D et al (2016) Genome-wide prediction models that incorporate *de novo* GWAS are a powerful new tool for tropical rice improvement. Heredity 116:395-408. https://doi.org/10.1038/hdy.2015.113

Spindel JE, Dahlberg J, Colgan M et al (2018) Association mapping by aerial drone reveals 213 genetic associations for Sorghum bicolor biomass traits under drought. BMC Genom. https://doi.org/10.1186/s12864-018-5055-5

Stephens M (2013) A unified framework for association analysis with multiple related phenotypes. PLoS ONE 8:e65245. https://doi.org/10.1371/journal.pone.0065245

Storey JD (2002) A direct approach to false discovery rates. J R Stat Soc Ser B Stat Methodol 64:479-498. https://doi.org/10.1111/1467-9868.00346

Storey JD, Bass AJ, Dabney A, Robinson D (2021) qvalue: Q-value estimation for false discovery rate control

Sun J, Rutkoski JE, Poland JA et al (2017) Multitrait, random regression, or simple repeatability model in high-throughput phenotyping data improve genomic prediction for wheat grain yield. Plant Genome 10:2016-2111. https://doi.org/10.3835/plantgenome2016.11.0111

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. https://doi.org/10.1093/jxb/err269

Technow F, Messina CD, Totir LR, Cooper M (2015) Integrating crop growth models with whole genome prediction through approximate bayesian computation. PLoS ONE 10:e0130855. https://doi.org/10.1371/journal.pone.0130855

Thomas H, Howarth CJ (2000) Five ways to stay green. J Exp Bot 51:329-337. https://doi.org/10.1093/jexbot/51.suppl_1.329

Thompson AM, Yu J, Timmermans MCP et al (2015) Diversity of maize shoot apical meristem architecture and its relationship to

plant morphology. G3 Genes Genomes Genet 5:819-827. https://doi.org/10.1534/g3.115.017541

Tollenaar M, Lee E (2006) Physiological dissection of grain yield in maize by examining genetic improvement and heterosis. Maydica 51:399

Touzy G, Rincent R, Bogard M et al (2019) Using environmental clustering to identify specific drought tolerance QTLs in bread wheat (T. aestivum L.). Theor Appl Genet 132:2859-2880. https://doi.org/10.1007/s00122-019-03393-2

Trachsel S, Sun D, SanVicente FM et al (2016) Identification of QTL for early vigor and stay-green conferring tolerance to drought in two connected advanced backcross populations in tropical maize (Zea mays L.). PLoS ONE 11:e0149636. https://doi.org/10.1371/journal.pone.0149636

Trachsel S, Burgueno J, Suarez EA et al (2017) Interrelations among early vigor, flowering time, physiological maturity, and grain yield in tropical maize (Zea mays L.) under multiple abiotic stresses. Crop Sci 57:229-242. https://doi.org/10.2135/cropsci2016.06.0562

Unterseer S, Bauer E, Haberer G et al (2014) A powerful tool for genome analysis in maize: development and evaluation of the high density 600 k SNP genotyping array. BMC Genom 15:823. https://doi.org/10.1186/1471-2164-15-823

van Eeuwijk FA, Bustos-Korts D, Millet EJ et al (2019) Modelling strategies for assessing and increasing the effectiveness of new phenotyping techniques in plant breeding. Plant Sci 282:23-39. https://doi.org/10.1016/j.plantsci.2018.06.018

Verbeke G, Molenberghs G (2009) Linear mixed models for longitudinal data. Springer, Cham. https://doi.org/10.1007/b98969

Wang Y, Xu J, Deng D et al (2016) A comprehensive meta-analysis of plant morphology, yield, stay-green, and virus disease resistance QTL in maize (Zea mays L.). Planta 243:459-471. https://doi.org/10.1007/s00425-015-2419-9

Wang H, Liang Q, Li K et al (2017) QTL analysis of ear leaf traits in maize (Zea mays L.) under different planting densities. Crop J 5:387-395. https://doi.org/10.1016/j.cj.2017.05.001

Wang B, Zhu Y, Zhu J et al (2018) Identification and fine-mapping of a major maize leaf width QTL in a re-sequenced large recombinant inbred lines population. Front Plant Sci. https://doi.org/10.3389/fpls.2018.00101

Webber H, Ewert F, Olesen JE et al (2018) Diverging importance of drought stress for maize and winter wheat in Europe. Nat Commun 9:4249. https://doi.org/10.1038/s41467-018-06525-2

Wolfe DW, Henderson DW, Hsiao TC, Alvino A (1988) Interactive water and nitrogen effects on senescence of maize. I. leaf area duration, nitrogen distribution, and yield. Agron J 80:859-864. https://doi.org/10.2134/agronj1988.000219620080000600004x

Wu R, Lin M (2006) Functional mapping—how to map and study the genetic architecture of dynamic complex traits. Nat Rev Genet 7:229-237. https://doi.org/10.1038/nrg1804

Würschum T, Liu W, Busemeyer L et al (2014) Mapping dynamic QTL for plant height in triticale. BMC Genet 15:59. https://doi.org/10.1186/1471-2156-15-59

Yang N, Lu Y, Yang X et al (2014) Genome wide association studies using a new nonparametric model reveal the genetic architecture of 17 agronomic traits in an enlarged maize association panel. PLOS Genet 10:e1004573. https://doi.org/10.1371/journal.pgen.1004573

Yang Z, Li X, Zhang N et al (2017) Mapping and validation of the quantitative trait loci for leaf stay-green-associated parameters in maize. Plant Breed 136:188-196. https://doi.org/10.1111/pbr.12451

Yang Y, Xu W, Hou P et al (2019) Improving maize grain yield by matching maize growth and solar radiation. Sci Rep 9:3635. https://doi.org/10.1038/s41598-019-40081-z