Carbon Isotope Composition and the NDVI as Phenotyping Approaches for Drought Adaptation in Durum Wheat: Beyond Trait Selection

High-throughput phenotyping platforms provide valuable opportunities to investigate biomass and drought-adaptive traits. We explored the capacity of traits associated with drought adaptation such as aerial measurements of the Normalized Difference Vegetation Index (NDVI) and carbon isotope composition (δ13C) determined at the leaf level to predict genetic variation in biomass. A panel of 248 elite durum wheat accessions was grown at the Maricopa Phenotyping platform (US) under well-watered conditions until anthesis, and then irrigation was stopped and plot biomass was harvested about three weeks later. Globally, the δ13C values increased from the first to the second sampling date, in keeping with the imposition of progressive water stress. Additionally, δ13C was negatively correlated with final biomass, and the correlation increased at the second sampling, suggesting that accessions with lower water-use efficiency maintained better water status and, thus, performed better. Flowering time affected NDVI predictions of biomass, revealing the importance of developmental stage when measuring the NDVI and the effect that phenology has on its accuracy when monitoring genotypic adaptation to specific environments. The results indicate that in addition to choosing the optimal phenotypic traits, the time at which they are assessed, and avoiding a wide genotypic range in phenology is crucial.


Introduction
Global warming makes agricultural productivity more unpredictable, especially in the semiarid regions of the world. Indeed, in the southern and eastern parts of the Mediterranean Basin where durum wheat is the main annual crop and of great cultural importance, the frequency and severity of drought events (understood as the combination of water and heat stresses) show an upward trend [1]. Drought is a major factor limiting wheat yield, but, despite decades of research, it remains a challenge for plant breeders [2]. Drought stress significantly affects growth and photosynthetic plant performance and, consequently, agronomical grain yield components such as the number of fertile spikes per unit Under field conditions, the δ 13 C has been negatively associated with biomass and yield, not only under well-watered conditions [36] but also in environments exhibiting moderate to severe drought [32,37]. As a lower δ 13 C implies a lower WUE under drought conditions, the yield advantage may result from constitutive (e.g., phenology, plant size, etc.) or stress-adaptive (e.g., osmotic adjustment, root architecture, etc.) traits that allow the crop to avoid the negative effects of drought stress. However, positive relationships between δ 13 C and biomass and yield [38][39][40] have also been reported, particularly under severe drought conditions [30,35,41,42]. δ 13 C is a good indicator of the water strategy that crop plants adopt in the field and is a desirable trait for breeding programs [3] because it integrates plant transpiration [43] and the photosynthetic performance of the plant over time [31] and can be easily sampled in a large number of genotypes [35] with high repeatability [43]. In addition, although the application of δ 13 C has been mainly related to plant water status, it has been used as a secondary trait for selecting genetic loci associated with drought tolerance and also as a measure of increasing yield under good agronomical conditions in wheat [44] and other cereals [45,46].
Although this is often forgotten, timing is key for phenotyping as well as defining derived actions such as selection indices and ideotypes and even the genetic dissection and formulation of quantitative trait loci (QTL) [47]. The objective of the present study was to explore the capacity to predict genotypic variability in green biomass during grain filling by using aerial measurements of NDVI at different phenological stages and evaluate the δ 13 C in the flag leaves sampled initially under well-watered conditions and, subsequently, under drought. In addition to the assessment of these phenotypic traits, another important objective was to determine the best time to measure such traits in terms of crop phenology and growing conditions. The study was performed in a panel of 248 genotypes of durum wheat grown in the Maricopa phenotyping platform, one of the best-integrated agricultural technology platforms and the largest open-air phenotyping robotic scanner in the world. The rationale of our work was to illustrate the confounding effects phenology may have when phenotyping wheat. We addressed phenology under two premises. Firstly, genotypic variability in phenology, which results in different genotypes being sampled at different phenological stages at a given moment. This may affect, for example, the sign and strength of the relationships between δ 13 C and plant biomass. Secondly, the phenological stage at which plants are measured may strongly affect the performance of any phenotypic trait. As an example, the sign and strength of the relationship between NDVI and crop biomass may be affected by the crop stage when NDVI is measured. We might even add a third factor, which may somehow interact with the previous ones: Besides phenology, the growing conditions during evaluation will likely influence phenotypic values. This is illustrated again by the relationship between carbon isotope composition and biomass, where the strength of the relationship changes from well-irrigated plants to plants last irrigated two weeks before, a scenario frequently experienced by wheat under Mediterranean conditions where terminal drought (i.e., during grain filling, the last part of the cropping cycle) occurs. In other words, even analyzing the same organ (the flag leaf in our study), the information derived from δ 13 C is not similar but is affected by growing conditions.

Germplasm Used and Experimental Conditions
The field experiment was conducted during 2016 at the Maricopa Agricultural Center (33.070 • N, long. 111.974 • W, elev. 360 m) on a Casa Grande sandy loam soil (fine-loamy, mixed, superactive, hyperthermic Typic Natrargids). The plant material included 248 accessions of durum wheat (Triticum turgidum L. subsp. durum (Desf.) Husn.) from the association mapping population assembled by Maccaferri and co-workers at the University of Bologna (UNIBO), namely the UNIBO-Durum Panel representing an Elite Durum Panel (EDP) that comprises a large portion of the genetic diversity present in the most important improved durum wheat gene pools. EDP encompassed Mediterranean-adapted accessions from different private and public international breeding programs and included elite cultivars and lines originally collected from Mediterranean countries, Mexico, Australia, Europe, and the Southwestern USA [23]. The accessions were planted on 20 December 2016 according to a Randomized Complete Block Design (RCBD) with two replicates and border plots (consisting of durum wheat cv. Orita). Each accession was evaluated in two-row plots (3.5 m long, 0.76 m apart) with a final density of 22 plants m −2 . Prior to sowing, the experiment received 112 kg ha −1 of nitrogen and 56 kg ha −1 of phosphorus (P 2 0 5 ). Phenology was assessed visually using Zadoks decimal code and flowering time was assessed as the number of days to reach anthesis between the latest and earliest accessions, assigning zero to the earliest accession (Supplementary Table S1). Disease and insect pest pressure were negligible throughout the crop cycle.
Irrigation started 28 days after sowing, with the initial irrigation provided by sprinklers and further irrigations delivered via subsurface drip irrigation. The field received a total of 115 mm of water distributed from 11 irrigations. Drip irrigation was stopped on 16 March 2017, which corresponded to growth stage (GS) 45 (averaged across the 248 accessions) according to the Zadoks decimal code [48], and from that date the accessions were subjected to a progressive drought stress until 3-4 April 2017 (averaged GS 62) when plants were harvested to measure biomass.
The accumulated precipitation from planting to crop harvest was 50 mm, while the accumulated potential evapotranspiration from planting to the last irrigation was 222 mm and from planting until crop harvest was 320 mm (Figure 1). Environmental conditions during growth are detailed in Figure 1.  Table S1). Disease and insect pest pressure were negligible throughout the crop cycle. Irrigation started 28 days after sowing, with the initial irrigation provided by sprinklers and further irrigations delivered via subsurface drip irrigation. The field received a total of 115 mm of water distributed from 11 irrigations. Drip irrigation was stopped on 16 March 2017, which corresponded to growth stage (GS) 45 (averaged across the 248 accessions) according to the Zadoks decimal code [48], and from that date the accessions were subjected to a progressive drought stress until 3-4 April 2017 (averaged GS 62) when plants were harvested to measure biomass.
The accumulated precipitation from planting to crop harvest was 50 mm, while the accumulated potential evapotranspiration from planting to the last irrigation was 222 mm and from planting until crop harvest was 320 mm (Figure 1). Environmental conditions during growth are detailed in Figure  1.

NDVI Measurements
UAV-based NDVI data were extracted from georeferenced orthomosaic Geostationary Earth Orbit Tagged Image File Format (GeoTIFF) files generated from images captured by autopiloted flights of a Parrot Sequoia (Parrot, Paris, France) multi-spectral camera (referred to as SE), carried on an eBee (SenseFly, Lausanne, France) fixed-wing aircraft, and a MicaSense RedEdge (MicaSense, Seattle, WA, USA) multi-spectral camera (referred to as RE), carried on a hexacopter. Differences between the two multispectral cameras were mainly based on band centers and bandwidths as described in [23]. The purpose for using cameras with slightly different band centers and bandwidths was to explore the phenotypic performance of the NDVI generated with the two multispectral sensors. The flights were performed at a height of 40 m to 42 m from ground level where it was possible to achieve a ground sampling distance of 4.4 cm/pixel for the Sequoia and ~3 cm/pixel for the RedEdge.

NDVI Measurements
UAV-based NDVI data were extracted from georeferenced orthomosaic Geostationary Earth Orbit Tagged Image File Format (GeoTIFF) files generated from images captured by autopiloted flights of a Parrot Sequoia (Parrot, Paris, France) multi-spectral camera (referred to as SE), carried on an eBee (SenseFly, Lausanne, France) fixed-wing aircraft, and a MicaSense RedEdge (MicaSense, Seattle, WA, USA) multi-spectral camera (referred to as RE), carried on a hexacopter. Differences between the two multispectral cameras were mainly based on band centers and bandwidths as described in [23]. The purpose for using cameras with slightly different band centers and bandwidths was to explore the phenotypic performance of the NDVI generated with the two multispectral sensors. The flights were performed at a height of 40 m to 42 m from ground level where it was possible to achieve a ground sampling distance of 4.4 cm/pixel for the Sequoia and~3 cm/pixel for the RedEdge. Latvia) for the RedEdge camera. An 80% image overlap along the flight corridors was achieved for all flights. Both the Sequoia and RedEdge cameras use global shutters. Orthomosaics for each camera band were generated with Pix4DMapperPro desktop software (Pix4D SA, Prilly, Switzerland, http://pix4d.com), as explained elsewhere [23]. Real-time kinematic (RTK) survey precision was used to geolocate six to eight ground control points (GCPs) and to georeference the orthomosaics. At the beginning of each flight, manufacturer-supplied reflectance panels were imaged in order to calibrate camera images. QGIS software version 2.18.3 (QGIS, Zurich, Switzerland, http://www.qgis.org) was used to create plot-level NDVI means from UAVs. Spatial adjustment analysis of the raw NDVI plot data from aerial platforms was conducted with the lme4 package (r-project) and custom R scripts using a mixed procedure including random effects of columns and rows and a moving mean of variable size for optimizing spatial adjustment. Annotated, single-plot polygons enclosed by shape files were generated with an R (r-project.org) script. Based on RTK survey-grade measuring devices, shape files with GCPs as features (points) were also employed. NDVI orthomosaics generated from Pix4D and the GeoTIFF were combined with the GCP shape files and the plot polygon in a single QGIS project. Verification of proper geolocations of the Pix4D orthomosaics was performed by visually confirming the alignment of the visible GCPs with the corresponding points in the feature shape file. The Zonal Statistics function in QGIS was used to generate NDVI plot means. NDVI was measured at multiple times using RE and SE sensors. NDVI based on SE sensor was measured on 13 and 21 February and 3 and 13 March 2017, whereas NDVI based on RE sensor was measured at two different time points on 21 and 28 March 2017.

Dry Biomass Evaluation
At the end of the field trial, plants within the entire two-row plots were harvested mechanically (Carter Manufacturing Company, Brookston, IN, USA) and the weight was adjusted to 0% moisture. Plant moisture content (%) at harvest was estimated from a subsample of biomass from 2-3 plants, either placed directly in a drying oven or stored temporarily in an uncooled greenhouse that reached a diurnal high temperature of 60 • C before being transferred to an oven at 60 • C for final drying.

Total Nitrogen and Carbon Content and Carbon Isotope Composition
On 17 March (the day after the final irrigation, i.e., before water stress induction) and 30 March (about two weeks after drought induction) six representative leaves (the flag leaf, if present, or the most recent fully developed leaf) from main tillers were collected per plot, oven dried at 70 • C for 48 h, weighed, and finely ground for carbon isotope analysis. The stable carbon isotope composition (δ 13 C) together with the total carbon and nitrogen concentrations of the control (17 March) and stressed (30 March) leaves were determined using an elemental analyzer (EA; Flash 1112 EA, Thermo Fisher Scientific, Bremen, Germany) coupled with an isotope ratio-mass spectrometer (IRMS; Delta C with CONFLO III interface, Thermo Fisher Scientific, Bremen, Germany) operating in continuous-flow mode in order to determine the stable carbon ( 13 C/ 12 C) isotope ratios of the same samples. Samples of approximately 1 mg of total dry matter (DM) for the dry biomass and reference materials were weighed into tin capsules, sealed, and then loaded into an automatic sampler (Thermo Fisher Scientific, Bremen, Germany) before EA-IRMS analysis. The 13 C/ 12 C ratios of the plant material were expressed in δ notation [49]: δ 13 C = ( 13 C/ 12 C) sample /( 13 C/ 12 C) standard − 1, where "sample" refers to plant material and "standard" to international secondary standards of known 13 C/ 12 C ratios (International Atomic Energy Agency (IAEA) CH7 polyethylene foil, IAEA CH6 Sucrose, and the United States Geological Survey (USGS) 40 l-glutamic acid) calibrated against Vienna Pee Dee Belemnite calcium carbonate with an analytical precision (SD) of 0.10% . Total carbon and nitrogen contents were expressed as a percentage of the dry matter (%). Measurements were carried out at the Scientific Facilities of the University of Barcelona. The δ 13 C of the samples collected before drought induction (17 March) and two weeks after (30 March) are referred to as δ 13 C control and δ 13 C stress, respectively.

Statistical Analysis
NDVI measurements, phenology, δ 13 C, and total carbon and nitrogen contents were subjected to one-way analysis of variance (ANOVA) using the general linear model to calculate the effect of crop growth prior to anthesis and progressive drought from anthesis to mid grain filling. For the NDVI, the date of measurement was included as a fixed factor in the ANOVA, whereas in the δ 13 C it was the absence of irrigation. Means were compared by Tukey's honestly significant difference (HSD) test. A bivariate correlation procedure was constructed to analyze the relationships between the measured traits and plant dry matter. Descriptive statistics of agronomic variables were estimated, and analysis of principal components (PCA) and stepwise analysis were calculated using the SPSS 18.0 statistical package (SPSS, Chicago, IL, USA).

General Characteristics of the Elite Durum Panel (EDP)
The EDP ranged from 0 to 35 days between extreme accessions in terms of date of anthesis. The dry biomass harvested two weeks after the last irrigation and averaged across the genotypes of the EDP was 2.67 Mg ha −1 ( Table 1). The Normalized Difference Vegetation Index (NDVI) increased as crop development progressed, reaching its maximum around anthesis (Table 1). Thus, the lowest NDVI value, averaged across the EDP, was 0.548 (measured with the SE sensor) and was measured on 13 February, corresponding on average to GS 31 on the Zadoks phenological scale [48]. The highest NDVI average value was 0.879, using the SE sensor, which corresponded to a value of 0.829 (when measured with the RE sensor). This measurement took place on 21 March, five days after the last irrigation. Twelve days after the last irrigation (28 March), which corresponded on average to the GS 58 stage, the NDVI (RE) decreased to 0.772. Control plants showed lower values for δ 13 C (more negative δ 13 C) compared to water-stressed plants ( Table 2). Total C and N concentrations were higher after drought induction, although in the one-way analysis of variance the effect of drought was more significant in C than N, as observed from the higher values of the mean squares in the former (Table 2). Table 2. Average (AVG) and standard error (SE) for the stable isotope composition of carbon (δ 13 C) and total carbon (Total C) and nitrogen (Total N) concentrations of the dry matter of aerial biomass sampled before (control) and after (stress) drought was imposed on the 248 EDP accessions. The outputs of the ANOVA for treatment effect (water availability) are included to the right of the variable. Level of significance ** p < 0.01, *** p < 0.001.

Control
Stress

Relationship between Studied Parameters and Biomass
Principal component analysis (PCA) showed that the δ 13 C of stressed plants was negatively correlated with plant biomass (Figure 2, upper panel), whereas the δ 13 C of control plants did not show a clear relationship to the biomass (because a perpendicular angle was observed between the two parameters). The δ 13 C (either control or stressed plants) showed a positive relationship with the flowering time, suggesting a sensitivity of δ 13 C to the wide range of phenologies present in the EDP (Figure 2, upper panel). Biomass was also sensitive to phenology, as observed by its negative relationship to flowering time in the PCA, with early accessions being related to higher biomass, whereas late accessions showed lower biomass (Figure 2, upper panel).
Considering the NDVI measurements at different time points, the NDVI (SE) measured at early time points during crop development and before drought induction (13 and 21 February) showed a closer and positive relationship to the biomass (Figure 2, lower panel). In contrast, the NDVI measured after the last irrigation (i.e., under increasing drought conditions) was poorly correlated with biomass, as observed by the perpendicular angle between the biomass and the NDVI (RE) measured at the later time points (21 and 28 March). However, the NDVIs measured after drought induction were closely correlated with flowering time in the PCA (Figure 2, lower panel), suggesting that these NDVIs primarily represented the 'greenness' of late-phenology accessions, which are associated with a lower biomass (as observed by the negative relationship between flowering time and biomass). In fact, when considering the accessions for each flowering time independently, no correlations were observed between biomass and NDVI in either the accessions with early or late phenologies (Supplementary  Table S2). However, the NDVIs measured in accessions with early phenology (10 to 15 days to flowering time) tended to be positively correlated with biomass (meaning that early phenologies are associated with higher biomass, as observed in the PCA), whereas correlations in accessions with delayed phenologies (16 to 35 days to flowering time) tended to be negative (delayed phenologies are associated with a lower biomass).
Moreover, in accordance with the PCA, the linear regressions between the NDVI measured before anthesis in the absence of water stress (i.e., control plants) and biomass (Figure 3, upper panel) were positive and significant: r = 0.13 (p < 0.05) and r = 0.18 (p < 0.01), for NDVI (SE) measurements performed on 13 and 21 February, respectively. However, no correlation with biomass was observed when the NDVI was measured five days after the last irrigation (21 March), and the relationship tended to be negative and almost significant (p < 0.05) when the NDVI was measured nearly two weeks after the last irrigation (28 March). Moreover, linear correlation coefficients between flowering time and the NDVI changed from being non-significant before anthesis to positive and highly significant (p < 0.001) around anthesis, and further increased as drought advanced. Therefore, the positive relationship between NDVI and flowering time during the late stages of the crop (13, 21, and 28 March) was substantiated by the 'greenness' of the most delayed accessions (Figure 3, middle panel). Similarly, the δ 13 C of stressed plants showed positive correlations with NDVI (p < 0.001) and their strength increased progressively as NDVI was measured through drought induction (13,21,and 28 March), whereas correlation between the δ 13 C of control plants and the NDVI did not clearly increase over time in the presence of drought (on average r 0.12, p < 0.05) (Figure 3, lower panel).
The δ 13 C in the stressed plants seemed to be the best parameter correlated with biomass as it was the first variable chosen in the multilinear model (stepwise analysis) combining the whole set of genotypes with all flowering times (global) and considering the biomass as the dependent variable (Table 3) and all NDVI and δ 13 C measurements as independent variables. However, when subsets of genotypes were considered on the basis of their flowering time alone, the results were different. When a flowering time range of 10 to 15 days was considered, the δ 13 C of the control plants was the first variable chosen by the multilinear model, whereas the δ 13 C of stressed plants was the first variable chosen for the subset of genotypes with a flowering time of 17 days (Table 3). For the subset of late genotypes, with a flowering time of 25 days, the variables chosen by the multilinear model were NDVIs measured after drought (28 March, followed by 21 March) ( Table 3). The δ 13 C in the stressed plants seemed to be the best parameter correlated with biomass as it was the first variable chosen in the multilinear model (stepwise analysis) combining the whole set of genotypes with all flowering times (global) and considering the biomass as the dependent variable (Table 3) and all NDVI and δ 13 C measurements as independent variables. However, when subsets of genotypes were considered on the basis of their flowering time alone, the results were different.    Stepwise analysis for the entire 248 panel of durum wheat elite advanced lines and cultivars from around the world with biomass as a dependent variable, and the following independent variables: Carbon stable isotope composition (δ 13 C) and total nitrogen (Total N) concentration of the flag leaf dry matter sampled before stress conditions were imposed in the field (control) and two weeks after the last irrigation (stress), flowering time, and the Normalized Difference Vegetation Index (NDVI) using Sequoia (SE) and RedEdge (RE) sensors measured on different days after planting. The "global" stepwise analysis represents values obtained by considering accessions from all flowering times together. The 10, 12, 15, 17, and 25 days of flowering represent the stepwise analyses obtained with accessions within each flowering time classification.

Effect of Phenology on the Relationship of δ 13 C to Biomass
Biomass was negatively correlated with the δ 13 C of stressed and control plants ( Figure 4). However, the correlation between biomass and the δ 13 C was slightly improved in stressed plants (r = 0.54, p < 0.001) compared to control plants (r = 0.46, p < 0.01). Even so, considering the sensitivity of the δ 13 C to phenology (as observed by the close relationship between flowering time and the δ 13 C in the PCA, see above), correlations among biomass and δ 13 C for the subset of genotypes with specific flowering times were studied ( Figure 5 and Supplementary Table S3). Such correlations still remained negative and significant at different flowering times (Supplementary Table S3), especially in accessions with earlier phenologies (flowering times of 10 and 17 days). Moreover, days to flowering was positively correlated with the δ 13 C of control plants and even more strongly with the δ 13 C of stressed plants ( Figure 6). Thus, considering that one of the driving factors linked to the δ 13 C is the intrinsic photosynthetic capacity of the leaf, total N (used as a proxy for the intrinsic photosynthetic capacity) was correlated with δ 13 C (Figure 7). Both parameters showed no significant correlations in the control or stressed plants.

Effect of Phenology on the Relationship of δ 13 C to Biomass
Biomass was negatively correlated with the δ 13 C of stressed and control plants (Figure 4). However, the correlation between biomass and the δ 13 C was slightly improved in stressed plants (r = 0.54, p < 0.001) compared to control plants (r = 0.46, p < 0.01). Even so, considering the sensitivity of the δ 13 C to phenology (as observed by the close relationship between flowering time and the δ 13 C in the PCA, see above), correlations among biomass and δ 13 C for the subset of genotypes with specific flowering times were studied ( Figure 5 and Supplementary Table S3). Such correlations still remained negative and significant at different flowering times (Supplementary Table S3), especially in accessions with earlier phenologies (flowering times of 10 and 17 days). Moreover, days to flowering was positively correlated with the δ 13 C of control plants and even more strongly with the δ 13 C of stressed plants ( Figure 6). Thus, considering that one of the driving factors linked to the δ 13 C is the intrinsic photosynthetic capacity of the leaf, total N (used as a proxy for the intrinsic photosynthetic capacity) was correlated with δ 13 C (Figure 7). Both parameters showed no significant correlations in the control or stressed plants.  . Linear regression of the relationship between the plant biomass at mid grain filling and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control) and two weeks after the last irrigation (δ 13 C stress). Level of significance, ** p < 0.01, *** p < 0.001.
Agronomy 2020, 10, x FOR PEER REVIEW 12 of 20 Figure 5. Linear regression of the relationship between the plant biomass at mid grain filling and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control, right panels) and two weeks after the last irrigation (δ 13 C stress, left panels) across the subset of genotypes exhibiting the same date of anthesis: 10 days after the earliest genotype (n = 76) and 17 days (n = 63). Level of significance, ns p > 0.05, ** p < 0.01, *** p < 0.001. Figure 5. Linear regression of the relationship between the plant biomass at mid grain filling and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control, right panels) and two weeks after the last irrigation (δ 13 C stress, left panels) across the subset of genotypes exhibiting the same date of anthesis: 10 days after the earliest genotype (n = 76) and 17 days (n = 63). Level of significance, ns p > 0.05, ** p < 0.01, *** p < 0.001. Figure 5. Linear regression of the relationship between the plant biomass at mid grain filling and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control, right panels) and two weeks after the last irrigation (δ 13 C stress, left panels) across the subset of genotypes exhibiting the same date of anthesis: 10 days after the earliest genotype (n = 76) and 17 days (n = 63). Level of significance, ns p > 0.05, ** p < 0.01, *** p < 0.001.   . Linear regression of the relationship between the total nitrogen (Total N) concentration and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control, right panels) and after two weeks of the last irrigation (δ 13 C stress, left panels). Relationships were studied by considering the subset of genotypes with a flowering time of 10 (n = 76) and 17 days (n = 63). Level of significance, ns p > 0.05.

Effect of Drought in the EDP Based on the NDVI and δ 13 C
To the best of our knowledge, this is the first study that estimates crop growth capacity, in terms of biomass at mid grain filling, in a panel of hundreds of durum wheat accessions by simultaneously measuring the NDVI via UAV-based remote sensing platforms and the δ 13 C at a plant level and using these parameters as proxies for plant growth, stay-green, and performance under different water conditions. In the present study, the UAV-based NDVI measured at different time points in the crop cycle reflected the increase in biomass during crop development until anthesis and the degree of staygreen during grain filling in response to progressive drought ( Table 1). Regardless of the multispectral imager used (RE or SE), the NDVI increased from the early stages of crop growth along with increases in the green biomass and the ability of the crop to intercept radiation [50] and peaked Figure 7. Linear regression of the relationship between the total nitrogen (Total N) concentration and the stable carbon isotope composition (δ 13 C) in the flag leaf dry matter sampled before stress conditions were imposed (δ 13 C control, right panels) and after two weeks of the last irrigation (δ 13 C stress, left panels). Relationships were studied by considering the subset of genotypes with a flowering time of 10 (n = 76) and 17 days (n = 63). Level of significance, ns p > 0.05.

Effect of Drought in the EDP Based on the NDVI and δ 13 C
To the best of our knowledge, this is the first study that estimates crop growth capacity, in terms of biomass at mid grain filling, in a panel of hundreds of durum wheat accessions by simultaneously measuring the NDVI via UAV-based remote sensing platforms and the δ 13 C at a plant level and using these parameters as proxies for plant growth, stay-green, and performance under different water conditions. In the present study, the UAV-based NDVI measured at different time points in the crop cycle reflected the increase in biomass during crop development until anthesis and the degree of stay-green during grain filling in response to progressive drought ( Table 1). Regardless of the multispectral imager used (RE or SE), the NDVI increased from the early stages of crop growth along with increases in the green biomass and the ability of the crop to intercept radiation [50] and peaked (0.879, i.e., close to saturation) shortly after anthesis and five days after the last irrigation (21 March). Such peak on NDVI values has been observed in the past as a consequence of canopy closure [51,52]. In fact, on 21 March, the NDVI was recorded simultaneously with RE and SE multispectral cameras and the different band centers and bandwidths of the cameras did not significantly affect the NDVI correlation with biomass or the other traits, and this suggests that both sensor types produced comparable NDVIs. In other words, although absolute NDVI values were in some cases statistically different, pattern and significance of the relationships of NDVI with biomass and other traits were not affected by the sensor used. Fourteen days after anthesis (28 March), progressive canopy senescence (inferred from the decrease of NDVI values and the visual, on-ground observation of the lower plant leaves) during grain filling probably affected NDVI by reducing the reflectance in the near-infrared while the reflectance in the red part of the spectrum increased, most likely due to cell degeneration and chlorophyll loss [14,16,50] as well as decreases in cell water content [53]. Moreover, water-limiting conditions also affected the δ 13 C (Table 2). Thus, reductions in soil water content and increases in vapor pressure deficit (VPD) contributed to differential water use [35] increasing the δ 13 C of the plant biomass relative to control in our study as the season progressed. Higher δ 13 C values in the plant biomass under drought conditions were probably related to low stomatal conductance, causing a lower C i /C a , which is frequently associated with an increased WUE [30,33]. Nitrogen concentration values of the flag leaves were very high at the two sampling dates and even increased slightly after two weeks from the last irrigation, which suggests that senescence induced by the treatment was not severe and did not affect the upper leaves. In that sense, the decreased NDVI values observed two weeks after the last irrigation relative to the values obtained five days after the last irrigation might have an alternative explanation. A leaf-rolling effect resulting from water deficit stress [23,54] could have been an important factor that decreased the NDVI, as observed from the negative relationship between the NDVI and leaf rolling in a previous study in the same trial [23]. A methodological limitation of the NDVI is that this vegetation index tends to saturate beyond values of 0.65-0.70, which would be the case for both NDVI measurements taken on 21 sand 28 March. However, notwithstanding this limitation, a small decrease was recorded between five days and two weeks after the last irrigation, suggesting plants were progressively experiencing senescence, leaf rolling, or a combination of both factors, after last irrigation.

Biomass Predictions in Early Stages of Crop Development in the EDP Based on the NDVI
NDVI measurements performed during the relatively early stages of crop development before anthesis and full water conditions were positively correlated with biomass at harvest (mid grain filling). This is in keeping with previous studies where NDVIs measured in the early stages of crop development have been positively correlated with biomass [55]. This is not trivial since biomass is correlated with yield and, in turn, yield is strongly correlated with the number of grains per area, a component that is basically determined in the early stages of crop development [56]. Thus, the NDVI values, which are usually below the threshold for saturation when measured early during the crop cycle, make this vegetation index a good proxy of green biomass [10,11,55,57,58]. Therefore, differences in canopy development estimated during the early stages of crop growth according to the NDVI may help to forecast crop growth capacity [14,55,58] and thus yield. The positive relationship between biomass and NDVI was strongest early in the crop development (Figure 3), suggesting that accessions with greater canopy development before anthesis were the most productive in terms of biomass at mid grain filling.

Biomass Prediction During the Late Stages of Crop Development Based on the NDVI
Unlike the relationship of biomass to the NDVI measured in the early stages of crop growth, the relationship between the NDVI measured after drought induction and biomass was poor and even tended to be negative (Figure 3). Although it was not the objective of our study, maximum growth in wheat is usually achieved around anthesis and afterwards the canopy starts to senesce, more or less quickly, depending on the growing conditions, and starting with the older leaves in the lower part of the plant. In accordance with this, after two weeks of progressive drought induction, the NDVI measured at mid grain filling decreased compared to values observed just five days after the last irrigation. However, the negative correlation between the NDVI measured after drought induction versus biomass and the positive relationship between the NDVI measured after drought versus flowering time (Figures 2 and 3) led to other interpretations of our results. Thus, genotypes with delayed phenology are characterized by more vegetative growth, which accounts for the higher NDVI of such accessions, while earlier genotypes are characterized by less vegetative growth [59]. The larger transpirative biomass of the late accessions may lead to faster water exhaustion and, so, to a worse water status. These results suggest that, after drought, late phenologies with lower biomass (although greener than in early phenologies) were probably associated with a lower proportion of stems and spikes relative to leaves, hence, affecting plant water status negatively due to a faster exhaustion of soil moisture in comparison to early flowering genotypes. This would agree with the more positive δ 13 C values of the late phenology genotypes (Figures 2 and 3).
Nevertheless, while the general relationships between NDVI and biomass are well established, selection of the correct crop phenological stage for NDVI monitoring of genotypic adaptation to specific environments will depend on the growing conditions and also the genotypic range of phenologies. For example, contrary to positive correlations between NDVIs measured after anthesis and yield in bread [16] and durum wheat [14,60], the range in phenology across the panel of genotypes tested was negligible.

Biomass Predictions Based on δ 13 C
The accessions with higher biomass were those with better water status (and lower WUE), as observed by the negative correlation between biomass versus δ 13 C in both stress and control conditions (Figure 4). In addition, positive correlations between the δ 13 C and flowering time ( Figure 6) supported the idea of better water status in accessions with relatively early phenologies (e.g., flowering times of 10 to 17 days) in contrast to the worse water status (as shown by the higher δ 13 C values) of accessions with late phenology (e.g., with a flowering time of 35 days). These results suggest that the more positive δ 13 C values in late phenologies are driven by a decrease in stomatal conductance due to terminal stress [61]. However, when the potentially disturbing effect of phenology was removed by analyzing the relationship between δ 13 C and biomass within subsets of genotypes with similar dates of anthesis ( Figure 5), relationships were in many cases still negative and significant, suggesting genotypes with higher δ 13 C (and, thus, higher WUE) were those which accumulated less biomass. Negative genotypic relationships of δ 13 C with canopy biomass and even more with grain yield have been reported before [3,38,47,62]. Accordingly, an extensive review on wheat and other herbaceous crops [63] have highlighted the positive association between higher stomatal conductance and a more productive genotype. In our study, the increase in δ 13 C due to stomatal closure was supported by the absence of a correlation between δ 13 C and total N either within subsets with the same date of anthesis (Figure 7) or across the whole set. The absence of a correlation between δ 13 C and total N discarded the intrinsic photosynthetic capacity of the plant as the main cause of δ 13 C enrichment [64] in response to elevated N fertilization and differences in phenology and senescence during grain filling [42]. While a positive relationship between days to anthesis and the δ 13 C of mature kernels in wheat has been reported [32,41,65], our study is the first to show that this relationship is already present in the δ 13 C of the flag leaves. Genotypes with lower δ 13 C values were likely associated with greater transpirative rates that were triggered by an increase in stomatal conductance. Thus, early flowering accessions that escaped from terminal drought maintained higher levels of stomatal conductance (which decreased δ 13 C), being the most productive in terms of biomass. This finding suggests that under well-watered conditions it is worthwhile for the plant to invest in growth rather than survival [8,37] and to ensure the maximization of production per unit of transpiration [66,67]. In fact, the maximization of productivity of green biomass occurred not only in well-watered plants but also after drought because the correlation between biomass and δ 13 C under drought was stronger than in the control and remained negative ( Figure 3). Therefore, our results indicate that even under (moderate) water-limiting conditions, the plant strategy is likely to have an impact on growth rather than survival mechanisms [37]. More importantly, our results also suggest that environmental factors are not the only basis of such a relationship. While the positive relationship between the date of anthesis and the δ 13 C of mature kernels under Mediterranean conditions can be related to the later genotypes being exposed to more severe drought conditions during grain filling [41], in our study the δ 13 C was analyzed in flag leaves sampled simultaneously, and they were, therefore, exposed to the same environmental conditions. In fact, in the absence of water stress, a lower (i.e., more negative) δ 13 C has been associated with faster growth [37,68]. The mechanisms underlying the observed response behind this are unclear but different avenues should be explored in the future. For example, higher stomatal conductance as an indicator of better water status [30,33,35,37,41] or thinner leaves associated with faster growth [37,68] may be related to a lower δ 13 C.

Concluding Remarks on How and When to Use Phenotyping Traits
Although phenology affected the δ 13 C, this trait (rather than days to anthesis) was the first variable chosen by the multilinear model for biomass prediction both when we considered the whole set of genotypes as well as for the larger subsets with a given anthesis time (Table 3). In contrast, NDVI measurements based on UAV-mounted platforms did not predict biomass effectively, although a notable decrease in the NDVI was observed 14 days after drought induction. The results observed in this study not only reveal the importance of environmental conditions on NDVI measurements but also the significance of the developmental stage in which the NDVI is measured when monitoring genotypic adaptation to specific environments. Consequently, the crop developmental stage should be considered in deciding when to measure NDVI because phenology could bias interpretation of the results due to the negative association of biomass with NDVI during advanced developmental stages. Therefore, accessions undergoing NDVI or δ 13 C assessment should have similar phenology, otherwise the phenotypic performance of these traits may be biased. Our results concur with a recent report on wheat [69] that advocates the phenotyping of photosynthetic traits at multiple scales (leaf to canopy) and multiple development stages (vegetative, and pre-and post-anthesis). Additionally, our study highlights the confounding effect of phenology when phenotyping other agronomic traits, a well-known issue when interpreting the effects of QTL for drought-related proxies of biomass and yield productivity [47,70]. Accordingly, any given trait such as NDVI and δ 13 C may have a positive, negative, or null association with the target trait under selection, in accordance with previous reports [38,41,47,55,71]. Therefore, a wide variation in phenology should be avoided in studies that aim to elucidate the relationship between proxy traits and other agronomic traits, otherwise the interpretation of the results may be compromised.
Our study illustrates how a correct phenotyping is only feasible providing (1) a good understanding of phenotyping conditions in terms of traits measured, crop stages, and growing conditions when these traits are evaluated and (2) genotypic variability in intrinsic characteristics, such as phenology, is accounted for.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/11/1679/s1, Table S1: Number of accessions at different flowering times and phenological stages, Table S2: Linear regression of the relationship between the biomass at mid grain filling and the NDVI, Table S3: Linear regression of the relationship between the biomass and δ13C in the flag leaf dry matter sampled before imposition of stress and after two weeks of stress conditions, Figure S1: Image of the field trial at the Maricopa phenotyping platform.