{'fraction_non_alphanumeric': 0.036324486261292886, 'fraction_numerical': 0.02754762907831297, 'mean_word_length': 4.799782815257228, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantifying senescence in bread wheat using multispectral imaging from an unmanned aerial vehicle and QTL mapping

Abstract Environmental stresses from climate change can alter source–sink relations during plant maturation, leading to premature senescence and decreased yields. Elucidating the genetic control of natural variations for senescence in wheat (Triticum aestivum) can be accelerated using recent developments in unmanned aerial vehicle (UAV)-based imaging techniques. Here, we describe the use of UAVs to quantify senescence in wheat using vegetative indices (VIs) derived from multispectral images. We detected senescence with high heritability, as well as its impact on grain yield (GY), in a doubled-haploid population and parent cultivars at various growth time points (TPs) after anthesis in the field. Selecting for slow senescence using a combination of different UAV-based VIs was more effective than using a single ground-based vegetation index. We identified 28 quantitative trait loci (QTL) for vegetative growth, senescence, and GY using a 660K single-nucleotide polymorphism array. Seventeen of these new QTL for VIs from UAV-based multispectral imaging were mapped on chromosomes 2B, 3A, 3D, 5A, 5D, 5B, and 6D; these QTL have not been reported previously using conventional phenotyping methods. This integrated approach allowed us to identify an important, previously unreported, senescence-related locus on chromosome 5D that showed high phenotypic variation (up to 18.1%) for all UAV-based VIs at all TPs during grain filling. This QTL was validated for slow senescence by developing kompetitive allele-specific PCR markers in a natural population. Our results suggest that UAV-based high-throughput phenotyping is advantageous for temporal assessment of the genetics underlying for senescence in wheat.


Introduction
Consumed by half the world, wheat (Triticum aestivum) is a major contributor to food security, but its production needs to increase by 3 billion tons annually by 2050 to meet global needs (Tester and Langridge, 2010). To achieve sustainable wheat production as climate change increase the severity of stresses on plants, we need to integrate new genomics and phenomics technologies in studies to more deeply dissect the mechanisms of plant stress responses to various extremes . Breeding based on such crossdisciplinary knowledge that links important physiological and biochemical interactions to new genetic information will better enable us to overcome agricultural challenges . The physiochemical, morphological, and molecular changes in plants in response to climatic fluctuations such as drought and heat (Zhu, 2016;Sade et al., 2017) can alter source-sink relationships, which affect growth and yield of plant. Two mechanisms are mainly involved in these alterations: (1) premature senescence/yellowing and (2) inhibition of photosynthesis activity by producing assimilates in source organs. This helps to decrease consumption of resources within the source organs and increase their mobilization toward sink tissues (Albacete et al., 2014;Yolcu et al., 2017).
Although senescence is a genetically programmed system of plants for dynamic accumulation of nutrients to sink tissues, a "stress-induced senescence" can also be activated by external stress stimuli to tradeoff the losses in sink tissues under stress conditions (Schippers et al., 2015;Sade et al., 2017). Under normal conditions, the coordinated breakdown of chlorophyll is an integral process of senescence that is developmentally planned to facilitate the remobilization of nutrients from senescing tissues to newly growing organs/ tissues such as young leaves (Lim et al., 2007). Elucidating the progress of senescence by temporal remote sensing of green biomass and chlorophyll degradation is of key importance in developing climate-resilient genotypes (Christopher et al., 2016). However, such selection for slow senescence has been difficult and error prone in the field.
Adapting field-based advanced remote sensing technology to phenotyping physiological and biochemical traits beyond conventional traits that are assessed by eye is the future of crop breeding (Araus and Cairns, 2014). Physiological attributes such as canopy chlorophyll content (CCC) and photosynthesis and the impact of nutrients, water and heat stress on plant growth Yang et al., 2020) can only be quantified nondestructively by measuring variations in the spectral reflectance . For example, when growing conditions are optimal, healthy plants look greener than those under stress because the absorbance maximum for chlorophyll is in the red wavelength region and green wavelengths are reflected (Hatfield et al., 2008). Handheld active sensors can estimate chlorophyll levels and photosynthetic rates (Pleban et al., 2020), wheat physiology can be assessed using traits derived from light reflectance such as normalized difference vegetation index (NDVI) to select for genotypes tolerant to abiotic stress (Christopher et al., 2016). Operational and resolution limitations of sensors for capturing a wide range of light bands when light reflection fluctuates, however, contribute to the difficulty in resolving the complexity of senescence and underlying genetics.
The integration of new techniques for image (RGB (Red-Green-Blue), X-ray or hyper/multispectral)-based data acquisition of traits can improve the accuracy of plant phenotyping and accelerate the molecular breeding process by increasing the rate of gene discovery (Yang et al. 2014;Ruckelshausen and Busemeyer, 2015;Campbell et al. 2017;Zhou et al. 2019). Many computer vision-based tools and imaging from various types of sensors have become vital for phenotyping with high accuracy and throughput (Mochida et al., 2018). By converting these digital measurements into useful biological knowledge (crop traits), scientists can have "smart" eyes on disease, photosynthesis, chlorophyll status, senescence, and other physiochemical properties of crops (Tardieu et al., 2017;Su et al., 2019;Rasheed et al., 2020). Recent deployment of unmanned aerial vehicle (UAV) platforms for aerial surveillance of crops using hyper/multispectral imaging has increased the access and capability of scientists to cover large number of experimental trials at multiple times during growth over a limited time Rasheed et al., 2020). UAV-based multispectral traits derived from the visible and beyond-visible range of the light spectrum, such as NDVI, green and red edge chlorophyll indices, and normalized difference red edge have been used as smart indicators to differentiate senescence rates among genotypes exposed to drought and heat stresses (Duan et al. 2017;Hassan et al., 2018). However, no effort has been made to use this knowledge in the genetic dissection of complex traits such as senescence. Therefore, linking UAVbased VIs with high-density single-nucleotide polymorphisms (SNPs) will open a new avenue in quantitative genetic studies to discover new loci for breeding.
Some genes including NAC, WRKY, MYB, AP2/EREBP, and bZIP have been reported to play significant roles in the regulation of age-induced senescence in many plant species such as Arabidopsis thaliana (Woo et al., 2013), wheat (Gregersen and Holm, 2007), and rice (Liu et al., 2008). Interestingly, these same genes are involved in stress tolerance. Most attempts to identify genes that control variation in a target trait have used data from a single time point (TP), usually at maturity or a fixed number of days from planting. Precision phenotyping to monitor plant growth over time remains a bottleneck (Furbank and Tester, 2011). Advances in UAV and computer vision are overcoming barriers to score traits at multiple times during development (Salas Fernandez et al., 2017;York, 2018). To the best of our knowledge, no studies have yet reported on integration of VIs derived from UAV-based multispectral pixels to identify loci controlling dynamic senescence in wheat. Therefore, the aims of this study were to (1) digitally quantify senescence from heading (H) to maturation using UAV-based VIs, (2) identify highly heritable predictive traits to evaluate senescence, and (3) identify loci that contribute to variations in wheat senescence.

Results
Adaptation of field-based UAV system for multispectral scan of wheat physiology We previously optimized the UAV platform for highthroughput multispectral phenotyping of the wheat canopy and predicted grain yield (GY) by monitoring seasonal growth. The accuracy of the estimates of the UAV-based traits was highly consistent with values from conventional measurements (Hassan et al., 2018;2019a, 2019b. The multispectral scans by the UAV system could thus be used to quantify the wheat senescence and thus phenotype wheat doubled-haploid (DH) populations to identify loci controlling senescence (see Materials and methods). The average maximum focal length for each sensor was 1,445.20 pixels (5.40 mm) to capture four multispectral images of light reflection, that is, near-infrared (NIR), red, green, and red-edge bands of varying wavelengths, with optimized ground sampling distance (GSD) of 2.23 cm (i.e. 1 pixel indicated 2.23 cm distance at ground; Figure 1D). Dense point-clouds for making 2D mosaic maps, over five images were overlapped for each pixel with an average reprojection error of 0.21 pixels across the whole data collection campaign ( Figure 1E). In total, 240 GB data in the form of raw multispectral images were collected through aerial surveillance of 2,400 plots during the whole phenotyping campaign consisting of 2 years and sites within 5 h flight time (Supplemental Table S1). This raw data were converted into useful pixel information through machine learning to calculate UAV-based VIs data (400 kb) to estimate the physiological status of all five growth TPs (GTPs) (see process in "Materials and methods"). The combination of multispectral pixels from different light reflectance gives a better overall estimate of the physiological status of plants. The red band is reportedly more reliable at stages before maturation due to saturation issues in detecting high chlorophyll level after canopy closure (Hatfield et al., 2008). The NIR band is best for detecting a wide range of variations in green biomass and chlorophyll, and reflectance in the green and red-edge bands ranges is also sensitive to the whole range of variations in chlorophyll and green biomass (Hatfield et al., 2008). Senescence is a complex quantitative trait, rapid, precise temporal phenotyping of each individual beyond typical physiological measurements was possible because phenotypic correlations were high (up to r = 0.82) between UAV-based VIs and GY compared to ground NDVI measurements ( Figure 1I).

Phenotypic variations in UAV-based VIs, ground NDVI, and GY
All five UAV-based VIs derived from the pixels of the 2D multispectral images obtained by aerial surveillance of field trials and the ground NDVI from the ground GreenSeeker sensor were normally distributed at all GTPs. Significant variations (at P 5 0.05) among genotypes across environments and genotype Â environment (G Â E) factor were observed for five UAV-based VIs, senescence, ground NDVI at H, anthesis (Ans), early (TP90 at 10% senescence), mid (TP50 at 50% senescence), and late (TP10 at 90% senescence concluded) grain filling (GF) GTPs (Supplemental Table S2). Genotypes also varied significantly (P 5 0.05) for GY in both environments, and the G Â E interaction was also significant. The mean broad-sense heritabilities for UAV-based VIs and GY at each GTP ranged from 0.71 to 0.95 at both locations. Heritabilities for the ground NDVI data set ranged from 0.73 to 0.95 at various GTPs ( Figure 2; Supplemental Table S3). As expected, variances among TPs increased postAns for all traits (Supplemental Table S2). Two UAVbased VIs, CCC index (CCCI) and red NDVI (RNDVI), were selected to visualize the physiological status of wheat canopy using images generated by different combinations of multispectral pixels for NIR, red, green, and red-edge reflectance for four GTPs (Figure 3). Visualization results of CCCI and RNDVI also illustrated a temporal decline in canopy greenness and phenotypic variations for slow senescence in some DH lines at TP50 to TP10. The temporal pattern of the decline in UAV-based VIs was not linear in most cases, indicating slow senescence in some DH lines ( Figure 4). Several DH lines  had slowsenescence phenotypes with high GY. Multispectral pixel visualization results were in accordance with the graphical results for these DH lines. Senescence in Zhongmai 895 was delayed compared to Yangmai 16 across the postAns GTPs, and transgressive segregation was observed for all five UAVbased VIs ( Figure 5).

Relevance of UAV-based VIs with GY assessment
All UAV-based VIs were strongly correlated with GY at all GTPs in both environments (Supplemental Table S4). For example, genotypic correlations of RNDVI, CCCI, green chlorophyll index (GCI), and red-edge chlorophyll index (RECI) with GY ranged from r = 0.22-0.82 at TP90, r = 0.28-0.62 at TP50, and r = 0.22-0.69 at TP10 in E1, while r = 0.59-0.80 at TP90, r = 0.64-0.81 at TP50, and r = 0.61-0.69 TP10 in E2. Genetic correlations between ground NDVI and GY were low, ranging from r = 0.31-0.36 in E1 and from r = 0.19-0.29 in E2 at all GTPs. A similar trend was observed in phenotypic correlations of RNDVI, CCCI, GCI RECI, and CVI with GY. Correlations for E2 were higher, ranging from 0.55 to 0.79 at TP90, 0.61 to 0.79 at TP50, and 0.58 to 0.68 at TP10 from ground NDVI compared to E1. CVI had low correlations with GY, ranging from 0.12 to 0.23 in E1.
We also examined the relevance of GTPs TP90 to TP10 for predicting GY in terms of high genetic and phenotypic correlations. To further examine the relevance of UAV-based VIs to GY, we plotted the average UAV-based VIs for the 10% higher-and 10% lower-yielding genotypes against the mean-yielding genotypes for each environment (Figure 4). Clear differences were found between the highest and lowest-yielding genotypes under both environments, with the higher-yielding genotypes having high values for UAVbased VIs. The differences appear greatest in the UAV-based VIs and ground NDVI in both environments, except for CVI in E1. The smallest differences in UAV-based VIs and ground NDVI between the high-and low-yielding groups were at TP90 when the correlation between yield and UAV-based VIs was low. These results are in accordance with the genetic correlations observed between yield and UAVbased VIs ( Figure 6). UAV-based VIs of high-yielding genotypes declined slower than those of the low-yielding genotypes, but the difference in UAV-based VIs in most of the graphs decreased as senescence progressed, such that they attained full senescence at a similar time after Ans ( Figure 3). Some DH lines showed a significant delay in trait decline with high-yielding effects. It is unlikely that the yield contrast between high-and low-yielding groups resulted from differences in the date of Ans. Differences between groups in the mean period from sowing to Ans were small (ranging from 3 to 4 d) and nonsignificant for all environments (P 4 0.05).

Multispectral image-based senescence was high at TP10
Senescence derived from UAV-based VIs and ground NDVI at TP90, TP50, and TP10 was not linear for most of traits under both environments (Supplemental Figure S1). Among the DH lines, canopy senescence (CS) was slow from Ans to TP50 in many DH lines as illustrated in the CCCI and RNDVI multispectral pixels maps and graphical results ( Figure 3; Supplemental Figure S2). However, a sudden increase was observed from TP50 to TP90. For senescence to start slowly but completing at a similar time in the slowand the fast-senescence genotypes, senescence must progress faster at particular points, resulting in an increase in  overall green leaf area and higher yield. Therefore, a majority of the DH lines reached maturity at a similar time due to the trend for fast CS at TP10. The slow CS from TP90 to TP50 provided sufficient time for the mobilization of nutrients from source to sink nutrient, resulting in high GY. CS at TP50 was estimated to be higher for GCI and RECI in E1 compared to other UAV-based VIs and E2. The box plot in Figure 6A shows a significant deviation in UAV-based VIs  Assessment of senescence progress using CCCI visualization and impact on GY. A and B, CCCI was estimated for four high-yielding and four lowyielding DH lines compared with parent cultivars at TP90, TP50, and TP10. Green pixels indicate maximum VI values within the plot. C, GY comparison of slow and fast senescence DH lines and parents. Statistical significance was determined by t test: *P 5 0.05, **P 5 0.01. from their means level across the postAns GTPs especially at TP10. This kind of natural variation in CS at TP10 can be used to detect potential loci for slow senescence and stable GY.

VIs-based selection of low-senescence and highyielding genotypes
A principal component analysis (PCA)-biplot was used to examine relationships between CS calculated from UAVbased VIs and yield for selecting the low-senescence DH lines, while maintaining high yield compared to fast-senescence individuals ( Figure 6B). As anticipated, high CS was negatively correlated with GY, as discussed above, highyielding DH lines maintained higher UAV-based VIs level across the postAns GTPs compared with the low yielding. A strong negative correlation was observed between CS derived from UAV-based VIs and GY especially at TP50 and TP10 in both environments. CS measured from CCCI and CVI were negatively correlated with those from GCI, ground NDVI, RECI, and RDVI across the postAns GTPs in E1. Vectors of CS for CCCI, CVI, GCI, and RECI were in a similar direction at TP90 and TP50 but slightly opposite at TP10, for which vectors were close to ground NDVI and RNDVI under E2. As Figure 6B illustrates, CS was higher at TP10 than at TP90 and TP50; CS vectors at TP10 were plotted opposite of those for TP90 and TP50 in E2. But an exception was observed in case of CCCI and CVI in E1; all of their CS vectors were in a similar direction and had high correlations with each other.
In some cases, plant height QTL were also detected as senescence loci (Pinto et al., 2016;Christopher et al., 2018). Height QTL associated with Rht-B1 (4B) and Rht-D1 (4D) genes, were closely colocalized with QS on chromosomes 4B and 4D in all environments (Supplemental Table S5). In our DH population, Rht-B1b was contributed by Yangmai 16 and Rht-D1b by Zhongmai 895 (Hassan et al., 2019a(Hassan et al., , 2019b. QTL on chromosomes 2B, 3A, 3D, and 5D were detected as new QTL for senescence in wheat with 5.2%-18.1% of PVE. As QTL-caas.5D (315.5 Mb) showed high PVE for the senescence ( Figure 7A), two SNPs in the region of this QTL at 315.5 Mb and 320.5 Mb with 18.13% and 10.1% of PVE, respectively, were converted successfully to kompetitive allelespecific PCR (KASP) markers. To check the effectiveness of this genomic region for slow senescence in wheat, 207 accessions of a natural population were genotyped with these KASP markers. Results showed significant differences in chlorophyll levels between genotypes within full and limited irrigations based on SNPs linked with QTL-caas.5D (Figure 7b).

Discussion
To bridge the genome to phenome gap, predictive traits with high heritabilities in the field condition are urgently needed Rasheed et al., 2020). Traits derived from UAV-based multispectral imaging have been described for precise and rapid phenotyping of wheat growth dynamics (Hassan et al., 2018). The complex mechanism underlying senescence in plants involves dynamic and diverse responses that are controlled by many loci with moderate to minor effects (Sade et al., 2017;Christopher et al., 2018). Minor-effect loci can be difficult to discover using traditional assessments of senescence-related traits. UAVbased multispectral traits can accelerate the identification of both major and minor loci and advance our understanding of the genetic architecture behind the temporal responses of plants during senescence.
We extracted five different UAV-based VIs through integrating light of various wavelengths reflected from the canopy and compared with ground truth NDVI measurements. UAV-based VIs and ground NDVI were measured at vegetative GTPs (i.e. H, Ans) and senescence GTPs (i.e. TP90, TP50, and TP10) in multiple environments (i.e. two locations in 2 years) (Figure 1, A and F). This strategy allowed us to monitor the consistency of UAV-based image pixels for calculating UAV-based VIs to predict senescence and behavior (G Â E) of the genotypes in the varied environments. This approach could increase the efficiency of predicting complex physiological traits over that of traditionally limited approaches for trait acquisitions. Transgressive segregation for UAV-based VIs and yield suggests that alleles for improvements in each of these traits could be contributed by either of the parents, though they were more often contributed by Zongmai 895 (Supplemental Table S5). The heritability (H 2 ) of UAV-based VIs was generally high in all environments, highlighting high genetic variances within environments. UAV-based VIs showed higher genetic and phenotypic correlations with GY compared to ground NDVI, as in our previous study on GY prediction based on multispectral traits (Hassan et al. 2018). High correlations indicate that UAV-based VIs from image pixels can be used as secondary traits to predict GY and that variations in these traits can explain the underlying genetics. Multispectral pixel-based field visualization also accurately illustrated the transgressive segregations in DH lines at all postAns GTPs, as validated through selected high and low senescence DH lines. Our results suggest that field trials using multispectral visualization could speed up assessment of the overall status of genotypes. High genetic correlations between UAV-based VIs and GY especially CCCI, illustrated distinct average differences in UAV-based VIs of low-and high-yielding DH lines at senescence GTPs, especially at TP50 and T10 (Figures 3 and 4). Most of the DH lines started and completed senescence within a similar time span, with some exceptions. Thus, the low rate of senescence was not due to the duration between Ans and maturity but due to the slow degradation of chlorophyll at TP90 to TP50 in highyielding DH lines (Christopher et al., 2016(Christopher et al., , 2018. During this period, source to sink mobilization was the highest as evidenced by the retention of a high green area, resulting in a high yield (Senapati et al., 2018). TP50 to TP10 senescence was greatly accelerated due to internal and external factors with the least impact on GY in high-yield DH lines. In lowyield DH lines, the rate of senescence was rapid from right after Ans to TP10, resulting in lower yields due to the rapid degradation of chlorophyll (Senapati et al., 2018). The rate of senescence from TP90 to TP50 is more crucial in terms of yield losses and grouping the high-and low-yielding genotypes ( Figure 5). But TP10 is also vital for final refilling of the sink, that could also significantly increase GY through a reenforcement mechanism (Lim et al., 2007;Borrell et al., 2014;Christopher et al., 2016). For example, GY significantly differed among DH-082 and three other slow-senescence DH lines due to rapid senescence at TP10, whereas a significant difference for GY was also observed among four slow and four rapid senescence DH lines at TP10. We detected the same kind of phenotypic tradeoff in DH lines as reported by Christopher et al. (2016). Previous PCA-biplot results showed negative correlations between CS and GY at all postAns GTPs (Hassan et al., 2018), whereas all postAns GTPs are similarly involved in affecting the GY ( Figure 6B). Based on our results, multiple UAV-based VIs are more suitable for selecting slow senescence with high GY as compared to a single ground-based trait because CS through multi-angle data sets of light reflection plotted against GY instead of a traditional trait (ground NDVI) can provide a clear picture about the physiological performance of genotypes and increase the selection accuracy.
QTL that was identified consistently in both sites and years were considered credible, and QTL common to UAVbased VIs and ground NDVI indicated the higher accuracy and credibility of UAV-based data sets for QTL analysis. QTL detected across the GTPs were temporal loci, suggesting that the causative alleles have persistent effects at the various vegetative and senescence GTPs as reported in other time-series studies (Guo et al., 2018;Lyra et al., 2020). Clustering of the QTL for different UAV-based VIs (for biomass and chlorophyll) to the same chromosomal region indicates that a single or multiple genes at one locus affect multiple UAV-based VIs for plant physiology. Each UAVbased VI was equationally different, that is, they were derived from combinations or ratios of different spectral bands, but a surrogate of the same trait. For example, both RNDVI and CVI were proxy traits of green biomass, whereas CCCI, GCI, and RECI for chlorophyll content. Therefore, identifying QTL for any surrogate trait means that these QTL likely contributed to biomass or chlorophyll. Most of the QS detected in our study have been reported previously for stay-green or senescence-related traits (Vijayalakshmi et al., 2010;Pinto et al., 2016;Shi et al., 2017;Christopher et al., 2018). For example, Christopher et al. (2018) detected staygreen QTL on chromosomes 4B, 4D, 5B, and 7B based on analysis of a DH population in four environments. Similarly, Pinto et al. (2016) identified QTL on chromosomes 1B, 4B, 4D, 5B, 7A, and 7D in an RIL population in five heat-stress environments. In our results, loci on chromosomes 2B, 3A, 3D, and 5D were detected as new QS. QTL-caas.2B_2 (560.5Mb) on 2B temporally controlled CCCI and CVI at Ans and at TP90 and TP50 and colocalized with GY locus QTL-caas.2B_1 (567.5 Mb). A QS on chromosome 3D at 248.5 Mb also contributed to controlled stage-specific senescence, influencing chlorophyll and green biomass surrogates CCCI, RECI, and CVI from Ans to TP10. We found that the QS QTL-caas.5D (315.5 Mb) identified for contributing to slow senescence was not closely linked with a major vernalization response gene (Vrn-D1) previously reported on chromosome 5D (Eagles et al., 2009). Our analysis of candidate genes showed that SNPs of the 5D locus were linked with the TraesCS5D02G403900 gene, while the gene was annotated as encoding a glycosyltransferase (Supplemental Table  S6). The glycosyltransferase protein family has been reported to be involved in stress response regulation in A. thaliana (Rehman et al., 2018) and rice (Shi et al., 2020). Moreover, our validation results of 5D QTL using two KASP markers indicated that it could be used to select slow-senescence genotypes. Both QTL on 4B and 4D associated with Rht alleles were identified as QS across the GTPs as reported previously (Pinto et al., 2016;Christopher et al., 2018), whereas GY QTL was also colocalized with senescence on 4D chromosome. GY QTL on 5A and 5B were also identified near a QS. QTL-caas.6D_3 (455.6 Mb) for CCCI, RECI and RNDVI identified at postAns GTPs also colocalized with a GY QTL on chromosome 6D. In most of these cases, alleles conferring slower senescence also conferred improved yield. In those experiments where major QTL for yield are not coincident with QTL for senescence traits, genetic factors other than senescence have a greater influence on yield (Christopher et al., 2018).
In conclusion, we found that UAV-derived VIs could be used to predict senescence and GY at postAns GTPs. More importantly, most of the traits are very difficult to measure by conventional methods and now can be assessed temporally using the low-cost UAV system. Previous QTL mapping studies have been based on conventional phenotypic data at limited growth stages, but still QTL reported in our results have good overlap with QTL previously found in several stay-green and senescence-related studies (Vijayalakshmi et al., 2010;Pinto et al., 2016;Shi et al., 2017;Christopher et al., 2018). Around 75% of the QS was detected using UAVbased VIs rather than ground NDVI. Four QTL for senescence colocalized with GY, indicating that senescence influences low-and high-yielding DH lines. We found several new QTL for senescence on chromosomes 2B, 3D, 3A, and 5D that explained the high phenotypic variations (5.2%-18.1%). We also demonstrated the temporal nature of the 5D QTL control at preAns and postAns GTPs. However, KASP markers based on SNPs linked with the 5D QTL can be used to select for slow senescence genotypes in breeding programs. Combining high-throughput phenotyping technology and large-scale QTL analysis not only greatly expands our knowledge of the dynamics of wheat development but also provides a new strategy for breeders to optimize plant physiology toward ideotype breeding in wheat.

Germplasm and field trial
A DH population of 198 lines derived from Yangmai16/ Zhongmai895 cultivars was used to phenotype UAV-based multispectral traits to quantify senescence and discover loci that control natural variations in senescence. Yangmai16 is a spring wheat cultivar that was released in the Yangtze River area of Jiangsu Province, China in 2004. It is a popular cultivar because of its high-yield potential and stability under drought. The facultative wheat Zhongmai895 was released in the southern Yellow and Huai River valleys in 2012. It has a high yield potential with high GF rate, high nitrogen-use efficiency, and tolerance to drought and high-temperature fluctuations .
Field trials were conducted during two cropping seasons (2016-2017 and 2017-2018) at two sites in Henan Province, China: Xinxiang (environment 1 [E1], 35 18 0 0 00 N, 113 52 0 0 00 E) and Luohe (E2, 33 34 0 0 00 N, 114 2 0 0 00 E) ( Figure. 1A). Both sites had experienced fluctuations in temperature, rain and daylight during maturation in both years (Supplemental Figure S3), which can decrease yield in nonresilient genotypes. The experiment was designed with 600 plots at each site, with each plot (3.6 m 2 area [3 Â 1.2 m]) having one DH line in six rows at 0.20 m spacing. Randomized complete blocks were set up with three replications at both sites and in both years to minimize experimental error. Plots were equally irrigated across the trials at the same stages: tillering and right after Ans. Plots were agronomically managed according to local practices at each site.
A natural population of 207 genotypes was used to validate the KASP markers results using chlorophyll content data. These genotypes were grown under two water treatments (i.e. full and limited irrigations) with two replications at Xinxiang ( UAV platform, multispectral imaging, and pixel extraction A DJI Phantom drone (SZ DJI Technology Co., Shenzhen, China) carrying a Micasense RedEdge multispectral sensor (Micasense, Seattle, WA, USA) (https://support.micasense.com/ hc/en-us/articles/235402807-Getting-Started-With-RedEdge) was used for all multispectral imagery (Supplemental Figure 1B). User manual of RedEdge sensor can be downloaded using link given in online links section. The DJI Phantom can fly slowly and at low altitude for 16 min. The red-edge sensor consists of five sensors, one for blue, green, red, red-edge, and NIR spectral bands. A sunshine sensor and GPS device were connected to the multispectral sensor on the top of the Phantom to calibrate environmental irradiation during light reflectance measurement and geo-referencing. Band values before and after flight were standardized using a calibration board with known reflectance provided by Micasense. DJI pro version 3.6.0 (https://www.DJI.com) was used to design the flight mission. All flights were conducted at 30 m altitude at 2 m s -1 , maintaining 85% forward and side overlaps among images to generate a dense point cloud for good quality orthomosaic. Average GSD of sensor was between 2.0 and 2.5 cm.
Pix4D mapper software (version 1.4, PIX4d, Lausanne, Switzerland) was used for orthomosaic generation ( Figure 1C). Key steps for the orthomosaic generation using Pix4D mapper comprised camera alignment, georeferencing, point cloud creation, and orthomosaic generation as previously reported (Hassan et al., 2018(Hassan et al., , 2019a(Hassan et al., , 2019bFigure 1, D and E). QGIS version 3.2 (an open-source geographic information software that support viewing and editing of data) was used for image segmentation of each plot. For this, polygon shapes were generated with the specific plot ID for a particular DH line. Spectral values in the form of pixels were extracted from the segmented parts of orthomosaic TIFF images of five different bands using polygon shapes as mask through a computer vision approach in the program IDL (version 8.6, Harris Geospatial Solutions, Inc., Broomfield, CO, USA) ( Figure 1C). The plant density of each plot was well enough to ignore the data noise due to background soil. Multispectral data were collected by the UAV at five TPsthat is, H, Ans, 10 d (TP90), 17 d (TP50), and 24 d (TP10) after Ans (Figure 1f).
Five UAV-based VIs were quantified from multispectral image pixels captured from the reflectance of four spectral bands at five important GTPs ( Figure 1G): red NDVI (RNDVI), chlorophyll VI (CVI) for green biomass, CCCI (GCI, and RECI for chlorophyll content estimation. The equations for the estimation of all VIs are given in Figure 1H and Supplemental Table S1. VIs were visualized through the combinations of particular multispectral pixels of light reflectance bands using Pix4D calculator ( Figure 1H). The assessment of GY through VIs was done by correlation analysis ( Figure 1I).

Estimation of ground truth measurements
Ground NDVI was measured using a handheld GreenSeeker at five TPs that is, H, Ans, at 10 (TP90; 10% senescence), 17 (TP50; 50% senescence), and 24 (TP10; 90% senescence) days after Ans for DH lines. Agronomic and yield-related traits such as number of days to reach to H, flowering, physiological maturity, and GY were also recorded using standard procedures detailed previously (Gao et al., 2017). Chlorophyll level of natural population was assessed using SPAD-502 Plus (Konica Minolta, Japan) at mid GF (MGF) stage to validate the KASP markers. Chlorophyll was averaged from flag leaves of 10 plants of each genotype.

Estimation of CS
Estimated values for the UAV and ground-based VIs were highest at H stage compared with all other TPs; thus, the value at this stage was considered to be the maximum VI value (VI max ). The degree of senescence was estimated by subtracting the decreasing values of VIs measured at 10 (TP90; 10% senescence), 17 (TP50; 50% senescence), and 24 (TP10; 90% senescence) days after Ans. The mean senescence was estimated after each 100 degree-days since Ans using the following equation: where CS is mean CS, VI max indicates the maximum VI value (at H), GTP TP90, TP50, and TP10 indicate the VI values at 10%, 50%, and 90% senescence stages, respectively. These TPs were calculated according to thermal TPs after Ans (Lopes and Reynolds, 2012;Christopher et al., 2018). Genotypes with slow senescence were selected by plotting the UAV-based VIs derived senescence at TP90 to TP10 against the GY using a PCA-biplot analysis.

SNP genotyping and QTL mapping
The 198 DH lines and parents were genotyped using the commercially available Wheat 660K SNP array that was developed by Affymetrix and Prof. Jizeng Jia at the Institute of Crop Sciences, CAAS (Capital Bio Corporation, Beijing, China). Previously, our research group constructed a genetic map for mapping QTL for key yield-related traits. Markers with distorted segregation, no polymorphisms among parents, and missing at a rate 420% were removed in a subsequent linkage analysis. Around 10,242 markers, each representing a bin site, were selected to construct the linkage map for the DH population. The map comprises 25 linkage groups for all 21 chromosomes of the A, B, and D genomes. Inclusive composite-interval mapping was used for QTL analysis in IciMapping version 4.1 software (Meng et al., 2015). The averaged data of each trait and GTPs from 2 years at Xinxiang (E1) and Luohe (E2) separately, and best linear unbiased estimates (BLUEs) values of across environments (E1 and E2) were used for QTL detection. BLUEs values were calculated using a model explained by Alvarado et al. (2020) for a randomized complete block design (RCBD) experiment in META-R software (a multi-environment analysis tool developed at CIMMYT, Mexico City, Mexico). The genotype of Yangmai 16 was defined as A, that of Zhongmai 895 as B. Hence, alleles from Yangmai 16 reduced trait values when the additive effects were negative. Recombination frequencies were converted into map distance using the Kosambi mapping function. Locations of QTL for the traits were determined by ICIM-ADD using the same software as for the linkage analysis. The threshold for declaring the presence of a significant QTL for each trait was defined by 1,000 permutations at P 5 0.05, and the minimum LOD score of 2.5 was chosen; the walking speed was set at 1.0 cM.

KASP marker development
KASP markers for locus QTL-caas.5DL were developed based on two corresponding SNPs to validate the UAV-based VIs for genetic dissection of senescence (Supplemental Table  S5). The flanking sequences of the SNPs were used as queries in a blast search against the reference genome of wheat using IWGSC (CS Refseq version 1.0; IWGSC 2018). Chromosome-specific KASP primers were developed by alignment of homologous sequence. Allele-specific primers carrying FAM and HEX were designed with targeted SNPs at the 3 0 -end, and common reverse primers were designed with less than 200 bp of amplified sequences for chromosome-specific amplification. The detailed information about the two KASP markers is provided in Supplemental  Table S7. The KASP assay mixture consisted of 40 mL of common primer (100 mL), 16 mL of each tailed primer (100 mL), and 60 mL ddH 2 O. Each reaction mixture comprised 2.5 mL of 2ÂKASP master mixture (LGC Genomics, Hoddesdon, UK), 0.056 mL of KASP assay mixture, and 2.5 mL of DNA (30-50 ng/mL). PCR was performed in a 384-well plate with denaturation at 95 C for 15 min; 9 touchdown cycles (95 C for 20 s touchdown at 65 C initially then deceasing by 1 C per cycle for 1 min); and 32 cycles of denaturing, annealing, and extension (95 C for 10 s, 57 C for 1 min). KASP genotyping results were compared among 207 accessions of a natural population for chlorophyll level at mid to late GF (LGF) stage and original chip-based results using a t test in XLSTATE software to confirm the association of the KASP markers with the trait.

Analysis of putative candidate genes
The genes located in the physical intervals of the genomic region of identified QTL were screened based on annotations in the wheat reference genome (CS RefSeq version 1.0; IWGSC 2018), and those related to growth, development, stress resistance, and nutrient mobilization were considered as candidate genes. Gene annotations were retrieved using EnsemblPlant and EMBL-EBI (http://www.ebi.ac.uk/inter pro) databases. Gene annotations for putative proteins were done using BLAST2GO (https://www.blast2go.com/).

Statistical analyses
Phenotypic and genetic correlation matrices were calculated to evaluate the relationship between all observed parameters using META-R software (a multi-environment analysis tool developed at CIMMYT, Mexico City, Mexico; Alvarado et al., 2020). Logistic regression was done to check the relevance of UAV-based VIs for GY using a generalized linear model in XLSTAT software by Addinsoft. A mixed linear model was used to test the significance of variation among DH lines, and G Â E interaction for UAV-based VIs, ground truth NDVI and GY explained by Alvarado et al. (2020). The results were considered as significant at P 5 0.05. (2) Where Y ijk is the trait of interest, l is the overall mean effect, Loc i and Rep j is the random effect of the ith location and jth replicate, Gen k is the fixed effect of the kth genotype, Loc i and Loc i Â Gen k are the random effects of the ith environment and the G Â E interaction, respectively, and e ijk is the effect of the error associated with the ith location, jth replication and kth genotype, which is assumed to be independently and identically distributed (iid) normal with mean zero and variance r e 2 .
To confirm that phenotypic variations of traits were due to genetic diversity, we calculated broad-sense heritabilities for all traits in each and across the environments using entries as a random effect using the following equation (Sehgal et al., 2015): where r g 2 and r e 2 represent the genotypic and error variances, respectively. The term r ge 2 is the G Â E interaction variance component, e is number of environments (locations and years) and r is the number of replications for each genotype in each environment. To calculate the H 2 for each of two locations e is considered number of years. H 2 of each component at all growth stages provides an indication of the consistency of the trait in a particular environment.
PCA was used to assess the diversity among DH lines for senescence. The basic equation used for the PCA in matrix notation was where W is a matrix of coefficients that is determined by PCA and X is an adjusted data matrix consisting of n observations (rows) on p variables (columns).

Supplemental data
The following materials are available in the online version of this article. Supplemental Figure S1. Climate conditions (temperature and rain) during 2016-2017 and 2017-2018 growing seasons.
Supplemental Figure S2. CS estimated from VIs at three GTPs for two environments.
Supplemental Figure S3. Summary of QTL identified as contributing to temporal control of senescence at individual GTPs and GY.
Supplemental Table S1. Information of VIs, data size in raw and biological formats and total time for the entire field phenotyping campaign.
Supplemental Table S2. Significance test (ANOVA) results for all traits and GTPs.
Supplemental Table S3. H 2 results for all traits and growth time points (GTPs).
Supplemental Table S4. Average genetic and phenotypic correlations between GY and VIs.
Supplemental Table S5. QTL identified for vegetative growth, senescence and GY.
Supplemental Table S6. QTL with corresponded candidate genes and putative proteins.