Multispectral airborne imagery in the field reveals genetic determinisms of morphological and transpiration traits of an apple tree hybrid population in response to water deficit

Highlight This research successfully used image-based spectral indices acquired in the field to assess variability of response to drought in a tree mapping population and to detect the related genetic determinisms.


Introduction
According to current climate change models for the 21 st cen tury, an increase in global mean temperatures is expected, with longer or more frequent episodes of extreme temperatures and drought, notably in the Mediterranean basin (IPCC, 2014). Climate change will lead to reconsideration of bree ding programmes for many crops, and optimization of water use by improving the plant water use efficiency and/or the tolerance to drought will become an increasingly important issue (Hamdy et al., 2003;Condon et al., 2004). Although plant behaviour in response to drought can be analysed in terms of survival (McDowell et al., 2008), it more usually refers to the ability of one genotype to yield better than another under more or less severe water deficit. However, while breeding programmes in fruit species have not yet included drought tolerance among the targeted traits, some authors consider that tree response to water scarcity will become a crucial character to consider (Bassett, 2013).
Plants have developed various mechanisms to cope with drought that depend on the duration and intensity of the water deficit, and their responses occur at different temporal and spatial scales, from cell to whole tree level (Jones et al., 2002). One first response to soil drought is stomatal closure, an avoidance mechanism mediated by the hormone abscisic acid (Pantin et al., 2013). A main consequence of stomatal closure is the decrease in CO 2 influx and assimilation, which can lead to carbon depletion. When transpiration is reduced by stomatal closure, the outgoing water vapour flux and the latent heat dissipation are also reduced. Stomatal closure thus induces an increase in leaf temperature with a risk of heat stress (Tardieu, 2005). However, efficiency of stomatal regulation is variable according to species and Tardieu and Simonneau (1998) have shown that plants display contrasting transpiration behaviours (isohydric vs anisohydric) in response to drought. At the intraspecific level, genetic vari ability of stomatal regulation has also been highlighted in apple (Massonnet et al., 2007;Liu et al., 2012) and grapevine (Marguerit et al., 2012;CoupelLedru et al., 2014).
As leaf or canopy temperature can be used as a proxy for stomatal conductance, thermal infrared (TIR) imagery appears as a powerful tool to reveal genetic variability of sto matal behaviour (Jones et al., 2009). Numerous indices have been developed to assess crop water stress from canopy sur face temperature (T s ) with data acquired in signal or imagery mode, from aerial platforms (satellites, aircrafts, unmanned aerial vehicles) or sensors installed directly in the field to observe crop canopies (White et al., 2012). T s minus air tem perature (T a ) is a raw variable that is easy to extract from images, but it is sensitive to rapid changes in environmental conditions such as radiative conditions, wind speed and air vapour pressure deficit (Maes and Steppe, 2012).
The presence of mixed soil/plant pixels is a recurring problem when TIR imagery is applied to phenotyping of heterogeneous covers (Hackl et al., 2012;Jones and Sirault, 2014). It is generally considered that using the vegetation sur face temperature directly is risky, because the weight of mixed or soil image pixels yielded in porous plant covers can create a shift towards the soil surface temperature (Jackson et al., 1981). To overcome the limitations of environmental and soil influence on T s , Moran et al. (1994) developed a Water Deficit Index (WDI) based on the Vegetation Index−Temperature (VIT) trapezoid concept, which is applicable to field crops with varying contributions of bare soil in the aggregated ther mal pixels. This index is particularly suitable for estimation of transpiration rates on heterogeneous vegetation cover. It has been successfully related to the soil moisture and to the plant midday stem water potential (Köksal, 2008;Virlet et al., 2014). Different authors indicated that the use of aerial vec tors (ultralight aircraft or unmanned aerial vehicle) coupled with high resolution sensors enables to distinguish the indi vidual trees within a plant grove, even in the TIR waveband where image resolution is low (Berni et al., 2009;Stagakis et al., 2012). Moreover, the intracrown T s variability has also been used in tree crops as complementary indicator of moderate water stress effect (GonzálezDugo et al., 2012), confirming previous work that considered leaf temperature distribution as better indicator of stress than its average (Fuchs, 1990).
Apart from TIR imagery, plant cover can be characterized by different vegetation indices based on the combination of spectral reflectances measured in visible and nearinfrared (NIR) wavebands (ZarcoTejada et al., 2005). These indices can possibly be acquired by broadband commercial sensors (Lebourgeois et al., 2008). In the remotelysensed image, reflectance in the Red band is affected by light absorption of leaf pigments (mainly chlorophyll a), while the NIR waveband is affected by the scattering in the medium (ZarcoTejada et al., 2005). Therefore, vegetation indices computed from Red and NIR, such as the normalized difference vegetation index (NDVI), can be related to canopy structure and biomass production (ZarcoTejada et al., 2005) and also considered as indicators of tree vigour. However, NDVI is sensitive to low chlorophyll concentration (Peng and Gitelson, 2011) and it also tends to saturation when leaf area index (LAI) is higher than 3 or 4. Two other indices only retrieved from visible bands were used: the visible atmospherically resistant index, VARI, which shows a better sensitivity to higher values of vegetation cover fraction (Gitelson et al., 2002) and the simple ratio pigment index, SRPI, which enables characteri zation of the crop nitrogen status, being sensitive to change in the pigment relative content (chlorophyll vs carotenoids) (Peñuelas et al., 1994(Peñuelas et al., , 1995. Recent studies on field crops, e.g. wheat (Babar et al., 2006;Comar et al., 2012), maize (Cairns et al., 2012) and cotton (AndradeSanchez et al., 2014) assessed potentiality of vegetation indices to be used for largescale phenotyping. More generally, plant phenotyping based on multispectral or hyperspectral imagery shows promise as a noninvasive method adapted for screening a wide range of individuals in a short period of time. Connecting genotype to phenotype on large datasets currently sustains the development of pheno mics (Furbank and Tester, 2011;Fiorani and Schurr, 2013).
To date, quantitative genetic analyses of tree features in fruit crops have mostly concerned disease resistance, yield and production regularity (Guitton et al., 2012;Celton et al., 2014), and plant architecture (Segura et al., 2008). Owing to lowthroughput techniques, few studies on genetic determi ni sms of traits related to water use have been undertaken in these crops except recently in grapevine (Marguerit et al., 2012). Other perennials like forest trees have been compared in natural environments (Brendel et al., 2008) and controlled environments (e.g. Salix: RönnbergWästljung et al., 2005;Populus: Street et al., 2006) to distinguish wellirrigated and Field multispectral imagery and QTL analysis of apple tree response to water deficit | 5455 water deficit conditions and to study the genetic/genomic bases of responses to drought and/or water use efficiency.
In this study, we assumed that a genetic analysis could be performed on an apple segregating population submitted to contrasting water regimes, considering different traits mainly issued from airborne multispectral imagery. An experiment was conducted in two successive growing seasons, during which imagebased phenotypic variables and agronomic traits such as fruit production or trunk diameter (a proxy for tree vigour) were analysed for both wellwatered and waterstress conditions, as well as the difference between the two for a given genotype. For each imagebased variable, we considered the mean value of a representative tree crown zone and the variation within this zone, on which mean broadsense heri tability was computed from genetic linear models. Using a genetic map, quantitative trait loci were detected. Altogether these results demonstrate the relevance of airborne imagery for highthroughput phenotyping of individual trees in the field for their response to water stress and provide the first demonstration that QTL detection could result from such methodology and plant descriptors.

Field experiments and meteorological measurements
The apple tree population studied consisted of progeny derived from a 'Starkrimson'×'Granny Smith' cross, characterized by variability in tree vigour, architectural traits (Segura et al., 2008), biennial bear ing (Guitton et al., 2012), hydraulic traits (Lauri et al., 2011) and stomatal regulation in response to vapour pressure deficit (Regnard et al., 2009). In February 2007, four replicates of 122 F1hybrids and their two parents were grafted onto M9 rootstock and randomly planted in an experimental field at the INRA Melgueil experimental station (Diaphen platform, southeast of France, 43°36ʹ N, 03°58ʹ E). Plantation consisted of 10 rows oriented northwest-southeast, with 5 × 2 m planting distances. The orchard management was per formed consistently with professional practices, throughout the study. Automated soil resistivity profiling conducted in March 2009 showed that the soil of the trial plot (at depths of 0−50 cm and 50−100 cm) could be considered spatially homogeneous for water holding capacity, and this was confirmed by soil profile descriptions. The field plot was irrigated using a system of microsprayers located in the rows, with one emitter per tree. During summer, contrasting hydric regimes were established. Full irrigation was ensured in half of trees [two replicates per genotype, wellwatered trees (WW)], while irrigation was withheld in the other half, resulting in pro gressive summer soil drought [two replicates per genotype, water stressed trees (WS)] since the summer rainfall was negligible. Trees submitted to water deficit during summer were the same during the 2010 and 2011 seasons, and three dates per year were studied, repre senting various water regimes in order to disentangle genotypic and environmental effects in the tree response. Water regimes developed in WW and WS treatments are illustrated by the soil hydric potential mean values (Ψ soil , Table 1A). Micrometeorological data acquired at field included global radiation (R g ), air temperature (T a ), air relative humidity (HR), air vapour pressure deficit (VPD), wind speed (u) and rainfall (Table 1A).

Image acquisitions
The image acquisition system from the ultralight aircraft consisted of two commercial digital cameras (either Canon EOS 400D or 500D, with 10.1 and 15.1 Megapixel CMOS sensors, respectively, Table 1B) equipped with 35mm lenses, and one FLIR B20HSV (FLIR Systems Inc., Wilsonville, USA) thermal infrared camera (320*240 matrix) (for details, see: Lebourgeois et al., 2008Lebourgeois et al., , 2012Virlet et al., 2014). One camera acquired visible images in red, green and blue bands (RGB). The second was modified according to Lebourgeois et al. (2008Lebourgeois et al. ( , 2012 to obtain images in nearinfrared (NIR). Three flights per year were performed during the summers of 2010 and 2011 (Table 1A, 1B). In 2010, flights were realized for low, intermediate and severe water constraints, respectively 8, 27 and 41d after the beginning of drought (Dates 1, 2 and 3). In 2011, the first flight (Date 4) occurred 17 d before the beginning of the drought period, before WW and WS differentiation, while the second and third flights (Dates 5 and 6) were performed respectively 14 and 34d after the beginning of the drought treatment. During the period of water deprivation (i.e. at Dates 1, 2, 3, 5 and 6) WS trees were not irrigated.

Spectral image preprocessing and indices computation
Image preprocessing was performed with Erdas Imagine 9.3 software (Intergraph Corporation, Huntsville, USA). Procedure of ortho rectification for RGB and NIR images and radiometric normaliza tion on invariant field targets between dates are fully described in Lebourgeois et al. (2008Lebourgeois et al. ( , 2012 and Virlet et al. (2014), as well as image geolocation. Thermal infrared images issued from the six acquisition dates were orthorectified on the base of both RGB and NIR images and geolocated as well. For each of the six dates, the dif ference between the surface and air temperature (hereafter referred to as TsTa) was obtained by subtracting from each pixel value of the TIR images the air temperature acquired at ground level. Spatial resolution of RGB and NIR images was lowered from initial reso lution (c. 3−5 cm) to that of TIR image (30 cm). From RGB and NIR bands, three vegetation indices were computed: NDVI, VARI and SRPI (Table 2). NDVI and TsTa were combined to compute the water deficit index (WDI) as described in Virlet et al. (2014).
For each tree, multispectralbased index values were extracted from a 60 cm radius buffer zone containing the central upper part of the tree crown. From each buffer (12−16 pixels), mean and standard deviation, SD, were retrieved and considered as two complementary variables characterizing the vegetation response of individuals. As SD characterized the variation occurring inside the buffer zone, it indicated the degree of heterogeneity of the crown structure for the vegetation index and the variability of transpiration rates for the stress indices.
In planta measurements Trunk circumference (TC) of each tree was tapemeasured 15 cm above the grafting point each year in February. On that basis, the trunk crosssectional area (TCSA) considered as representative of tree vigour was calculated (Table 2). Fruits were harvested each year between 22 August and 2 September before the resumption of irriga tion, irrespective of the real maturity picking date (September). The number of fruits per tree and the harvest fresh mass (kg per tree) were determined.

Data analyses
Statistical analyses were performed using R software v.2.13.2. (R Development Core Team, 2011). For each variable, phenotypic means were computed from each tree for (i) WW and WS conditions confounded, (ii) WW and WS considered separately and (iii) the dif ference between them, hereafter referred to as the differential index (DI). For example, the phenotypic mean value of NDVI indepen dent of water treatment is referred to as NDVI, while sdNDVI_WS refers to the phenotypic mean value of the SD in the WS condition.
For each variable, two mixed linear models were built. The first one was used to analyse the response of variables in both WW and WS trees. The second one was used to analyse the drought response of each genotype, through the DI obtained. The first mixed linear models included the irrigation modality (M), date (D) or year (Y), which were considered as fixed effects, while the genotype (G), the interactions between genotype and irrigation modality (G×M), and genotype and date or year (G×D or G×Y) were considered as random effects. For each variable, a selection of the best model was performed through minimization of the Bayesian Information  Peñuelas et al., 1994, 1995Lebourgeois et al., 2012   Field multispectral imagery and QTL analysis of apple tree response to water deficit | 5457 Criterion (BIC). For the DI, effects considered in the mixed linear model were the same as mentioned above, but without M and G×M. For each trait, best linear unbiased predictors (Blups) were extracted for estimation of the G effect, which was considered as independent from the irrigation modality and the date (or year) of experimenta tion and is hereafter referred to as GBlup. The Blups corresponding to G×M and G×D effects were computed for each irrigation regime (WW or WSBlup), and date (or year) considered separately. For each variable, when G and interaction effects were included, broadsense heritability of the mean (h 2 b ) was estimated as follows: where n is the number of trees per genotype (two in the present case), and n M the number of irrigation modalities (two in the pre sent case). When a G×D (or G×Y) effect was included in the model, the denominator also integrated G×D (or G×Y) variances and was divided by the number of dates (six) or years (two). The residual variance σ ε 2 was divided by the product of the number of trees per genotype and per irrigation modality and the number of dates (or years). This led us to estimate the broadsense heritability value of the mean of phenotypic values which accounts for the repetitions of each genotype that were present in the experimental plot, according to Gallais (1989). Phenotypic variables were considered heritable if h 2 b values were greater than 0.2.

QTL mapping
The QTL analysis was performed using means and Blups extracted per genotype (Gmeans, GBlups) for each variable. A consensus genetic map of STK and GS, which integrated 177 SSR and SNP genetic markers, was used for QTL mapping (Guitton et al., 2012). QTL analyses were carried out using MapQTL ® 6.0 (Van Ooijen, 2009). First, a permutation test was performed to determine the loga rithm of the odds (LOD) threshold at which a QTL was declared significant, using a genomewide error rate of 0.05 with 1000 permu tations of the data (Van Ooijen, 2009). In a second step, an interval mapping analysis was carried out with a step size of 1 cM, with a LOD score higher than the threshold. Finally, the nearest marker to each QTL peak was selected as a cofactor to perform a multiple QTL mapping (MQM) (Van Ooijen, 2009). Each significant QTL was characterized by its LOD score, its percentage of explained phe notypic variation, and its confidence interval (in cM) corresponding to a LOD score drop of 1 or 2 on either side of the likelihood peak. QTLs that showed clearly overlapping confidence intervals, close LOD peaks and similar allelic effects, were considered to colocalize. When a QTL was detected with at least two cofactors, models considering markers and their interactions as cofactors were cons tructed using a backward procedure under R software v2.13.2. Models were selected based on the minimum Akaike Information Criterion values (AIC). In the selected model, the global percentage of phenotypic variation (global R 2 ) was then estimated. When one marker was derived from only one of the parents, the nearest maker included in the QTL and deriving from both parents was chosen. The location of QTLs on the genetic was finally illustrated using MapChart® (Voorrips, 2001).

Variance analysis and heritability
Models selected for vegetation indices were similar whether means or SDs were considered (Table 3A). All vegetation indices were significantly impacted by G, D and M effects.
For NDVI, the model included only the G×M interaction, whereas for VARI, SRPI, TsTa and WDI variables G×M and G×D interactions were also taken into account. Concerning tree vigour and fruit production, the models selected included G and Y effects. For TCSA only the G×M interaction was retained in the mixed linear model, while the G×Y interac tion was retained for fruit number (NbFr) and both G×M and G×Y interactions were retained for fruit yield biomass (BmFr). For all DI variables (Table 3B) the models selected included G and D (or Y) effects. G×D was also included for sdVARI_DI, sdSRPI_DI, sdTsTa_DI while G×Y was included for NbFr_DI only. It is noticeable that the random interaction effects (G×M and/or G×D) were generally lower than the G effects.
Broadsense heritability h 2 b for both WW and WS (Table 3A) showed medium to high values (0.49 to 0.77) except for sdWDI and BmFr, whose heritability was low (0.31). For DI variables (Table 3B), fairly high h 2 b values (0.50 to 0.70) were found for NDVI, sdVARI, SRPI, sdSRPI, TsTa, sdTsTa and TCSA. In contrary, h 2 b for the other variables, including WDI_DI, was much lower (0.17 to 0.38) than that found for WW and WS. Moreover, higher h 2 b were found in VARI_DI, SRPI_DI, and TsTa_DI for SD values than for mean values.

Correlations between variables
High pairwise positive correlations were observed between NDVI, VARI and SRPI for the Gmean, WW, WS and for DI (Pearson's r from 0.45 to 0.88, Table 4) even though lower r values were found between VARI and SRPI for Gmean and WW (Tables 4A, 4B). These three variables were significantly and negatively correlated with sdNDVI and sdSRPI for the Gmean, WW, WS and for DI (from −0.46 to −0.66, Table 4), despite much lower correlation being found between sdNDVI and these variables for DI (from −0.06 to −0.32, Table 4D). Moreover, sdNDVI, sdVARI and sdSRPI presented pairwise positive correlations for the Gmean and WW (0.33 to 0.88, Table 4A, 4B). For WS and DI, only sdSRPI was significantly correlated with sdNDVI and sdVARI (0.74 and 0.42 respec tively, Table 4C, 4D). WDI was highly and positively corre lated with sdWDI (0.49 to 0.52). The trunk diameter variable, TCSA, presented generally moderate to high positive correla tions with NDVI, VARI and SRPI. The highest correlation was observed with NDVI, either for the Gmean, WW, WS or for DI (0.55 to 0.67). Variables relative to fruit production, NbFr and BmFr, were highly intercorrelated (0.75 to 0.85), and a high positive correlation of these variables with TCSA was also observed, particularly for BmFr (from 0.50 to 0.55). Moreover, BmFr was positively and more highly correlated to NDVI, VARI and SRPI (0.23 to 0.52) than NbFr (0.13 to 0.32) for the Gmean, WW and WS.

QTL detection
Seventyfour QTLs were detected, mapping over 16 of the 17 linkage groups (LGs) of the consensus STK×GS genetic map. As 56 of these QTLs were found only at specific dates, they are not detailed in the following text. The complete list of QTLs detected is presented in Supplementary Table S1 and Supplementary Fig. S1. The results exposed hereafter ( Fig. 1; Table 5) are focusing on the 18 most reliable QTLs that were mapped over nine LGs. These QTLs were detected for GBlup or for the Gmean, and in some cases for both, whatever the date.

QTLs for traits related to vegetation indices NDVI, VARI and SRPI.
For GBlup of NDVI and SRPI variables, three QTLs were detected independently of the environment (water regime or date): two QTLs concerned sdNDVI and sdSRPI and pre sented a common zone located on LG08. They explained 14.2% and 15.8% of the variability respectively, and were both characterized by female allelic effects. One QTL for SRPI was detected on LG09 and explained 14.3% of the variability. It showed female, male and dominance effects.
Four QTLs were detected for specific Gmeans of WW on four different LGs. Two of these QTLs, related to NDVI_ WW and SRPI_WW, were detected on LG14 and LG09, respectively. They explained 14.9% and 15.7% of the varia bility, and both of them resulted both from female, male and dominance effects. Two other QTLs, related to sdVARI_WW and sdSRPI_WW, were mapped on LG01 and LG08, and resulted from dominance and female effects, respectively. The QTL for sdVARI_WW was also identified for GBlup.
For WS, one QTL was identified for sdNDVI_WS at the top of LG06. It explained 14.1% of the variability and mainly resulted from female effects. For DI, one QTL was detected for VARI_DI and mapped at the bottom of LG03 for both GBlup and Gmean. It explained 14.5% of the variability and resulted from female, male and dominance effects.

QTLs for traits related to tree foliage transpiration
One QTL was detected for WDI in WW condition, at the top of LG03 for both GBlup and Gmean. It explained 14.7% of the variability and was characterized mainly by female and male effects. Finally, two QTLs were detected for TsTa_2_DI (at Date 2) on LG05 and LG06, respectively. The global linear model indicated an interaction between these two QTLs, which together explained 20.8% of the variability. They Table 3. Description of fixed (M, modality; D, date; Y, year) and random (G, genotype) effects used in selected mixed linear models For each variable, models related to phenotypic values in (A) WW and WS, and (B) models related to DI (differential index: difference of the variable between WS and WW trees) were built. Percentage variances of each random effect and of the residuals (Res), and broad-sense heritability values (h 2 b ) are indicated.   LG14 Fig. 1. Genomic positions of the QTLs detected on the consensus 'Starkrimson'×'Granny Smith' (STK×GS) map. QTLs are represented by boxes, in which length represents the LOD-1 confidence interval and extended lines represent the LOD-2 confidence interval. Boxes relative to QTLs for mean values of variables are in white, and those relative to QTLs for standard deviations SD are in black. QTL detected for G-Blups are in bold type and * stand for QTLs detected for G-Blups and G-means. Field multispectral imagery and QTL analysis of apple tree response to water deficit | 5461 mainly resulted from female and male effects. The QTL which mapped on LG06 colocalized with the QTL for sdNDVI_WS.

QTLs for traits related to tree vigour and fruit production
Two QTLs were detected in relation to the tree vigour (TCSA), the first one for TCSA_WS, and the second one for TCSA_DI. Both mapped on LG10 and colocalized near the MdVRN1b_S marker. These QTLs explained 14.5% and 15.4% of the variability, respectively, and mostly resulted from female allelic effects. For fruit number, one QTL was detected for WW (NbFr_WW) at the bottom of LG13. It explained 19.6% of the variability. Another QTL was detected for NbFr_DI at the same position, explaining 17.8% of varia bility. These two QTLs mainly resulted from male and female allelic effects. Finally, two QTLs were found at the top of LG05 for NbFr_WS and BmFr_WS. They explained 15.6% and 15.9% of the variability, respectively, and both resulted from female and dominance effects. Interestingly, these two QTLs colocalized with the QTL identified for TsTa_2_DI.

Variability of the phenotypic traits
Due to, on the one hand, the changes in environmental condi tions and imagery flight parameters and, on the other, the dif ficulty to apply comparable water constraints from one year to the next, the analysis of the population behaviour undertaken through linear mixed models did take into account the large variability of phenotypic traits (Supplementary Table S2). Among vegetation indices, NDVI appeared the most stable vegetation index, independent of the environment and acqui sition conditions (no G×D interaction) whereas other spec tral indices chosen proved to be sensitive to drought, with a G×D interaction revealed for VARI, SRPI, TsTa and WDI. By contrast, the absence of G×D interaction for Gmeans for the DI variables suggested that the use of the difference between waterstressed and wellirrigated trees somewhat smoothed out the interdate variations.
To our knowledge, this study is the first one to make use of spectral indices to assess genetic variability in perennial plants in response to drought, and to analyse the related determinisms. As a consequence, comparisons in the ensuing discussion are often referring to results obtained in annual crops. In the present study, the vegetation indices used were extracted from a buffer zone of the same size located in the central zone of each tree crown. As such, the value of these vegetation indices must be considered as more related to the vegetation cover fraction and biomass production rather than to the foliage physiological properties. Whatever the traits considered, either relative to vegetation or to transpi ration, moderate to high values of broadsense heritability were found, indicating an interesting contribution of multi spectral imagery for genetic analysis of these traits in a tree population. Concerning the vegetation indices, our results were consistent with those found in annual crops, where high  NDVI, VARI, SRPI, T s T a , WDI, TCSA, NbFr and BmFr in well-watered (WW) and/or water-stress (WS) conditions and for the differential index DI (WS−WW) QTLs detected for G-Blups are in bold type and * stand for QTLs detected for G-Blups and G-means.

Traits
LG heritability values were also found (e.g. in wheat: Babar et al., 2006;in cotton: AndradeSanchez et al., 2014). Heritability values found for agronomic traits in apple, TCSA, NbFr and BmFr, were in the same order of magnitude as those for spec tral indices. However, TCSA, which is an integrative variable, exhibited the highest h 2 b value whereas the lower values for NbFr and BmFr likely resulted from the influence of Y, and Y and M effects, respectively. When heritability was computed for DIs, values of most variables stayed high except for WDI_ DI, probably because of the composite nature and complexity of this trait. Concerning TsTa and TsTa_DI, high heritability values were found, consistent with previous results found in wheat (Mason et al., 2013;Rebetzke et al., 2013).

Trait correlations, QTL detection and co-localization
Positive correlations between the vegetation indices, TCSA and fruit production were observed. In particular NDVI exhibited the highest correlation with the trunk diameter, TCSA, and with fruit yield biomass, BmFr. As well as NDVI, TCSA is generally related to plant size, leaf area and light interception (GonzálezTalice et al., 2012), which suggests that NDVI is a good indicator of vigour and development in apple trees. However, QTLs that were found for vegetation indices, fruit yield and TCSA did not colocalize, and this indi cated that genetic determinisms controlling these traits likely differ. Moreover, while the three vegetation indices presented high and positive intercorrelations, related QTLs did not colocalize. Indeed, QTLs for NDVI and SRPI were mostly detected on two different LGs: LG14 and LG09, respec tively. This could be explained by the spectral bands used: the former vegetation index, making use of NIR for compu tation, is more related to canopy structure than the second one, only computed from visible bands that are more related to light absorption (ZarcoTejada et al., 2005;Lebourgeois et al., 2012). Otherwise, the QTL found for VARI_DI, which highlights differences between WS and WW tree responses, suggests that drought could affect the radiation absorption capacity, by lowering the fractional vegetation cover.
Intracrown variations of vegetation indices, sdNDVI and sdSRPI, presented high and negative correlations with means of NDVI and SRPI, and moderate and negative correlations with TCSA. This indicates that the highest density of vegeta tion was less variable than the lowest, but also that the vegeta tion heterogeneity was lower where tree vigour increased. As the intracrown variability could be an indicator of branching patterns or leaf clumpiness (Da , this needs further investigation. Similarly, as WDI and sdWDI were posi tively correlated, this could be due to spatial heterogene ity in stomatal conductance within the tree crown in response to moderate stress (GonzálezDugo et al., 2012). The QTL detected for sdNDVI_WS could be attributed to the variation of leaf rolling over genotypes in response to drought, a phe nomenon also observed in other species and limiting plant cover fraction (e.g. maize: Lu et al., 2012). Furthermore, sdNDVI and sdSRPI were strongly intercorrelated and QTLs for these variables colocalized at the middle of the LG08 (Table 5; Supplementary Table S2). The location of these QTLs also matched with a QTL zone for traits involved in gas exchange, xylem conductance and fruit production on the STK×GS population (Regnard et al., 2009;Lauri et al., 2011, Guitton et al., 2012. These colocalizations could be explained by an increased capacity of the plant to transport water, carbohydrates and sugar to the growing organs, as sug gested by Lauri et al. (2011). Nonetheless, these colocaliza tions might also be explained by a pleiotropic effect of these QTLs, or by clustering of functionally related genes (Cai and Morishima, 2002). Gene clusters have already been reported in apple for various traits such as resistance to pathogens (Xu and Korban 2002;Baldi et al., 2004). However, discrimina ting between linked and pleiotropic QTLs was not practicable in the present study, considering the limited population size and the density of the genetic map available.
Among QTLs detected for fruit production, two of them, NbFr_WW and NbFr_DI, were located at the bottom of LG13. This zone was adjacent to the one found for bien nial bearing on the same STK×GS population (Guitton et al., 2012). Otherwise, a yearspecific QTL, NbFr_11, was detected at the same location as NbFr_WW and NbFr_DI ( Supplementary Fig. S2), which could confirm the importance of this zone in the control of biennial bearing. In addition, QTLs for fruit production in WS trees (NbFr_WS and BmFr_ WS) and for leaf temperature (TsTa_2_DI) colocalized on LG05. Although matching for those traits was only tempo rary, i.e. when the difference between WW and WS treatments was considered at Date 2 for TsTa, it illustrates the negative relationship between yield and leaf temperature. Indeed, as stated by Naor and Girona (2012), a positive link between plant yield and evapotranspiration is generally observed, and the increase in leaf temperature (here in WS trees) is an indi cator of lower transpiration rate and likely of lower carbon assimilation. Reduction of stomatal conductance in WS trees could likely be invoked for limitation of fruit production in this case because water constraint was severe (no irrigation occurred during the summer period). Such a causal relation between water withholding and its effect on yield reduction is nevertheless not straightforward in fruit trees (Bassett, 2013), particularly when a moderate water deficit occurs during stages of low fruit growth (Goodwin and Boland, 2002).
Although one QTL was detected on LG03 for WDI_WW, no colocalization with fruit production variables was observed. Moreover, genetic correlations of WDI with fruit production variables during the two years of study (Table 4) were not significant, whereas significant negative correlations of this variable were observed with NbFr and BmFr (r value of −0.53 and −0.55 respectively, data not shown) when only year 2011 was considered. This negative correlation in a specific year could be attributed to the propensity to biennial fruit bearing already shown on this progeny (Guitton et al., 2012).
It has been recently shown that wild germplasm of Malus is exhibiting a certain range of tolerance to drought (Bassett, 2013), but available information on the commercial paren tal genotypes used in the present study is scarce. González Talice et al. (2012) suggested that the 'Granny Smith' cultivar has a lower hydraulic conductance and/or stomatal conduc tance than that of the 'Gala' cultivar, indicating a stronger sensitivity of the former to high evaporative demand and/or soil drought. In contrast, the 'Starkrimson' cultivar is well known for its photosynthetic efficiency (Yang and Wang, 1993), while its response to drought is not documented. In the present study, prevailing female allelic effects were found for almost all variables, suggesting a larger polymorphism and allelic contrast in the 'Starkrimson' parent than in 'Granny Smith'. 'Starkrimson' could thus provide interesting alleles for adaptation to drought, even though other Malus genetic backgrounds need to be explored in the next future.
Partial homology between LGs that has been described in the 'Golden Delicious' apple genome (Velasco et al., 2010) led us to examine homologous regions in which main QTL zones were detected. The median zone on LG08, around CH02g09 in which many QTLs were detected, matches the top of LG15, above CH03b6, where a QTL was detected for SRPI_1_DI. Similarly, the QTL zone on the top of LG03 detected for NDVI_4_WW corresponds to the region on LG11 on which a QTL was detected for sdNDVI_1_DI. By contrast, the QTL zones found in LG05 and LG10 for NbFr_WS and NbFr_11_DI, respectively, were located on chromosomal fragments that are inverted on their respective chromosome and therefore did not matched. Similarly, the QTLs detected on LG06 that were located either above or below CH03d07 were compared to those detected on LG14, without finding evident homology in those regions. These findings suggest that further investigations of QTL zones homologies will be required, along with identification of candidate genes in the zones of highest interest.
To summarize, this work is an important step in the study of tree field phenotyping for response to abiotic stress. It confirmed-if proof was needed-the strong potential of remote sensing tools as a method for screening a large panel of genotypes. Airborne imagery proved relevant to acquire simultaneous information on a tree population, notably for characterizing transpiration behaviour at the individual tree scale as a result of images yielded in the thermal infrared domain. Indices derived from highresolution airborne field imagery appeared to be highly heritable and enabled detec tion of a large number of QTLs, for vegetation and water stress indices, and the tree response to water deficit. This study opens future avenues for analysis of candidate genes related to foliage response to drought, and may contribute to future selection of new plant woody plant material bred for its response to drought and/or water use efficiency.
Supplementary Table S1. Total list of QTLs detected for all phenotypic variables in wellwatered and/or waterstress conditions and for the differential index.
Supplementary Table S2. Values of traits (mean and SD) investigated in wellwatered or waterstressed conditions and for the differential index.