Spectral reflectance indices as proxies for yield potential and heat stress tolerance in spring wheat: heritability estimates and marker-trait associations

The application of spectral reflectance indices (SRIs) as proxies to screen for yield potential (YP) and heat stress (HS) is emerging in crop breeding programs. Thus, a comparison of SRIs and their associations with grain yield (GY) under YP and HS conditions is important. In this study, we assessed the usefulness of 27 SRIs for indirect selection for agronomic traits by evaluating an elite spring wheat association mapping initiative (WAMI) population comprising 287 elite lines under YP and HS conditions. Genetic and phenotypic analysis identified 11 and 9 SRIs in different developmental stages as efficient indirect selection indices for yield in YP and HS conditions, respectively. We identified enhanced vegetation index (EVI) as the common SRI associated with GYunder YP at booting, heading and late heading stages, whereas photochemical reflectance index (PRI) and normalized difference vegetation index (NDVI) were the common SRIs under booting and heading stages in HS. Genomewide association study (GWAS) using 18704 single nucleotide polymorphisms (SNPs) from Illumina iSelect 90K identified 280 and 43 marker-trait associations for efficient SRIs at different developmental stages under YP and HS, respectively. Common genomic regions for multiple SRIs were identified in 14 regions in 9 chromosomes: 1B (60–62 cM), 3A (15, 85–90, 101– 105 cM), 3B (132–134 cM), 4A (47–51 cM), 4B (71– 75 cM), 5A (43–49, 56–60, 89–93 cM), 5B (124–125 cM), 6A (80–85 cM), and 6B (57–59, 71 cM). Among them, SNPs in chromosome 5A (89–93 cM) and 6A (80–85 cM) were co-located for yield and yield related traits. Overall, this study highlights the utility of SRIs as proxies for GY under YP and HS. High heritability estimates and identification of marker-trait associations indicate that SRIs are useful tools for understanding the genetic basis of agronomic and physiological traits.


Introduction
Bread wheat (Triticum aestivum L.) is one of the staple food crops in the world providing 20% of the proteins and calories to the world population; increasing grain yield (GY) of spring wheat is therefore important for world food security [1] . Short episodes of extreme heat are expected to become more frequent, which could have adverse impacts on yield [2] . Conventional breeding based on the selection of yield per se has been very efficient in raising wheat production under irrigated conditions in the past [3] . However, the estimation of GY requires the harvest of the experimental plots, which is expensive and timeconsuming at a breeding scale [4] . The identification of proxies for GY in early stages of crop growth can increase the efficiency of selection and reduce the number of plots harvested [5] . Use of morphological and physiological traits as indirect selection criteria for GY is an alternative breeding approach [6] . In recent years, spectroradiometry has offered a quantum leap forward in the efficiency of screening physiological traits that have been typically laborious and time-consuming [7][8][9] . The prospect of future genetic improvements using spectral reflectance to identify and track physiological traits, offers crop breeding programs new opportunities to increase yield potential (YP) and to improve wheat to abiotic stresses such as heat stress (HS) [10,11] .
The spectral reflectance of a canopy is directly linked to a plant's biophysical and biochemical properties, which determine absorption and reflection patterns at specific wavelengths [12] . Pigments in leaves absorb light strongly in the visible spectrum (VIS, 380-740 nm) rather than the near-infrared region (NIR, 740-2500 nm) [10] , resulting in a different reflection of radiation in the two regions. Spectral reflectance indices (SRIs) represent the simplest way to obtain information from spectral reflectance measurements and were developed using simple mathematical formulae such as ratios or differences between reflectances observed at different wavelengths. The simple ratio (SR) index [13] and the normalized difference vegetation index (NDVI) [14] were the first SRIs to be described and are commonly used to predict green biomass and green leaf area index. Some SRIs use only reflectance at the VIS wavelentghs, such as the photochemical reflectance index (PRI) which is used to assess radiation-use efficiency [15] ; or only NIR, such as the water index (WI) that is related to water status of the canopy [16,17] . SRIs have been widely used in assessing different physiological conditions of the canopy such as the estimation of total dry matter, leaf area index, photosynthetic capacity, and pigment concentration [18][19][20][21] . GY has also been frequently estimated using SRIs measured at different growth stages of the wheat crop under irrigated and water-limited environments [22][23][24][25][26] .
SRIs may be related to complex quantitative traits controlled by multiple genetic factors. Recent studies identified the genomic regions for several SRIs under fully irrigated, rainfed, and drought-stressed environments in bread wheat [27,28] . However, little is known about the genomic regions associated with SRIs in spring wheat under HS. In the present study, a spring wheat association mapping initiative (WAMI) panel was phenotyped under YP and HS conditions. The objectives of this study were to (1) identify efficient SRIs as proxies for GY under YP and HS environments, and (2) identify the genomic regions associated with SRIs for YP and HS through genome-wide association study (GWAS).

Germplasm, genotyping, and population structure
The WAMI panel consisting of 287 elite spring wheat lines was assembled at the International Maize and Wheat Improvement Center (CIMMYT) for dissecting the genetic basis of complex traits. This is a collection of elite lines assembled from the international nurseries that CIMMYT distributes to developing countries through the International Wheat Improvement Network (IWIN). The material was genotyped using the Illumina iSelect 90K single nucleotide polymorphism (SNP) array. A total of 21871 assays showed three distinct clusters corresponding to the AA, AB, and BB genotypes expected for a biallelic SNP. Of the remaining assays with poorer cluster separation, manual clustering was applied. Overall, 36133 of the 81587 functional SNPs visually revealed polymorphism among the WAMI population. We were able to locate 28614 (35.1%) of the SNPs on the published genetic map [29] . After removing SNPs with minor allele frequency < 5% and missing value > 10%, 18704 SNPs were used for GWAS. The population structure of the panel was associated with pedigree of the lines. A detailed description of the population structure, linkage disequilibrium, and minor allele frequency can be found in earlier studies of the WAMI population [30][31][32] .

Field experiment and agronomic practices
Field trials were conducted during the 2015-2016 crop season under field conditions at Campo Experimental Norman E. Borlaug (CENEB), Ciudad Obregon, Yaqui Valley, Mexico (27°25′ N, 109°54′ W, 38masl). The WAMI population was phenotyped in a randomized lattice design with two replications under YP and HS. The YP environment was simulated by planting in late November and full irrigation (total water supply > 700 mm) was provided. HS was applied by late sowing (late February) with full irrigation (total water supply > 700 mm) to avoid the effect of drought; a typical method for applying HS at CENEB that has been demonstrated to be successful in generating germplasm for HS environments such as ME5 [33] . Seeds were sown on raised beds 2 m long with two rows (25 cm between rows) and 80 cm between the beds. Appropriate fertilizers, pesticides, and insecticides were applied according to normal practices prevalent in this region.

Phenotypic traits measurement
A high-resolution spectroradiometer (FieldSpec® Hand-Held 2, ASD Inc., Longmont, CO, USA) was used to measure the spectral reflectance of each plot. The HandHeld 2 is a hand-held spectro-radiometer that uses a FieldSpec 4 VNIR spectrometer for accurate analysis in the 325-1075 nm spectral range. To collect one spectral measurement, the sensor was located 50 cm over the canopy pointing in nadir position. A 12.7 cm Â 12.7 cm calibrated white panel (99% Spectralon®, Labsphere Inc., North Sutton, NH, USA) was used for the estimation of the incoming radiation and the spectral reflectance. An average spectrum per plot was estimated using four measurements collected at different locations across the plot and avoiding the first 50 cm at each extreme of the plot. The spectral measurements were conducted at 11:00-14:00 to avoid the possibility of confounding effects from variations in solar radiation. The spectral measurements were conducted at booting (Zadok 50), heading (Zadok 60), and 7 days after heading (heading plus 7 days, H7) stages in the YP environment, while in HS they were conducted at booting (Zadok 50) and heading (Zadok 60) stages. Several SRIs were derived from the collected data. The names, formulas and references can be found in Table 1.
In addition to SRIs, the following agronomic traits were measured: GY (g$m -2 ), thousand-grain weight (TGW, g), grain number [GN, estimated as GY (g$m -2 )/TGW (g) Â 1000], days to heading (DTH, i.e., the number of days from emergence when 50% of the spikes had emerged); days to maturity (DTM, i.e., the number of days from emergence to when the kernel became hard); plant height (PH, cm, from the soil surface to the tip of the spike without awns).

Statistical analysis
The best linear unbiased prediction (BLUP) values for agronomic traits and SRIs were calculated using META-R [50] . Correlations between the traits were estimated using SAS v9.1 (SAS Institute, Cary, CA, USA). Broad sense heritability (H 2 ) was estimated as:
where 2 g and 2 e were genotype and error variance, and r was the number of replications.
For the selection of efficient SRIs that could predict GY, two steps of analysis were conducted: (1) to identify the SRIs that were significantly correlated with GY (P < 0.01) at different crop stages; (2) to remove redundant SRIs based on correlations and only keeping those with relative highest H 2 in the same correlation cluster.

Marker trait associations (MTAs)
GWAS analysis was performed in TASSEL v5.0 software using the generalized linear model and mixed linear model [51] . Models were compared by fitting simple model, familial relatedness matrix (K) as random factor, population structure (Q 1 -Q 10 ) from STRUCTURE, principal components (PC 1 -PC 10 ), (Q 1 -Q 10 ) + K, and (PC 1 -PC 10 ) + K, and coefficient of parentage as a random factor in combination with the PC and Q matrix and observing the Q-Q plots of P-values. Manhattan plots were created using the 'qqman' package in R [52] . Once the MTAs were identified, comparisons were made between the traits and also with previously detected MTAs in the WAMI population [30][31][32] .

Phenotypic variation and broad-sense heritability
Agronomic traits showed a substantial reduction in trait values under HS compared to YP. GY, TGW, and GN decreased by 31%, 17%, and 17%, respectively. The length of the growth cycle also decreased under HS compared to YP, where the average DTH and DTM were shortened by 14 and 38 days, respectively. PH reduced by 16 cm on population average under HS (Table S1). Similar to the results found with agronomic traits, the values of the SRIs decreased in HS compared to YP, while the water indices (WI, NWI_1, NWI_2, NWI_3, and NWI_4), TCARI (transformed chlorophyll absorption in reflectance index), AB ratio (Chlorophyll ratio), ARI1 (anthocyanin reflectance index 1), ARI2, and NDPI (normalized difference pigment index) increased under HS, compared to YP ( Table 2).
Under YP, the SRIs related to light use efficiency (e.g., PRI, SIPI2 (structure insensitive pigment index 2)) showed an increasing trend from booting to late heading (H7), whereas SRIs related to vegetation (SR, NDVI_670, NDVI_705, and MTCI1 (MERIS terrestrial chlorophyll index 1)) and water status (WI, NWI_1, NWI_2, NWI_3, and NWI_4) showed a decreasing trend during the same period. Under HS, the SRIs related to chlorophyll (RARSa (ratio analysis of reflectance spectra for chlorophyll a), RARSb, ARI1, and ARI2) and carotenoid (NDPI and PSRI (plant senescence reflectance index)) content increased from booting to the heading stage, and SRIs related to vegetation (SR, NDVI_670, NDVI_705, EVI (enhanced vegetation index), MTCI1, OSAVI (optimized soil adjusted vegetation index), and TCARI) decreased during the same period, while indices related to water status (WI, NWI_1, and NWI_4) did not change much between the two stages. Most of the SRIs had similar trends under YP and HS, while TCARI, PRI, NWI_2, NWI_3, and RARSb showed contrasting patterns between the two environments.
The highest H 2 for most of the SRIs were observed at the heading stage under YP. Under HS, moderate to high H 2 was observed at both booting and heading stage ( Table 2).

Correlations
Under YP, half of the studied SRIs showed significant correlations with GY at booting and heading, whereas only four SRIs showed good correlations at late heading (P < 0.01, Table 2). EVI, OSAVI, and ARI1 were correlated with GY across the three measured stages. Under HS, eight and 19 SRIs were correlated significantly with GY at booting and heading stage, respectively. NDVI_705, MTCI1, WI, NWI_1, NWI_2, NWI_3, NWI_4, and PRI were correlated with GY at both stages.
Significant correlations between the SRIs were observed in both YP and HS environments (Table S2). Most of the high correlations were observed in the same SRIs group (e.g., NDVI_670, and NDVI_705; WI, NWI_1, NWI_2, NWI_3, and NWI_4), while some strong correlations were found between SRIs from different groups (e.g., SR with RARSc; SIPI2 with NDVI_670; and PRI with NDQI).

Efficient SRIs for indirect selection of GY
According to the heritability estimates and correlations, several efficient SRIs were identified as indirect indices for GY at different growth stages: under YP, EVI was an efficient selection index from booting until late heading stage; NDPI, NPQI (normalized phaeophytinization index), PRI, ARI2, MTCI4, and NWI_1 were efficient SRIs at both booting and heading stages; NDVI_705 was efficient at booting stage; CRI2 (carotenoid reflectance index 2) was efficient from heading to late heading; while ARI1 and OSAVI were more efficient after heading ( Fig. 1(a)). Under HS, NDVI_705 and PRI were efficient in predicting GY at both booting and heading stages; while NW1_2 was efficient only at the booting stage, and RARSb, ARI2, AB ratio, PSSRa (pigment specific simple ratio a), MTCI2, and WI were efficient at the heading stage ( Fig. 1(b)). NDVI_705 and PRI were common efficient SRIs for YP and HS, while RARSb, AB ratio, PSSRa, and WI were specific SRIs for HS.  (Table S3).

Comparison with the previously identified marker-trait associations on the WAMI population
The WAMI population has been widely used for association mapping studies. In the present study, the locus in chromosome 1B at 70 cM for GN was previously reported for peduncle length [31] and spike ethylene [53] (Table S4). The other locus in chromosome 1B at 148 cM for PH was earlier reported for harvest index and TGW [31,54] . The locus in 2A (143 cM) for EVI was co-located QTL for TGW [54] . The locus on 3B (61 cM) for PH was earlier identified for SPAD at the grain-filling stage [31] . On 4B, the locus at 66 cM for DTM has been reported for GY [54] , and the locus at 71 cM for EVI and OSAVI was reported for spike dry weight [53] . On 5A, the locus at 60 cM for PRI was previously reported for TGW [54] , while the locus at 90 cM was reported for DTH and flowering time [1,30,31] and was identified as the Vrn-1 locus. On 6A, the locus at 80 cM for TGW was reported to co-located with QTL for PH [31] , and the locus at 85 cM for GY was reported for GY, biomass, TGW, and SPAD at grain-filling [31,54] . The locus on 7A (35 cM) for NPQI was previously identified for GY [54] .

Discussion
Identification of indirect selection indices for GY in target environments is important in wheat breeding programs with large sets of germplasm and genotypes. An efficient indirect index should be significantly correlated with GY, showing high heritability and genetic variability [55] . In the case of physiological traits, the phenological stage at which the measurements are taken also may affect the ability to predict GY. In the current study, most SRIs were associated with GY at booting and heading stages under YP, but high heritability estimates of SRIs was observed at the heading stage. Under HS, significant correlations and high heritability of SRIs were observed at both booting and heading stages. These results suggest that the heading stage could be the most optimum stage to conduct the measurement of SRIs for early selection/prediction of GY under both conditions. Similar observations were also reported in irrigated spring wheat [23,56] and durum wheat [57,58] . In addition, we also identified several efficient SRIs at booting and late heading stages, which could be used for phenotyping large populations.  2 Manhattan plots for spectral reflectance indices (SRIs) of wheat association mapping initiative (WAMI) population at booting, heading, and heading plus 7 days under yield potential (YP) ((a)-(f)) and heat stress (HS) ((g)-(j)) environments. NDPI, normalized difference pigment index; EVI, enhanced vegetation index; CRI, carotenoid reflectance index; ARI, anthocyanin reflectance index; NDVI, normalized difference vegetation index; NWI, normalized water index; AB ratio, chlorophyll ratio.
In the present study, common efficient SRIs such as NDVI_705, PRI, and ARI2 were identified for both YP and HS conditions, highlighting their values as proxies for GY. The water indices (WI and NWI) also showed significant correlations with GY under YP and HS conditions. The water-based indices estimate water content of the canopy, and the higher associations between these indices and GY indicated that canopy water content plays a vital role in the productivity of a given genotype. Similar results were observed by Prasad et al. [56] under rain-fed conditions. Further, the simple ratio index (WI) and normalized indices (NWI_1, NWI_2, NWI_3, and NWI_4) of water status showed similar predictive ability for GY in this study, confirming the findings in previous studies of spring wheat [12,23,56] .
To the best of our knowledge, this is the first study to report genomic regions associated with SRIs under HS. Recently, Gizaw et al. [27,28] identified the genomic regions for SRIs through GWAS in fully irrigated, rainfed, and drought-stressed environments in wheat. They reported that under irrigation, QTLs for NDVI were located on chromosome 3A, 3B, 4A, 5A, 5B, 7A, and 7B in winter wheat, and on 2B, 3B, 4B, 4D, 5B, and 7A in spring wheat. In contrast, in this study, we identified QTLs for NDVI_705 on 1B, 3A, and 5A, which were far away from the QTLs reported by Gizaw et al. [27,28] . The QTLs for NWI_1 under irrigation were reported on 6A, 6B [28] , and 2A, 5A, 5B [27] in winter and spring wheat, respectively; while QTLs for NWI_1 were located at a different locus on 5B in this study. Previously QTLs for Fig. 3 Genomic regions associated with agronomic traits and spectral reflectance indices (SRIs) in the wheat association mapping initiative (WAMI) population at booting (purple), heading (green), and heading plus 7 days (blue) under yield potential (YP) and at booting (orange) and heading (red) stages under heat stress (HS) environments. NDVI, normalized difference vegetation index; ARI, anthocyanin reflectance index; GN, grain number; RARS, ratio analysis of reflectance spectra for chlorophyll; PH, plant height; MTCI, MERIS terrestrial chlorophyll index; DTH, days to heading; GY, grain yield; TGW, thousand-grain weight; DTM, days to maturity; NDPI, normalized difference pigment index; PRI, photochemical reflectance index; CRI, carotenoid reflectance index; PSSR, pigment specific simple ratio; OSAVI, optimized soil adjusted vegetation index; EVI, enhanced vegetation index; NWI, normalized water index; NPQI, normalized phaeophytinization index; AB ratio, chlorophyll ratio; WI, water index. PRI were reported on 2A, 3A, 4A, 6B, and 7A in winter wheat under irrigation [28] , while in this study, QTLs for PRI were identified on 3B, 4B, 5A, and 6A. Also, the marker wsnp_Ex_c8884_14841846 (3A, 105 cM) accociated with CRI2 under YP in this study was reported to associate with PRI under drought, rainfed, and irrigation conditions [28] . The markers wsnp_Ex_c22727_31934296 and wsnp_Ex_c7729_13177883 (5A, 89-91 cM) for DTH in this study were previously reported associated with SR and NDVI [28] . Among the SRIs, NDVI is the most used and consistent index for in-season selection and yield prediction [25,26] . The QTLs for NDVI_705 on 1B (60-62 cM) under YP and on 2B (27 cM) and 6B (71 cM) under HS were very close to the reported loci in durum wheat under drought stress [59] . The locus for NDVI_705 on 7B (152 cM) under HS was about 10 cM away from a NDVI QTL detected in a bi-parental population under heat [60] . Our study identified several genomic regions for SRIs under YP and HS, while few common loci for both environments were observed. These observations again confirmed that different stresses on crops at the same time, or at different times during the same growing season, may produce different metabolic and physiological responses [61] . The WAMI population was widely used for GWAS studies, and some loci identified in our analysis were located closely to previously reported QTLs for GY and yield components in the same population. For example, the locus on chromosome 1B for NDVI and ARI2 was 2 cM away from the adaptation to density locus (at 64 cM) [32] ; the locus on 2A for EVI was co-located with that for TGW [54] ; the locus on 4B for EVI and OSAVI was reported for spike dry weight [53] ; the locus on 5A for PRI was previously reported for TGW [54] ; the locus on chromosome 6A for PRI and NWI_2 was previously reported to be associated with GY, biomass, PH, and chlorophyll content under irrigated YP conditions [31] ; and the locus on 7A for NPQI was previously identified for GY [54] . Also, the locus for ARI2 and PRI was co-located with DTH and PH on the Vrn-1 locus on 5A, verifying the common sense assumption that major developmental genes modulate key physiological processes. BLAST results of the significant markers associated with multiple SRIs indicated that most of the candidate genes were related to transportation, metabolism, senescence, and abiotic and biotic resistance in plants (Table S3) associated with NDPI and ARI2, was a FCS-Like Zinc finger gene regulated by sugars, cellular energy level, and abiotic stress [62] . A gene coding ABC transporter B family protein (TraesCS3A01G285200.1) was also associated with NDPI and ARI2, and this protein is a powerful transporter driving the exchange of compounds across biological membranes [63] . A blue copper protein gene (TraesCS4B01G271800.1) associated with EVI and OSAVI, can shuttle electrons from a protein acting as an electron donor to another acting as an electron acceptor in photosynthesis [64] . A Myb transcription factor gene (TraesCS3A01G367600.1) associated with NDPI, ARI2, and CRI2 was involved in controlling various processes, such as responses to biotic and abiotic stresses, development, differentiation, metabolism, and defense in plants [65] . A candidate gene (TraesCS3A01G009200.1) associated with NDVI_705 and MCI2 was characterized as a disease resistance gene encoding the NBS-LRR-like resistance protein [66] . Although our measurements cover phenological stages that are highly relevant for adaptation to stress conditions, further measurements of SRIs after anthesis would be needed to fully understand their relationships with yield. Nevertheless, our results have great potential value since the SRIs reported here were analyzed in elite populations under HS and in terms of their genomic regions for the first time. Further research supporting these findings in similar environments using different germplasm will help to reinforce our results. Analysis of the same genomic regions under different severity of HS will provide useful data for breeders and physiologists to help underpin the genomic bases of the traits. Additionally, the exploration of these traits under drought conditions, or combined heat and drought conditions would enable the identification of commonalities between different stresses.

Conclusions
Spectral reflectance indices are valuable tools to screen agronomic and physiological traits in large populations. In this study, we identified EVI as the SRI that was correlated with GY under YP condition at the booting, heading, and heading plus 7 days stages. PRI and NDVI were efficient under HS at booting and heading stages. Several SRIs showed high heritability estimates and marker-trait associations were identified for them. This information can be further used in developing molecular markers for selecting high yielding lines under YP and HS conditions.