QTL Mapping of Leaf Area Index and Chlorophyll Content Based on UAV Remote Sensing in Wheat

Abstract

High-throughput phenotypic identification is a prerequisite for large-scale identification and gene mining of important traits. However, existing work has rarely leveraged high-throughput phenotypic identification into quantitative trait locus (QTL) acquisition in wheat crops. Clarifying the feasibility and effectiveness of high-throughput phenotypic data obtained from UAV multispectral images in gene mining of important traits is an urgent problem to be solved in wheat. In this paper, 309 lines of the spring wheat Worrakatta × Berkut recombinant inbred line (RIL) were taken as materials. First, we obtained the leaf area index (LAI) including flowering, filling, and mature stages, as well as the flag leaf chlorophyll content (CC) including heading, flowering, and filling stages, from multispectral images under normal irrigation and drought stress, respectively. Then, on the basis of the normalized difference vegetation index (NDVI) and green normalized difference vegetation index (GNDVI), which were determined by multispectral imagery, the LAI and CC were comprehensively estimated through the classification and regression tree (CART) and cross-validation algorithms. Finally, we identified the QTLs by analyzing the predicted and measured values. The results show that the predicted values of determination coefficient (R²) ranged from 0.79 to 0.93, the root-mean-square error (RMSE) ranged from 0.30 to 1.05, and the relative error (RE) ranged from 0.01 to 0.18. Furthermore, the correlation coefficients of predicted and measured values ranged from 0.93 to 0.94 for CC and from 0.80 to 0.92 for LAI at different wheat growth stages under normal irrigation and drought stress. Additionally, a linkage map of this RIL population was constructed by 11,375 SNPs; eight QTLs were detected for LAI on wheat chromosomes 1BL, 2BL (four QTLs), 3BL, 5BS, and 5DL, and three QTLs were detected for CC on chromosomes 1DS (two QTLs) and 3AL. The closely linked QTLs formed two regions on chromosome 2BL (from 54 to 56 cM and from 96 to 101 cM, respectively) and one region on 1DS (from 26 to 27 cM). Each QTL explained phenotypic variation for LAI from 2.5% to 13.8% and for CC from 2.5% to 5.8%. For LAI, two QTLs were identified at the flowering stage, two QTLs were identified at the filling stage, and three QTLs were identified at the maturity stage, among which QLAI.xjau-5DL-pre was detected at both filling and maturity stages. For CC, two QTLs were detected at the heading stage and one QTL was identified at the flowering stage, among which QCC.xjau-1DS was detected at both stages. Three QTLs (QLAI.xjau-2BL-pre.2, QLAI.xjau-2BL.2, and QLAI.xjau-3BL-pre) for LAI were identified under drought stress conditions. Five QTLs for LAI and two QTLs for CC were detected by imagery-predicted values, while four QTLs for LAI and two QTLs for CC were identified by manual measurement values. Lastly, investigations of these QTLs on the wheat reference genome identified 10 candidate genes associated with LAI and three genes associated with CC, belonging to F-box family proteins, peroxidase, GATA transcription factor, C₂H₂ zinc finger structural protein, etc., which are involved in the regulation of crop growth and development, signal transduction, and response to drought stress. These findings reveal that UAV sensing technology has relatively high reliability for phenotyping wheat LAI and CC, which can play an important role in crop genetic improvement.

1. Introduction

Wheat is one of the most important staple crops for 4.4 billion people in the world [1,2]. Sustainable wheat production is an important guarantee for food security, which depends on the ability to adapt to environmental stress and extreme climate changes [3,4]. It was reported that drought and high temperatures have reduced the world’s annual wheat output by 9% to 10% [4]. As important phenotypic traits of crops, LAI and CC can effectively reflect crop growth, plant health, etc. and are also an important basis for predicting yield [5]. Therefore, real-time, rapid monitoring and accurate estimation of wheat LAI and CC are helpful to understand the growth information of wheat in a timely manner and guide field management [6].

As more and more crop genomes have been deciphered, the lack of phenotype information corresponding to the genome has become the main bottleneck restricting the exploration of the genotype–phenotype relationship [7]. The lack of efficient methods for phenotype collection and analysis has weakened the field of crop breeding [8]. Low throughput and efficiency, as well as heavy labor and time consumption, are the main problems when collecting phenotypic information in a large area using traditional methods. Using remote sensing technology to estimate vegetation has been widely developed especially for large-scale and long-term crop monitoring. Compared with traditional methods, UAV remote sensing technology has more advantages, e.g., convenient transportation, high flexibility, short operation time, and low investment, which provide a new method for fast, nondestructive, and high-throughput acquisition of field crop phenotypic information [9,10,11,12]. At the same time, many studies found that crop canopy spectral information was significantly related to crop LAI and CC, and they used spectral reflectance and the vegetation index to build mathematic models for estimating LAI and CC [8,13,14,15]. Gong et al. [14] developed a simple method to remotely estimate LAI with UAV imaging for a variety of rice cultivars throughout the entire growing season. Hunt et al. [16] used the UAV multispectral imagery-based green normalized difference vegetation index (GNDVI) to estimate wheat LAI by remote sensing and obtained good prediction (R² = 0.85). Singhal et al. [17] evaluated an advanced machine learning algorithm kernel ridge regression combined with spectral information and ground-truth chlorophyll data to model the chlorophyll estimation, resulting in high correlation (R² = 0.7452). Jin et al. [13] used hyperspectral reflectance and the vegetation index to establish wheat CC estimation models (R² > 0.94) under two growth environments, irrigated land and dry land. These studies indicate that remote sensing technology has become an alternative for traditional phenotype collection and an important tool for crop growth monitoring [10,18].

With the application of UAV remote sensing technology in agriculture, phenotypic data have been recently collected from wheat for use in genetic studies. Combining high-throughput phenotypic data with QTL mapping and a genome-wide association study (GWAS), several agronomic traits in wheat were efficiently investigated. Condorelli et al. [19] used UAV multispectral remote sensing and manual measurement to obtain the normalized difference vegetation index (NDVI) for 248 durum wheat lines, performed GWAS analysis on the population, and found that the number of QTLs detected with UAV-based NDVI increased by 25% compared to the QTLs detected with manual measurement, with 46 QTLs overlapping for the two phenotypic data sources. The individual loci explained 2.69–8.43% of the phenotypic variance. Adeel et al. [20] used UAV remote sensing technology and artificial measurement to collect the plant height of a wheat double haploid (DH) population consisting of 198 wheat lines, and they found that the values of plant height based on imagery were significantly related to the measured values, with five height-related QTLs mapped by the UAV-predicted value, i.e., two more QTLs than detected by the measured value. Among them, three QTLs were coincidently detected by both manual measurement and UAV sensing-based data.

The LAI and CC of wheat are quantitative traits regulated by multiple genes and determined by the interaction between genes and the environment [21]. The genetic characteristics of LAI and CC were investigated in previous reports [5,13,22], but QTL mapping for LAI and CC on the basis of UAV remote sensing technology in wheat is rare. In our previous study, we constructed linkage maps with a 50K SNP array for an RIL population of 309 lines from wheat cross Worrakatta × Berkut, and we conducted QTL mapping of drought resistance at the germination stage [23].

Thus, the major objectives of the present study were to (1) standardize a rapid method for LAI and CC estimation using a UAV platform, (2) identify quantitative trait loci for LAI and CC using UAV and ground-based measurements, and (3) assess genomic prediction accuracy for LAI and CC in wheat.

2. Materials and Methods

2.1. Experimental Materials

The RIL population consisting of 309 lines from the cross of spring wheat cultivars Worrakatta × Berkut was originally from the International Maize and Wheat Improvement Center (CIMMYT). The RIL and parents were planted at the Manas experimental station (86°12′52.2″ E, 44°18′15.7″ N) of the Xinjiang Academy of Agricultural Sciences in 2020. The Manas experimental station has a mid-temperate continental semiarid climate with severe cold winters, hot summers, dryness and low rainfall, sufficient sunshine, high evaporation, and low precipitation. The experiment consisted of two water supply treatments, normal irrigation (NI) and drought stress (DS), and each treatment was repeated twice. Approximately 100 seeds of each RIL and parents were sown in a single row with a row length of 2 m and a space between two rows of 20 cm. Field planting adopts a random block design. Drip irrigation was applied to both treatments; however, watering was stopped after the wheat heading stage in the drought stress treatment. Other field management approaches, such as the fertilization, insect control, and weeding followed local regulations.

2.2. Ground Data Collection

2.2.1. Determination of LAI

The LAI value was determined using an LAI-2200C plant canopy analyzer (LI-COR, Lincoln, NE, USA) at the flowering, filling, and maturation stages of wheat. The value of LAI was calculated using one sky value and four target values for each family and the parents. The sky value was used as the calibration value, and the average of the four target values was taken as the LAI value for each family and the parents. When measuring, the probe of the analyzer was located at the root of the wheat plant and the lens was kept level, avoiding direct sunlight as much as possible.

2.2.2. Determination of CC

Soil plant analysis development (SPAD) can directly reflect the relative chlorophyll content in leaves [24]. The CC was measured with a chlorophyll meter (SPAD-502 Plus, Minolta, Japan). For each RIL family and parents, five wheat plants at the same growth stage were selected and marked for measurement. The CC value on the flag leaf in the main stem was measured at the wheat growth stages of heading, flowering, and grain filling. The average value of five plants in each RIL family or parent was used as the CC value for the subsequent analysis.

2.3. UAV Image Acquisition and Preprocessing

Wheat plant images were collected using an agriculture drone DJI Phantom 4 (Shenzhen Dajiang Technology Co. Ltd., Shenzhen, China) mounting multispectral camera sensors with customer-specified bandpass filters at wavelengths of 840 nm ± 26 nm, 717 nm ± 16 nm, 668 nm ± 16 nm, 560 nm ± 16 nm, and 475 nm ± 16 nm, corresponding to single-color near-infrared (NIR), red-edge, red, green, and blue bands, respectively (Figure 1). Wheat images were collected during a sunny, windless, and cloudless day on 29 May 2020 (heading), 2 June (flowering), 10 June (filling), and 24 June (maturity). During the data acquisition process, the drone flew autonomously according to the set route and recorded images at the same time. The parameters for flight were set as an altitude of 12 m, speed of 1 m/s, head overlapping of 75%, and side overlapping of 75%. The acquisition time was synchronized with the ground data collection.

Pix4Dmapper software (Version 1.4, Pix4d, Lausanne, Switzerland) (https://pix4d.com/) (accessed on 1 June 2021) was used to stitch and process the acquired UAV multispectral images. Briefly, we used the corresponding ground control point data to correct the multispectral images to generate a digital orthophoto map (DOM). Then, the gray board was used to perform reflectance correction on the multispectral images to obtain the reflectance images of the test site, which were saved in TIF format. Finally, the single-band spectral reflectance images of the four growth periods were cropped by ENVI software (Version 5.1, ENVI, Fort Collins, CO, USA, Exelis Visual Information Solutions) (http://www.exelisvis.com) (accessed on 17 June 2021), and the single-band average reflectance of each center of the wheat plot in the study area was extracted as the spectral reflectance of the plot for the further analysis.

2.4. Selection of Vegetation Index

The combination of changes in the reflectance of different bands constitutes vegetation indices, which can reduce the degree of influence of factors such as background soil on vegetation spectra to a certain extent and improve the accuracy of estimating LAI and CC. In this study, five kinds of spectral reflectance (RNIR, Rred-edge, Rgreen, Rred, and Rblue) were used. Two vegetation indices, normalized difference vegetation index (NDVI) [25] and green normalized difference vegetation index (GNDVI) [15], were used for LAI and CC estimation. NDVI can detect crop vegetation growth status and vegetation cover, while GNDVI can evaluate the photosynthetic activity of crops.

NDVI = (RNIR − Rred)/(RNIR + Rred).

(1)

GNDVI = (RNIR − Rgreen/(RNIR + Rgreen).

(2)

2.5. Algorithm Development for LAI and CC Estimation

The minimize squared error as a criterion was used for the classification and regression tree (CART) [26] to predict LAI and CC. The algorithmic model was built as follows: X and Y were assumed to be the input and output variables, respectively, for which, Y was a continuous variable for measurement of LAI or CC, and X was {GNDVI, NDVI}. Given a training dataset of $D=\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}$, where $x_i=({x_i}^{(1)},{x_i}^{(2)},\dots ,{x_i}^{(n)})$ is the eigenvector with n number of features, $i=1,2,\dots ,N$, and N is the total sample size. The partitioning of the feature space was done heuristically, i.e., each partition examined all the values of all the features in the current set one by one and selected the optimal one as the cutoff point according to the squared error minimization criterion. For example, for the j-th feature variable in the training set with value S as the cutoff variable and cutoff point, $R_1(j,s)=\left\{x\right|x^{(j)}\le \mathrm{s}\}$ and $R_2(j,s)=\left\{x\right|x^{(j)}>\mathrm{s}\}$ were used to find the optimal j and s. Thus, j and s were obtained through Equation (3) by minimizing the sum of squared errors of the two regions divided, where $c1=average(y_i|x_i\in R_1(j,s)),c2=average(y_i|x_i\in R_2(j,s))$. After finding the optimal cutoff point (j, s), the input space was divided into two regions in turn, and then the above division process was repeated for each region until the stopping condition was satisfied. Finally, a regression model was generated.

$$\underset{j,s}{\min }\left[\underset{c1}{\min }{\displaystyle \sum _{x_i\in R_1(j,s)}{(y_i-c1)}^2+\underset{c2}{\min }{\displaystyle \sum _{x_i\in R_2(j,s)}{(y_i-c2)}^2}}\right].$$

(3)

In the modeling process, 70% of the samples were selected as the modeling set, and 30% of the samples were selected as the validation set. The Kennard–Stone (KS) [27] method was used in the modeling set. The accuracy of the cross-validation procedure was evaluated by the coefficient of determination (R²), root-mean-square error (RMSE), and relative error (RE). The formulas of these statistics are shown below.

$${{\displaystyle R}}^2=\frac{{\displaystyle \sum _{i=1}^n{\left(\widehat{{{\displaystyle y}}_i}-\stackrel{-}{{\displaystyle y}}\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left({{\displaystyle y}}_i-\stackrel{-}{{\displaystyle y}}\right)}^2}},$$

(4)

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(\hat{{\displaystyle {{\displaystyle y}}_i}}-{{\displaystyle y}}_i\right)}^2}}{n}},$$

(5)

$$\mathrm{RE}=\frac{{\displaystyle \sum _{i=1}^n\left|\frac{{{\displaystyle y}}_i-\hat{{\displaystyle {{\displaystyle y}}_i}}}{{{\displaystyle y}}_i}\right|}}{n},$$

(6)

where $\widehat{y_i}$ is the predicted value of LAI or CC, $y_i$ is the measured value of LAI or CC, $\stackrel{-}{y}$ is the average value of LAI or CC, and n is the number of verification samples.

2.6. Data Analysis

For the experimental area, we had two replicates, one repetition for model building and evaluation and the other for application of the model. The statistics of all data were implemented in Spyder using Python 3.8.8 on a workstation with an Intel i7-6800 K 3.40 GHz CPU, 16 GB memory, and Nvidia GeForce GTX 2080Ti graphics, running the Win10 operating system. Software Pandas 1.3.2 (Greenwich, CT, USA), Matplotlib 3.4.2 (Washington, DC, USA), Scikit-Learn 0.24.2 (Paris, France), and SPSS 21.0 (Chicago, IL, USA) were applied for statistical analysis, correlation analysis, and the significance of differences test for wheat LAI and CC. In addition, the generalized heritability was calculated using the following equation:

$$h^2=\frac{{{\delta }_G}^2}{{{\delta }_G}^2+\frac{{{\delta }_{GE}}^2}{r}+\frac{{{\delta }_E}^2}{nr}},$$

(7)

where ${{\delta }_G}^2$ is the genetic variance, ${{\delta }_{GE}}^2$ is the genotype × treatment interaction variance, ${{\delta }_E}^2$ is the error variance, n is the number of treatments, and r is the number of replications.

2.7. QTL Mapping of LAI and CC and Candidate Gene Identification

The linkage maps of this population were constructed using a wheat 50K SNP array in the previous study [23], which consisted of 28 linkage groups with 11,375 SNPs and 2220.26 cM in length, with an average of 1.38 cM per bin marker pair. For each trait of LAI and CC, the mean values of each F6 RIL from two replicates were used for descriptive statistics. QTL analysis was performed using inclusive composite interval mapping (ICIM) in the software QTL IciMapping V4.1 (http://www.isbreeding.net/) (accessed on 5 June 2021). The walking speed for genome scanning was 1.00 cM, with p = 0.001 in stepwise regression. Significant LOD thresholds were evaluated using 1000 permutations tests with type I error set at 0.05. QTLs were named following the rules of International Rules of Genetic Nomenclature (https://wheat.pw.usda.gov/ggpages/wgc/98/Intro.htm) (accessed on 12 June 2021); ‘xjau’ represents Xinjiang Agricultural University.

To predict candidate genes involved in QTL for LAI- and CC-related traits from this study, we searched the QTL flanking markers in the wheat reference genome. First, the sequences of the closest markers of the QTL were blasted against the genome sequences of wheat Chinese Spring in CSS database (IWGSC RefSeq v2.0, https://urgi.versailles.inra.fr/blast_iwgsc/blast.php) (accessed on 7 June 2021), and their physical intervals on the reference genome sequences were identified. Then, for each QTL, the annotations of genes delimited by SNP marker intervals were obtained in EnsemblePlants (http://plants.ensembl.org/index.html) (accessed on 11 June 2021), and the wheat transcriptional expression database was investigated in the Triticeae Multi-omics Center (http://202.194.139.32/expression/index.html) (accessed on 14 June 2021).

3. Results

3.1. Multispectral Analysis of Reflectance

The reflectance of five lights (NIR, red-edge, green, red, blue) on wheat plants was compared at the heading, flowering, filling, and maturity stages under normal irrigation (NI) and drought stress (DS) conditions (

Table 1. Mean and standard deviation (SD) of NIR, red-edge, green, red, and blue reflectance for wheat at different growth stages under normal irrigation (NI) and drought stress (DS) conditions.

Water Condition	Bands	Parameters	Heading	Flowering	Filling	Maturity
DI	NIR	Mean	0.37	0.37	0.37	0.26
		SD	0.03	0.03	0.03	0.03
	Red-edge	Mean	0.28	0.28	0.28	0.26
		SD	0.03	0.03	0.03	0.02
	Green	Mean	0.07	0.07	0.07	0.09
		SD	0.01	0.01	0.01	0.01
	Red	Mean	0.05	0.05	0.05	0.10
		SD	0.01	0.01	0.01	0.02
	Blue	Mean	0.04	0.04	0.04	0.04
		SD	0.01	0.01	0.01	0.01
DS	NIR	Mean	0.33	0.32	0.30	0.28
		SD	0.03	0.03	0.02	0.02
	Red-edge	Mean	0.28	0.29	0.27	0.30
		SD	0.02	0.02	0.02	0.02
	Green	Mean	0.08	0.09	0.10	0.11
		SD	0.01	0.02	0.02	0.01
	Red	Mean	0.07	0.08	0.11	0.14
		SD	0.02	0.02	0.04	0.03
	Blue	Mean	0.04	0.06	0.05	0.05
		SD	0.01	0.03	0.01	0.01

).

Different growth stages for LAI and CC were studied. We focused on the flowering, filling, and maturation stages of wheat for LAI and heading, flowering, and grain filling stages for CC. The reflectance changes of the five bands corresponding to wheat LAI were obviously different. Under two water supply treatments, the order of single-spectral reflectance of each growth period was as follows: NIR > Red-edge > green > red > blue. For the normal irrigation treatment, the spectral reflectance of single light was 0.26–0.37 (SD 0.03) in the NIR band, 0.25–0.28 (SD 0.02–0.03) in the red-edge light band, 0.07–0.09 (SE 0.01) in the green light band, 0.05–0.10 (SD 0.01–0.02) in the red light band, and 0.04–0.04 (SD 0.01) in the blue light band. For the drought stress treatment, the spectral reflectance of single light was 0.28–0.32 (SD 0.02–0.03) in the NIR band, 0.29–0.30 (SD 0.02) in the red-edge light band, 0.09–0.10 (SD 0.01–0.02) in the green light band, 0.08–0.14 (SD 0.02–0.04) in the red light band, and 0.05–0.06 (SD 0.01–0.03) in the blue light band. Comparing the two water supply treatments, the values of LAI under drought stress were extremely lower than those under normal irrigation, and the five spectral reflectance values under drought stress were generally greater than those under normal irrigation. The reason could be drought stress reduced effective leaf area for light exposure.

The single-spectral reflectance of wheat CC under normal irrigation was 0.36–0.37 (SD 0.03) for NIR light, 0.25–0.29 (SD 0.03) for red-edge side light, 0.07 (SD 0.01) for green light, 0.05–0.06 (SD 0.01) for red light, and 0.04 (SD 0.01) for blue light. The single-spectral reflectance of wheat CC under drought stress was 0.30–0.33 (SD 0.02–0.03) for NIR light, 0.27–0.28 (SD 0.02) for red-edge light, 0.08–0.10 (SD 0.01–0.02) for green light, 0.07–0.11 (SD 0.02–0.04) for red light, and 0.04–0.06 (SD 0.01–0.03) for blue light. Under the two water supply treatments, the values of spectral reflectance were gradually decreased in the following order: NIR > red-edge > green > red > blue; however, the reflectance of the blue, red, and green light bands in the normal irrigation treatment were higher than that under drought stress. With the increase in CC under normal irrigation, the reflectance of the green light, red light, and blue light bands showed a slight upward trend, whereas the red-edge and NIR bands increased first and then decreased. Under drought stress, the red-edge band showed an upward trend, while the NIR band showed an increase first and then decreased. Under the two water treatments, the reflectance of the NIR band and the blue band was the highest and the lowest, respectively.

3.2. Estimation of LAI and CC by CART Modeling Verification Analysis

LAI and CC were estimated using the CART model, and its accuracy was evaluated by the determination coefficient (R²), root-mean-square error (RMSE), and relative error (RE) (

Table 2. Determination coefficient (R²), root-mean-square error (RMSE), and relative error (RE) of the algorithm for modeling the estimation of leaf area index (LAI) and chlorophyll content (CC) of wheat at different growth stages under normal irrigation (NI) and drought stress (DS) conditions.

Traits	Growth Stage	Treatment	R2	RMSE	RE
LAI	Flowering stage (LAI1)	NI	0.70	0.36	0.08
		DS	0.78	0.32	0.10
	Filling stage (LAI2)	NI	0.81	0.30	0.06
		DS	0.73	0.39	0.12
	Mature stage (LAI3)	NI	0.84	0.44	0.11
		DS	0.64	0.52	0.18
CC	Heading stage (CC1)	NI	0.88	0.99	0.01
		DS	0.86	0.97	0.01
	Flowering stage (CC2)	NI	0.88	1.02	0.01
		DS	0.88	0.94	0.01
	Filling stage (CC3)	NI	0.88	1.00	0.01
		DS	0.87	1.05	0.02

, Figure 2 and Figure 3). GNDVI and NDVI with high correlation coefficients were employed for LAI and CC estimation modeling, and we randomly selected 70% of the data for model building and 30% of the data for model validation. For LAI, the estimated values and measured values were positively correlated with a high coefficient (Table 2, Figure 2). Under normal irrigation conditions, the correlation was increased from the flowering to maturity stage, and a higher correlation was observed at the maturity stage with measured LAI > 2.0 (Figure 2). Under drought stress conditions, a similar positive correlation was found between predicted and measured values of LAI at the three growth stages, although the coefficient was slightly decreased after the flowering stage (Figure 2).

The relationship between predicted and measured values of CC also had a positive correlation (Table 2, Figure 3). High coefficients (between 0.86 and 0.88) and low relative errors (between 0.01 and 0.02) were found between the predicted and measured values at the three growth stages under both normal irrigation and drought stress conditions (Table 2).

Overall, the high R² and low RE values from the above results indicated the feasibility of predicting LAI and CC under different water treatments and different growth periods using CART algorithm modeling.

3.3. Genetic Variation of LAI and CC in Wheat RIL Population

The RIL populations showed some differences in LAI and CC at different growth stages under normal irrigation and drought stress (

Table 3. Genetic variation of measured LAI and CC of wheat cultivars Worrakatta and Berkut and their progeny recombinant inbreeding lines (RILs).

Traits	Growth Stage	Treatment	Parent		RILs
			Worrakatta	Berkut	Mean	SD	Range	CV (%)	h2 (%)
LAI	Flowering stage	NI	4.57	4.39	4.21	0.63	2.66~5.86	15.4	70
		DS	4.6	3.58	3.07 **	0.64	1.62~5.08	21	70
	Filling stage	NI	5.32	4.57	4.39	0.66	2.79~6.29	14.9	73
		DS	4.42	3.8	2.94 **	0.65	1.31~4.67	22.2	69
	Mature stage	NI	4.57	3.95	3.85	0.77	2.24~6.05	19.9	73
		DS	4.26	3.65	2.91 **	0.82	1.01~5.00	28.2	75
CC	Heading stage	NI	53.38	46.11	51.63	2.32	44.06~57.68	4.5	81
		DS	53.41	48.9	50.76 **	1.99	42.41~55.08	3.9	75
	Flowering stage	NI	57.68	48.28	53.1	2.33	47.43~60.86	4.4	78
		DS	55.54	47.68	51.19 **	2.16	44.52~57.08	4.2	76
	Filling stage	NI	54.82	46.48	51.31	2.11	45.62~56.99	4.1	74
		DS	54.08	47.06	50.34 **	2.07	39.57~55.67	4.1	74

** p < 0.001; CV, coefficient of variation; h², broad-sense heritability; SD, standard deviation.

). Compared with the normal irrigation treatment, the LAI and CC of parent Worrakatta showed a decreasing trend under drought stress at different fertility stages, while parent Berkut showed a relatively stable performance. This shows that parent Berkut was better able to maintain a more stable LAI, as well as a good green holding capacity, under drought stress. Compared with the normal irrigation treatment, the LAI and CC of the RIL population showed a highly significant decrease in different fertility periods under drought stress, indicating that both LAI and CC were sensitive to drought stress. Under both treatments, LAI and CC of the RIL population showed a continuous distribution and wide range at different fertility periods and in the presence of super parental segregation, indicating that this population was suitable for QTL localization analysis. The broad-sense heritabilities of LAI and CC were both high, ranging from 69% to 75% and from 74% to 81%, respectively. This indicated that the phenotypic differences of these two traits were strongly influenced by genetic factors.

Under normal irrigation, the mean LAI of the RIL population ranged from 3.85 to 4.39, with the variation ranging from 2.24 to 6.29, and the CV ranging from 14.9% to 19.9%; the value of mean LAI reached its maximum (4.39) at the filling stage. The mean value ranged from 2.91 to 3.07, with the variation ranging from 1.01 to 5.08, and the CV ranging from 21.0% to 28.2%; the mean value LAI reached its maximum (3.07) at the flowering stage under drought stress. This result showed that drought stress could promote an increase in LAI in a short period of time; however, if the wheat in the field was under long-term drought stress, wheat growth and development were hindered and LAI decreased, which reduced wheat yield.

Under normal irrigation, the mean CC values of the RIL population ranged from 51.31 to 53.10, with the variation ranging from 44.06 to 60.86, and the CV ranging from 4.1% to 4.5%. The mean values ranged from 50.34 to 51.19, with the variation ranging from 39.57 to 57.08, and the CV ranging from 3.9% to 4.2% under drought stress. Under both treatments, the mean CC value reached its maximum at the flowering stage.

3.4. Genetic Variation of Predicted LAI and CC in Wheat RIL Population

The variations of predicted LAI and CC among parents and RIL populations were observed in different growth periods and under different water treatments (

Table 4. Genetic variation of predicted LAI and CC of wheat cultivars Worrakatta and Berkut and their progeny recombinant inbreeding lines (RILs).

Traits	Growth Stage	Treatment	Parent		RIL
			Worrakatta	Berkut	Mean	SD	Range	CV (%)	h2 (%)
LAI	Flowering stage	NI	5.24	4.41	4.1	0.47	1.63–5.70	11.5	57
		DS	3.02	3.53	3.06 **	0.44	1.85–4.87	14.5	58
	Filling stage	NI	5.25	5.05	4.4	0.48	3.31–5.87	10.8	65
		DS	3.34	4.52	3.06 **	0.5	1.71–4.84	16.4	58
	Mature stage	NI	4.39	4.36	4.01	0.63	0.60–5.85	15.6	93
		DS	3.76	4.27	2.92 **	0.62	0.95–4.51	21.4	62
CC	Heading stage	NI	53.59	47.2	51.62	1.54	45.31–56.33	3	80
		DS	52.69	48.96	50.50 **	1.36	44.13–53.92	2.7	64
	Flowering stage	NI	57.06	54.58	53	1.52	48.61–57.41	2.9	77
		DS	55.29	50.32	50.68 **	1.39	43.76–55.06	2.8	64
	Filling stage	NI	53.5	49.66	51.44	1.41	47.67–56.09	2.7	70
		DS	50.6	48.46	50.20 **	1.43	41.38–53.70	2.8	67

** p < 0.001; CV, coefficient of variation; h², broad-sense heritability; SD, standard deviation.

, Figure 4). The values of LAI and CC for parental cultivar Berkut were significantly higher (p < 0.001) than those for Worrakatta at all growth stages under drought stress treatments (Table 4). The estimated values of LAI and CC in the RIL population showed a normal distribution at all growth stages and under different water treatments (Figure 4), indicating that LAI and CC are quantitative traits controlled by multiple genes. The frequency distributions of LAI at the flowering and maturity stages were consistent under normal irrigation or drought stress; however, at the filling stage, the mean values under drought stress were lower than those under normal irrigation conditions (Figure 4, LAI2). For CC, the frequency distributions under drought stress conditions tended to shift to the left, with lower values than those under normal irrigation conditions at all three growth stages (Figure 4). The values of LAI and CC in the RIL population decreased at all growth stages under drought stress compared with normal irrigation conditions, and the heritability of LAI ranged from 57% to 93%, while that of CC ranged from 64% to 80% (Table 4).

A high positive correlation was detected between predicted and measured values for both LAI and CC at different wheat growth stages under two different treatments (

Table 5. Correlation analysis between measured and predicted values for leaf area index (LAI) and chlorophyll content (CC) of 309 recombinant inbreeding lines (RILs) from wheat cross Worrakatta × Berkut.

Traits	Growth Stage	Treatment	Predicted Values		Measured Values		Correlation
			Range	Mean	Range	Mean
LAI	Flowering stage	NI	1.65–5.70	4.1	2.66–5.86	4.21	0.83 **
		DS	1.80–5.00	3.06	1.62–5.08	3.07	0.88 **
	Filling stage	NI	2.63–5.87	4.4	2.79–6.29	4.39	0.90 **
		DS	1.69–4.84	3.06	1.31–4.67	2.94	0.85 **
	Mature stage	NI	0.60–5.85	4.01	2.24–6.05	3.85	0.92 **
		DS	0.95–5.13	2.92	1.01–5.00	2.91	0.80 **
CC	Heading stage	NI	45.31–56.33	51.62	44.06–57.68	51.63	0.94 **
		DS	44.13–53.92	50.5	42.41–55.08	50.76	0.93 **
	Flowering stage	NI	48.61–57.41	53	47.43–60.86	53.1	0.94 **
		DS	43.76–55.06	50.68	44.52–57.08	51.19	0.94 **
	Filling stage	NI	47.67–56.09	51.44	45.62–56.99	51.31	0.94 **
		DS	41.38–53.70	50.2	39.57–55.67	50.34	0.93 **

** Highly significant correlation at the 0.01 level.

). The correlation coefficient of CC was consistently 0.94 at three wheat growth stages under the normal irrigation treatment, while it was between 0.93 and 0.94 under the drought stress treatment (Table 5). The correlation coefficient of LAI was between 0.83 and 0.92 from the flowering to mature stage under normal irrigation, while it was between 0.80 and 0.88 under drought stress (Table 5).

3.5. QTL Mapping of LAI and CC

The genetic map of the RIL population consisted of 28 linkage groups including 11,375 SNP markers [23]. Using the complete interval QTL mapping method, eight QTLs for LAI and three QTLs for CC were mapped on nine chromosomal regions, with three QTLs on chromosome 2BL and one each on 1BL, 1DS, 3AL, 3BL, 5BS, and 5DL. Each QTL explained the phenotypic variation for LAI from 2.5% to 13.8% and for CC from 2.5% to 5.8%. Eight QTLs had positive additive effect values, indicating that the synergistic alleles at these loci derived from the parent Berkut, while two QTLs had negative values, indicating that the synergistic alleles at these loci derived from the parent Worrakatta (

Table 6. Information on QTLs related to leaf area index (LAI) and chlorophyll content (CC) based on RIL population of cross Worrakatta × Berkut.

Traits	Source	Growth Stage	Treatment	QTL Name	Chr.	Marker Interval	Position (cM)	LOD	R2 (%)	Additive
LAI	Predicted value	Flowering	NI	QLAI.xjau-2BL-pre.1	2BL	AX-89402509–AX-95134753	101	5.6	3.1	0.13
		Filling	NI	QLAI.xjau-5DL-pre	5DL	AX-94670615–AX-111652649	108	5.7	3.6	0.13
		Filling	DS	QLAI.xjau-3BL-pre	3BL	AX-179558689–AX-94903264	56	5.3	3.3	0.16
		Maturity	NI	QLAI.xjau-5DL-pre	5DL	AX-94670615–AX-111652649	108	4.1	2.5	0.13
		Maturity	DS	QLAI.xjau-2BL-pre.2	2BL	AX-95631217–AX-94667571	54	5.3	3.3	0.16
	Measured value	Flowering	NI	QLAI.xjau-2BL.1	2BL	AX-111757968–AX-110468087	96	3.4	6.8	0.17
		Flowering	NI	QLAI.xjau-5BS	5BS	AX-179559227–AX-179561411	127	4.2	8.2	0.18
		Maturity	NI	QLAI.xjau-1BL	1BL	AX-111775778–AX-109367690	38	2.8	7.6	−0.19
		Maturity	DS	QLAI.xjau-2BL.2	2BL	AX-94407412–AX-111743644	56	4.6	13.8	0.3
CC	Predicted value	Heading	NI	QCC.xjau-1DS-pre	1DS	AX-109983565–AX-89314186	27	4.1	2.5	0.31
		Heading	NI	QCC.xjau-3AL-pre	3AL	AX-110449574–AX-109301028	103	3.9	2.5	-0.3
	Measured value	Heading	NI	QCC.xjau-1DS	1DS	AX-110335177–AX-109983565	26	2.5	5.8	0.52
		Flowering	NI	QCC.xjau-1DS	1DS	AX-110335177–AX-109983565	26	3.2	5.3	0.57

, Figure 5).

Among the eight QTLs for LAI, four QTLs, QLAI.xjau-2BL-pre.1, QLAI.xjau-2BL-pre.2, QLAI.xjau-3BL-pre, and QLAI.xjau-5DL-pre, were detected on the basis of predicted LAI values and four QTLs, QLAI.xjau-1BL, QLAI.xjau-2BL.1, QLAI.xjau-2BL.2, and QLAI.xjau-5BS, were detected on the basis of ground-measured values. QLAI.xjau-5DL-pre was detected at both filling and maturity stages, accounting for 3.6% and 2.5% phenotypic variation (Table 6). Three QTLs for LAI, QLAI.xjau-2BL-pre.2, QLAI.xjau-2BL.2, and QLAI.xjau-3BL-pre, were identified under drought stress treatment, and five QTLs were identified under normal irrigation treatment. QLAI.xjau-2BL-pre.2 and QLAI.xjau-2BL.2 were very closely linked, 2 cM apart, and QLAI.xjau-2BL.2 explained 13.8% of the phenotypic variation at the maturity stage (Table 6).

For chlorophyll content, QCC.xjau-1DS was detected at both heading and flowering stages and explained similar phenotypic variation (5.3% and 5.8%). QCC.xjau-1DS-pre was detected at the heading stage and mapped 1 cM apart from QCC.xjau-1DS, indicating that they are likely the same QTL. QCC.xjau-3AL-pre was detected at the heading stage, explaining 2.5% of the phenotypic variation (Table 6). All QTLs related to CC were detected under normal irrigation treatment.

3.6. Candidate Gene Mining and Comparative Analysis of Predicted and Measured Values

On the basis of the flanking SNP markers of the QTL, the DNA sequences corresponding to these loci were searched in the Chinese Spring reference genome; 10 candidate genes were identified for LAI-related QTLs, and three candidate genes were identified for CC-related QTLs (

Table 7. The candidate genes for QTLs mapped with predicted and measured values of LAI and CC.

Source	Traits	QTL Name	Marker Interval	Physical Location (Mb)	Gene	Function Description
Predicted value	LAI	QLAI.xjau-2BL-pre.1	AX-89402509– AX-95134753	797.209401–797.210138	TraesCS2B01G623600	F-box family protein
	LAI	QLAI.xjau-2BL-pre.2	AX-95631217– AX-94667571	555.852677–555.853828	TraesCS2B01G391300	F-box family protein
	LAI	QLAI.xjau-3BL-pre	AX-179558689– AX-94903264	81.813957–81.818009	TraesCS3B01G115100	Kelch-like protein
	LAI	QLAI.xjau-5DL-pre	AX-94670615– AX-111652649	473.052014–473.053148	TraesCS5D01G408900	Peroxidase
	LAI	QLAI.xjau-5DL-pre	AX-94670615– AX-111652649	474.813752–474.814753	TraesCS5D01G411400	C2H2-type zinc finger protein
	CC	QCC.xjau-1DS-pre	AX-109983565– AX-89314186	382.257141–382.258109	TraesCS1D01G284100	GATA transcription factor
	CC	QCC.xjau-3AL-pre	AX-110449574– AX-109301028	711.402978–711.407172	TraesCS3A01G480600	BTB/POZ domain- containing protein
Measured value	LAI	QLAI.xjau-1BL	AX-111775778– AX-109367690	652.781869–652.786461	TraesCS1B01G427400	Wuschel homeobox protein
	LAI	QLAI.xjau-2BL.1	AX-111757968– AX-110468087	794.07325–794.082431	TraesCS2B01G617700	MYB-related transcription factor
	LAI	QLAI.xjau-2BL.2	AX-94407412– AX-111743644	691.187478–691.1904	TraesCS2B01G493800	BTB/POZ domain-containing protein
	LAI	QLAI.xjau-2BL.2	AX-94407412– AX-111743644	712.002368–712.004414	TraesCS2B01G517100	F-box family protein
	LAI	QLAI.xjau-5BS	AX-179559227– AX-179561411	19.455584–19.456986	TraesCS5B01G019700	F-box family protein
	CC	QCC.xjau-1DS	AX-110335177– AX-109983565	382.257141–382.258109	TraesCS1D01G284100	GATA transcription factor
	CC	QCC.xjau-1DS	AX-110335177– AX-109983565	392.69691–392.697515	TraesCS1D01G293600	Abscisic acid receptor

). Four candidate genes for LAI belonged to F-box family proteins. On chromosome 1DS, candidate gene TraesCS1D01G284100 encoding a GATA transcription factor was identified from the two closely linked QTLs QCC.xjau-1DS and QCC.xjau-1DS-pre according to predicted and measured CC values, indicating that the two loci are likely same.

4. Discussion

4.1. Spectral Reflectance Difference Analysis

Plant leaves play important roles in the process of photosynthesis and respiration [28]. The formation of a reasonable leaf area in a crop provides an area to guarantee receiving light energy. Chlorophyll is one of the important pigments for light absorption, and its content is significantly related to the photosynthesis rate and growth vigor [29]. Since colored substances, such as plant carotenoids, chlorophylls, and brown pigments, absorb sunlight in the wavelength range of 400–760 nm, the light spectral curve of reflectance in this range is low [30]. In this study, the reflectance of the five spectral bands on wheat decreased in the following order: NIR light > red-edge light > green light > red light > blue light. Because plants mainly absorb red light and blue light, their reflectance is the lowest. The reason for the rapid increase in the reflectance of light in the range from NIR to red-edge light is due to the fact that the colored substances (mainly chlorophyll) in wheat leaves can absorb most of the visible light, but it is difficult to absorb radiation with wavelengths above 700 nm [31].

In this study, we found that the NIR reflectance of wheat LAI under drought stress was higher at all tested growth stages, compared with the NIR reflectance under normal irrigation conditions. LAI is affected by plant height, leaf size, and degree of curling. Drought stress accelerates plant water loss and causes leaf curling and plant height shortening, which decreases the leaf coverage in the wheat population and increases the exposed area of the ground and infrared reflectance. In another aspect, the water in wheat leaves under normal irrigation conditions was higher than that under drought stress, because water can absorb light [32], which is also a reason that the NIR reflectance on wheat was higher under drought stress than under normal irrigation conditions. Compared with normal irrigation conditions, green and red light reflectance increased under drought stress. The reason is that drought stress accelerates plant water evaporation and chlorophyll degradation, reducing the leaf area and chlorophyll content for receiving visible light.

4.2. LAI and CC Modeling Verification

With the rapid application of UAV drone and remote sensing technology on agriculture, more imagery data of crops have been obtained, laying the foundation for investigating crop phenotypes in multiple growth periods and environments and accelerating crop genetic research [7]. Processing the imagery data to extract useful information of the corresponding phenotypes is critical. Machine learning (ML) methodologies can be used to effectively analyze and utilize information-rich datasets, as well as high-dimensional observation data [33]. They have been used in the analysis and modeling of remote sensing data, but the accuracies of different ML methods vary. In this study, we used the CART algorithm [26] to build a model through the inductive learning method. First, the partial sample data were used to construct the decision tree, then the remaining sample data were gradually added to adjust the tree, and finally the decision tree was optimized and completed [34]. According to the greedy algorithm, the whole sample data were used, and the statistical properties of all sample data were obtained. The above two steps effectively resisted the noise pollution of the results. The corresponding predicted values were obtained after modeling by the algorithm, and the predicted values and the measured values were validated by the cross-validation method (Table 2). The verification coefficient of determination R² of LAI under normal irrigation was 0.70–0.84, while that under drought stress was 0.64–0.78. The verification coefficient of determination R² for CC was relatively high and consistent under normal irrigation (0.876–0.883) and drought stress (0.86–0.88) conditions. This algorithm was more effective for the phenotype prediction of the normal irrigation treatment, with a higher prediction accuracy and reliability.

Various algorithms and models have been applied to predict wheat LAI and CC through collecting images with UAV and multispectral sensor technologies [12,35,36]. Gao et al. [36] used the UAV hyperspectral vegetation index to model 103 LAIs of multiple wheat single-growth periods, and the verification R² of the model was 0.783. Other reports of UAV-based multispectral estimated vegetation indices in wheat, as well as other crops, under drought conditions had similar coefficients [11]. The results in our current study demonstrate that both NDVI and GNDVI can be used to estimate LAI and CC in wheat from the flowering to maturity stages.

Although UAV remote sensing technology has been widely used to monitor crop traits recent years, most studies used small numbers of plant genotypes; the collection of high-throughput data in a large number of genotypes, such as a mapping population, is rarely reported. By combining the multispectral imaging with the vegetation index to construct related models in this study, we phenotyped LAI and CC in a mapping population of 309 RILs at three wheat growth stages under two different water treatments. A high correlation was obtained between the UAV imagery-predicted values and ground-measured values. The usefulness of this platform was further demonstrated by the QTL mapping of the drought tolerance traits in the wheat population.

4.3. QTLs and Candidate Genes Association with Predicted and Measured Values of LAI and CC

UAV is a promising platform to predict the timeseries development of crop canopies, and these data can be further used to understand the genetic basis of phenotypic variation [37]. In this study, we identified eight QTLs associated with LAI and two QTLs associated with CC. Three QTLs related to LAI were positioned on chromosome 2BL, and two of them (QLAI.xjau-2BL-pre.1 and QLAI.xjau-2BL.1) were 5 cM apart. For CC, two QTLs were identified on chromosome 1DS (QCC.xjau-1DS-pre and QCC.xjau-1DS), closely linked with 1 cM apart. The effects of these QTL clusters were detected from both predicted and ground-measured values, strengthening the QTL analysis. The positions of three QTLs on chromosome 2BL were so close that they are very likely to be the same locus, and this was also the case for two QTLs on 1DS. Similar findings were also reported by Marco et al. [38]; they proposed that, when two QTLs on the same chromosome of wheat are ≤15 cM apart, the two QTLs can be considered as the same locus. To infer whether there were common interacting QTLs or regions across genetic backgrounds, we compared the QTLs for LAI and CC with the QTLs reported in previous studies. Wang et al. [39] used an RIL population (137 lines) constructed from Shanghai 3/CBRD and Naxos combined with SSR markers to identify the flag leaf width-related QTL on chromosome 2B (90.8 cM), which is 10.2 cM away from the LAI locus QLAI.xjau-2BL-pre.1 detected on chromosome 2BL (101 cM) in this study. Another QTL (QSR.sicau-3B) [40] for spike formation rate (SR) located on chromosome 3B (47cM), which is 9 cM away from the LAI locus QLAI.xjau-3BS-pre detected on chromosome 3BS (56 cM) in the current study. Zhao et al. [41] mapped the RIL population (128 lines) from Ningchun 4 × Ningchun 27 with SSR markers and identified a QTL related to flag leaf length and width on chromosome 5D (107.4 cM), which is 0.6 cM away from QLAI.xjau-5DL-pre for LAI detected in the present study.

The CC-related QTL site QCC.xjau-3AL-pre (103cM) in this study is potentially novel. This QTL was detected in the interval of AX-110449574 and AX-109301028 on chromosome 3AL, explaining 2.5% of the phenotypic variation of CC. There is no CC-related QTL previously reported in this chromosomal region, indicating that this QTL is likely a new genetic locus regulating CC. Therefore, the combination of the UAV image technique with QTL mapping of the RIL population, as reported here, can greatly improve the identification of novel genomic factors underlying the genetic variation under normal irrigation and drought stress conditions and can be used in functional annotation of the wheat genome.

According to the chromosome position and genomic sequences of the QTLs on the reference wheat genome, the potential candidate genes of the QTLs for LAI and CC were annotated in this study. Among 10 candidate genes, TraesCS2B01G623600 and TraesCS2B01G391300 both encode F-box family proteins, which are involved in the regulation of plant stress resistance responses [42]. TraesCS3B01G115100 encodes a Kelch-like protein, which is one of the C-terminal domains of F-box proteins and mediates the F-box protein response to drought stress [43]. TraesCS5D01G408900 encodes peroxidase (POD), which can remove H₂O₂ produced in plants to reduce the H₂O₂ damage on plants [44]. POD activity is also significantly related to the drought resistance of plants [45]. TraesCS5D01G411400 encodes a C₂H₂ zinc finger protein, which can slow down the decrease in leaf water content and increase in antioxidant enzyme activity, leading to a decrease in active oxygen content in plants and an improvement in plant resistance to drought stress [46]. TraesCS1D01G284100 encodes the GATA transcription factor, which participates in the formation of plant chloroplasts [47]. TraesCS3A01G480600 encodes a BTB/POZ domain junction protein, which participates in plant-related signal networks and responds to adversity stress [48,49]. These genes have functions in plant resistance under environmental stress.

5. Conclusions

In this study, UAV multispectral images and manual measurements of LAI and CC from a population of 309 wheat RILs were collected in multiple growth periods under normal irrigation and drought stress conditions. GNDVI and NDVI were used to predict the LAI and CC through the CART algorithm. The predictive model for LAI and CC was judged though cross-validation with good estimations of LAI and CC at three wheat growth stages under two water supply treatments, with R² ranging from 0.64 to 0.84 for LAI and from 0.86 to 0.88 for CC, RMSE ranging from 0.30 to 0.52 for LAI and from 0.94 to 1.05 for CC, and RE ranging from 0.06 to 0.18 for LAI and from 0.01 to 0.02 for CC. In addition, high correlations (0.87 to 0.91) were found between predicted and measured values for LAI and CC at each growth stage of wheat under different water treatments. These results indicate that the predictive model for LAI and CC had high accuracy and reliability. According to the high-density SNP map of the RIL population, eight QTLs of LAI and three QTLs of CC were mapped, with each QTL explaining 2.5% to 13.8% of the phenotypic variation for LAI and 2.5% to 5.8% of the phenotypic variation for CC. The potential candidate genes of these QTLs were identified in the corresponding chromosomal regions in the wheat reference genome. In addition to three QTLs detected by both measured and estimated values, three QTLs were detected by estimated values. These findings indicate that the application of UAV remote sensing technology can significantly improve the efficiency of high-throughput phenotyping and crop genetic studies.