Identifying developmental QTL alleles with favorable effect on grain yield components under late‐terminal drought in spring barley MAGIC population

Abstract Barley is the fourth most cultivated cereal worldwide, and drought is a major cause of its yield loss by negatively affecting its development. Hence, better understanding developmental mechanisms that control complex polygenic yield‐related traits under drought is essential to uncover favorable yield regulators. This study evaluated seven above‐ground yield‐related traits under well‐watered (WW) and late‐terminal drought (TD) treatment using 534 spring barley multiparent advanced generation intercross double haploid (DH) lines. The analysis of quantitative trait loci (QTL) for WW, TD, marker by treatment interaction, and drought stress tolerance identified 69, 64, 25, and 25 loci, respectively, for seven traits from which 15 loci were common for at least three traits and 17 were shared by TD and drought stress tolerance. Evaluation of allelic effects for a QTL revealed varying effect of parental alleles. Results showed prominent QTL located on major flowering time gene Ppd‐H1 with favorable effects for grain weight under TD when flowering time was not significantly affected, suggesting that this gene might be linked with increasing grain weight by ways other than timing of flowering under late‐terminal drought stress. Furthermore, a desirable novel QTL allele was identified on chromosome 5H for grain number under TD nearby sucrose transporter gene HvSUT2. The findings indicated that spring barley multiparent advanced generation intercross population can provide insights to improve yield under complex condition of drought.

important causes of barley yield loss (Jamieson et al., 1995;Rollins et al., 2013) as it negatively affects plant development and yield components such as number of grains per m À2 and thousand grain weight (Pennisi, 2008). Drought tolerance is controlled by various complex mechanisms, and therefore, improving it is a challenging task which requires better understanding plant development under drought stress (Sallam et al., 2019).
Drought can cause yield loss during all phases of plant life cycle (Salekdeh et al., 2009). However, sensitivity to drought stress varies in different stages of crop development. Early-season drought can damage yield by reducing seedling survival during the vegetative development (Lelièvre, 1981). On the other hand, late-season drought which occurs during reproductive development can have devastating effect on grain number per unit area and thousand grain weight and is more likely to occur in field; it can be referred to as terminal drought or lateterminal drought (TD) (Farooq et al., 2014;Jamieson et al., 1995;Samarah et al., 2009;Samarah & Alqudah, 2011;Shavrukov et al., 2017). Loss of grain yield components including grain weight, grain number, and ear number is an indication of yield reduction in barley. However, because number of grains per m À2 is the most important component of cereal yield Slafer, 2003;Sreenivasulu & Schnurbusch, 2012), reduction of ear number and grain number per ear is reported to be strongly connected to yield loss under drought stress (Samarah et al., 2009). Drought stress can reduce number of grains per m À2 by causing loss of spikelets, increasing floret sterility, and loss of seed set (Dolferus et al., 2011;Svobodová & Míša, 2004). After anthesis and seed set, the main negative effect of drought is by reducing the grain filling rate and duration by restricting the final number of endosperm cells or limiting the rate and duration of starch accumulation in the endosperm (Radchuk et al., 2009;Setter & Flannigan, 2001;Wallwork et al., 1998).
Gaining more insights into role of genetic mechanisms controlling plant phenology and development and their influence on grain yield components under drought is mandatory for developing high-yielding drought-tolerant barley cultivars (Cattivelli et al., 2008). Among developmental events, flowering time is a key trait which is controlled by genetic mechanisms that interact with environment (Afsharyan et al., 2020). Impact of drought stress on regulation of flowering time genes can affect grain yield components (Gol et al., 2020;Wiegmann et al., 2019). Major flowering time loci were mapped in association with grain yield components under drought stress such as Vrn-H1 gene and denso/sdw1 gene for tiller number and Vrn-H3 (HvFT1) gene for fresh weight . However, little is known about the contribution of developmental mechanisms including flowering time to yield during TD, when drought has little or no influence on timing of flowering. One study found that photoperiod gene PSEUDO-RESPONSE REGULATOR Ppd-H1 (HvPRR37) locus is linked with plant biomass and thousand grain weight when drought treatment was implemented after flowering time (Honsdorf et al., 2017). Therefore, the contribution of development and flowering time genetic loci to yield during late-terminal drought stress is not well elaborated yet.
To empower quantitative trait loci (QTL) detection, past QTL mapping efforts in barley used various approaches such as advanced-backcross population (Sayed et al., 2012), introgression lines (Honsdorf et al., 2014(Honsdorf et al., , 2017, and association panels (Wehner et al., 2015). In addition to mapping QTL associated with traits under normal and drought conditions, QTL were also detected based on stress indices as phenotype input (Honsdorf et al., 2014;Oyiga et al., 2020). Fernandez et al. (1992 proposed stress tolerance (ST) index (STI) which is calculated based on phenotyping data from normal and stress conditions to identify genotypes that show better tolerance under stress condition.
This index was used previously to directly detect QTL associated with drought tolerance (Oyiga et al., 2020). However, QTL mapping under drought condition proved complicated due to large genotype by environment interactions which resulted in detecting smaller effect QTL (Abdel-Haleem et al., 2012;Honsdorf et al., 2014;Li et al., 2013;Wehner et al., 2015). QTL mapping for drought tolerance was also challenging due to its complex nature including contribution of many genes with small effects (Honsdorf et al., 2014;Sallam et al., 2019).
Reports for using a more advanced mapping approach such as multiparental populations are rare. One study used a nested association mapping population to study traits such as dry weight, fresh weight, plant height (PLH), tiller number under control, and drought condition . Multiparent advanced generation intercross (MAGIC) strategy was designed to increase the power of QTL mapping (Bandillo et al., 2013;Cavanagh et al., 2008;King et al., 2012). In this study, we used a spring barley MAGIC population constructed from eight-way cross of seven landraces known as "founders of German barley" and one elite cultivar Barke (Sannemann et al., 2015) which has been utilized recently in a series of flowering time studies. A recent study which used this MAGIC population showed very precise recognition of known flowering time loci as a proof of concept, in addition to identifying new loci (Afsharyan et al., 2020;Mathew et al., 2018;Maurer et al., 2017;Sannemann et al., 2015).  Sannemann et al. (2015). Ragusa is a six-rowed facultative barley and other parents are two-rowed spring barley ecotypes. To study the effect of genetic mechanisms involved in plant development and flowering time on grain yield components under TD, seven above-ground yield-related traits were evaluated under WW and TD treatments. MAGIC DH lines and a set of controls consisting of eight parents and spring barley varieties were sown into 19.5 Â 25.5 cm plastic pots filled with 5.5 L of Terrasoil ® . For each treatment, four seeds per genotype were sown in each pot on April 4 (mean temperature: 11.99 C) and April 3 (mean temperature: (detecting two tillers-beginning of stem elongation) in 2012 which was 9 and 13 days before start of heading for years 2011 and 2012, respectively, for TD treatment. The treatment was imposed for 5 weeks in three stages. During the first 21 days, the water content in the pots was reduced to the permanent wilting point (15% VWC); then for 8 days, the water content was stabilized at 15% VWC; at 65 DAS, the pots under TD were rewatered slowly to 30% VWC; and finally at 73 DAS, the pots were rewatered to 40% VWC. For WW, VWC was sustained at 40% throughout the experiment (Figure 1).

| Phenotypic traits
The phenotypic values for days to heading (DHE), grain filling period (GFP), PLH, above-ground biomass (AGB), number of ears (NE), number of kernels (NK), and thousand kernel weight (TKW) were measured as described in Table 1. For evaluating ST for traits that were affected by TD, STI was calculated for each DH line in each year using the following equation: where y p is the trait for genotype under WW, y s is the trait for genotype under TD, and ȳ p is the trait mean for genotype under WW (Fernandez, 1992).

| Analysis of phenotypic data
Descriptive statistics were calculated for seven traits across 2 years by using summary function from core generic functions of R software to determine minimum, maximum, and mean. Functions std.error, var, and sd were used to calculate standard error, standard deviation, and coefficient of variation (CV), respectively, using R software (R core team, 2015). Analysis of variance (ANOVA) was performed for each trait using the restricted maximum likelihood (REML) method by PROC MIXED procedure in SAS (9.4 version, SAS Institute Inc., Cary, NC, USA) to fit the following mixed model: where Yijk is the response variable; μ is the general mean; Li is the random effect of ith DH line; Tj is the fixed effect of jth treatment; Ck is the random effect of the kth calendar year; Li Â Tj is the random effect of interaction of ith DH line with jth treatment, and εijk is the residual.
F I G U R E 1 Soil moisture content for well-watered (black) and terminal drought (red) treatments as well as daily mean temperature (orange) in 2011 and 2012. The days after sowing (DAS) in the start of reduction in irrigation, permanent wilting point, and rewatering is indicated for lateterminal drought treatment.
Variance components were estimated for each treatment by taking genotype (nonreplicated MAGIC DH lines), year, and genotype Â year interaction as random effects using REML method by PROC VARCOMP in SAS. Then, broad sense heritability (H 2 ) for each trait was estimated as where V G is the variance of genotype; V GY is the variance of genotype Â year; V E is the variance of experimental error; and y is the number of years.
Effect of row and column was evaluated separately for WW and TD treatments in R software with lm function to fit a fixed model by taking nonreplicated MAGIC DH lines, row, column, and controls as fixed effects. LSmeans was calculated for each trait using lsmeans function from "emmeans" package in R and then was used to calculate correlations among traits for each treatment. Correlation coefficients were computed by the Pearson's coefficient (r) using rcorr function from package "Hmisc" in R software. The QTL analysis was performed separately for WW, TD, marker by treatment interaction (M Â T), and ST using the PROC MIXED procedure in SAS 9.4 using the following linear model:

| Marker-trait association analysis
where Yij is the response variable; μ is the general mean; Mi is the fixed effect of ith marker; Tj is the fixed effect of the jth treatment  QTL was addressed by implementing multilocus analysis (Bauer et al., 2009;Kilpikari & Sillanpää, 2003;Sillanpää & Corander, 2002).
This approach is a forward selection procedure that in each iterative cycle inserts the most informative SNP inside the model and uses it to reanalyze the remaining SNPs. The iteration continues until no other SNP is found. P-values were calculated by F-tests and probability of false discovery rate was used to adjust P-values by incorporating the control of the QTL false discovery rate (false discovery rate value ≤.05) inside the model using PROC MULTTEST procedure. The model defined the significance of the SNPs for main effects as well as M Â T by a threshold of P-value ≤ 0.001 with 1000 permutations and false discovery rate value ≤.05 as putative QTL for the next iteration (Doerge & Churchill, 1996). The significance of QTL was validated by calculating the mean P-value of 20 rounds of a "leave20%out" crossvalidation procedure. The codes and procedure regarding QTL analysis including M Â T analysis was described in detail by Afsharyan et al. (2020). The QTL intervals were defined in the model by clustering SNPs based on their significance in the first iteration of the multilocus procedure. The confidence interval was set as 3.5 cM on both sides of the most significant SNP marker based on linkage disequilibrium of the population (Sannemann et al., 2015). Genetic variance explained by a single SNP marker (R M 2 ) was conducted as where SQ M is the sum of squares of M and SQ g is the type I sum of squares of the barley MAGIC DH lines in an ANOVA model (Von Korff et al., 2006). Additionally, total explained genetic variance by all QTL was determined.
In silico analysis was performed using the IPK barley BLAST server (Colmsee et al., 2015). Analysis of allelic variation was performed by using haplotype data as input for marker-trait association analysis as described by Afsharyan et al. (2020).

| Evaluating traits under WW and TD treatments in MAGIC DH lines
To study seven yield-related traits under WW and TD, a set of 534 spring barley MAGIC DH lines were used in a pot experiment under foil tunnel in two consecutive years. The start of reducing irrigation was shortly before flowering at 35 DAS (Figure 1). ANOVA revealed that TD effect was significant (P < 0.001) for all traits except for DHE (Table S1) and, depending on the trait, showed a reduction ranging from 2% to 46% under TD (Figure 2

| Marker-trait associations under WW and TD treatments
Genetic analysis was performed to identify the associated genetic regions under WW and TD for each trait (Figure 3). Results showed that some of genetic regions harbored QTL for various traits or treatments such as the region on chromosome 2H (18.90-28.70 cM), but F I G U R E 2 Box-whisker plots describe the variation for six traits affected by treatment in 534 spring barley multiparent advanced generation intercross population under WW and TD treatments using mean values for 2011 and 2012. Trait names and units are indicated above their respective subplot. Significant treatment effect is indicated with red asterisks with *P < 0.05, **P < 0.01, and ***P < 0.001. Change in trait mean (%) under TD treatment compared with WW is shown below the asterisks. AGB, above-ground biomass; DHE, days to heading; GFP, grain filling period; NE, number of ears; NK, number of kernels; PLH, plant height; TD, late-terminal drought; TKW, thousand kernel weight; WW, wellwatered.
others contained trait and/or treatment specific QTL such as locus on chromosome 1H (71.03 cM) for flowering time and locus on chromosome 4H (34.4-35.9 cM) for TD which was associated with NK and AGB ( Figure 4). In total for the seven traits, 69 and 64 QTL were found for WW (Table S5) and TD (  (Figure 3). The results revealed overall 25 genetic regions interacting with treatment for four traits which ranged from 10 for GFP to four for PLH (Table S7). The F I G U R E 3 Manhattan plots describe quantitative trait loci analysis of seven traits for WW treatment, TD treatment, M Â T, and drought tolerance using the spring barley multiparent advanced generation intercross population. The y-axes denote the significance of SNP markers as Àlog10 (P); the chromosomes are denoted on the x-axes. The highlighted SNP markers above the cutoff line are significant by a threshold of P ≤ 0.001 with 1000 permutations plus 20 times cross-validation.  (Table S7).
F I G U R E 4 Genetic map of quantitative trait loci (QTL) detected for traits days to heading, grain filling period, plant height, above-ground biomass, number of ears, number of kernels, and thousand kernel weight under well-watered and late-terminal drought treatments as well as stress tolerance in spring barley multiparent advanced generation intercross DH lines. Barley chromosomes are exhibited with white bars. The position for the peak SNP marker represents each QTL which is colored according to increasing (red) or decreasing (blue) effect of the minor allele of QTL. The located QTL are detailed in Tables S5, S6, and S8. The names of known genes as described for the Barke Â Morex recombinant inbred lines (RILs) by Mascher et al. (2013) are italicized in black indicating their position on chromosomes.
Chromosomal regions associated with drought tolerance were evaluated for each trait that was affected by treatment ( Figure 3). In total, 25 QTL were identified (Table S8) (Mogensen, 1992;Samarah, 2005;Sánchez-Díaz et al., 2002) and loss in seed set (Del Moral et al., 2003;González et al., 1999;Mogensen, 1992). As a result of timing of treatment, number of sterile fully developed florets did not decrease; therefore, loss in grain number could be due to grain abortion (González et al., 2003;Miralles et al., 2000). Height reduction in MAGIC DH lines might be the result of decreased gross photosynthetic rate and osmotic potential (González et al., 1999;Hopkins & Wilhelmova, 1997;Taiz & Zeiger, 2002). TD most impacted AGB which was positively correlated with PLH (González et al., 1999;Sánchez-Díaz et al., 2002) and NE. Severe impact of drought treatment on above-ground biomass was reported in previous studies Samarah et al., 2009). Reduction of NK, NE, TKW, and GFP in MAGIC DH lines was strong indication of yield loss (Samarah et al., 2009). On the other hand, positive correlation of TKW with GFP was smaller and in opposite direction with NK and NE, suggesting rise in source-sink ratios (Serrago et al., 2013) led to increase in weight of single grain which partly compensated the negative effect of shorter GFP. Due to the possible effect of pot size on root growth, further studies in field condition are required.  (Tables S5-S8). Peak marker BK_12 on chromosome 2H is gene-specific for Ppd-H1 gene (Colmsee et al., 2015) and was detected for traits GFP (WW, TD, and M Â T) and TKW (TD and ST) which were reported to be associated with Ppd-H1 region (Honsdorf et al., 2017;Maurer et al., 2016). A previous study also reported detecting the same marker for flowering time,

| QTL for WW and TD treatments and drought tolerance
which is known to be regulated by Ppd-H1 gene, as an indication of QTL mapping precision in another multiparental population (Maurer et al., 2015).  (Tables S5, S6, and S8). Another well-known locus Vrs1 (Sixrowed spike 1) which controls grain number in inflorescence (Komatsuda et al., 2007)  and the causative allele is also inherited from Ragusa which is the only six-rowed parental line (Figure 4). Vrs1 allele inherited from six-rowed barley is expected to have a stronger increasing effect for NK compared with other two-rowed parental lines as well as the opposite effect on TKW as higher grain number is expected to induce competition among single grains for assimilates (Maurer et al., 2016).  -Haleem et al., 2012;Li et al., 2013;Sallam et al., 2019;Wehner et al., 2015) that could be missed (Honsdorf et al., 2014). However, this study successfully detected regions includ- and revealed treatment interaction for PLH and AGB (Tables S5, S6, and S8).
Flowering time loci reported for this population including Vrn-H1 showed desirable effect >0.5 trait unit and were commonly detected for drought tolerance (Figure 5a and Tables S5, S6, and S8). These QTL were positioned in loci harboring genes involved in developmental mechanisms. The QTL for TKW colocated with major flowering time locus Ppd-H1. However, this region was not detected for flowering time in this study (Afsharyan et al., 2020;Sannemann et al., 2015), and effect of the most prominent QTL for flowering time which colocated with another major gene Vrn-H3 did not correspond to effect detected for Ppd-H1 locus (Figure 5b). This region was reported to be mapped for thousand grain weight under drought treatment that was imposed after flowering time (Honsdorf et al., 2017). Therefore, the results suggest Ppd-H1 might improve drought tolerance and grain weight by ways other than timing of flowering. The allele underlying this locus was inherited from parental line Ragusa and showed 6.23 g increase in grain weight (Figure 4). The QTL allele for NK, located on 5H chromosome (0.14 cM), was not reported before. This region was also mapped for NE under WW, TD, and for ST which was in line with reports that found it linked with tiller number in control (Alqudah et al., 2016) and drought stress condition . The far distal region of 5H chromosome was reported to be associated with grain size parameters (area, width, and diameter) and grain filling rate in normal field condition (Du et al., 2019). NK and NE were negatively correlated which corresponded to opposite effect of this locus on grain number and ear number (Figure 5a and  (Weschke et al., 2000). It was reported that restricted carbohydrate supply might lead to grain abortion in waterstressed maize (Zinselmeier et al., 1999). Therefore, the results suggests that sugar-related genes might contribute to controlling grain number which in barley is reported to be majorly controlled by Vrs1 gene (Komatsuda et al., 2007). Vrs1 locus is known to control spike row number as well as affecting traits related to grain size and weight (Ayoub et al., 2002;Sharma et al., 2018;Wang et al., 2016Wang et al., , 2017. This locus showed that the causative allele from Ragusa had favorable effect to increase grain number per ear by 2.50. However, the increase in grain number might have a negative influence on grain weight (Calderini et al., 2021;Maurer et al., 2016), and the impact of the trade-off should be considered on the potential yield. The findings indicated that spring barley MAGIC population can provide insights to reveal favorable QTL alleles to improve drought tolerance and grain yield under drought stress in breeding programs.

AUTHOR CONTRIBUTIONS
NPA was involved in conceptualization, methodology, statistics and QTL mapping, visualization, interpretation and discussion, writingoriginal draft, and writing-review and editing. WS was involved in conceptualization, methodology, foil tunnel experiment and phenotyping, and writing-review and editing. AB was involved in conceptualization, methodology, and writing-review and editing. JL was involved in conceptualization, methodology, supervision, project administration, funding acquisition, and writing-review and editing.
All the authors have read and approved the final manuscript.

CONFLICT OF INTEREST STATEMENT
The Authors did not report any conflict of interest.

PEER REVIEW
The peer review history for this article is available in the Supporting Information for this article.

DATA AVAILABILITY STATEMENT
Data described in this manuscript are available upon request to the author.