Image Phenotyping of Spring Barley (Hordeum vulgare L.) RIL Population Under Drought: Selection of Traits and Biological Interpretation

Image-based phenotyping is a non-invasive method that permits the dynamic evaluation of plant features during growth, which is especially important for understanding plant adaptation and temporal dynamics of responses to environmental cues such as water deficit or drought. The aim of the present study was to use high-throughput imaging in order to assess the variation and dynamics of growth and development during drought in a spring barley population and to investigate associations between traits measured in time and yield-related traits measured after harvesting. Plant material covered recombinant inbred line population derived from a cross between European and Syrian cultivars. After placing the plants on the platform (28th day after sowing), drought stress was applied for 2 weeks. Top and side cameras were used to capture images daily that covered the visible range of the light spectrum, fluorescence signals, and the near infrared spectrum. The image processing provided 376 traits that were subjected to analysis. After 32 days of image phenotyping, the plants were cultivated in the greenhouse under optimal watering conditions until ripening, when several architecture and yield-related traits were measured. The applied data analysis approach, based on the clustering of image-derived traits into groups according to time profiles of statistical and genetic parameters, permitted to select traits representative for inference from the experiment. In particular, drought effects for 27 traits related to convex hull geometry, texture, proportion of brown pixels and chlorophyll intensity were found to be highly correlated with drought effects for spike traits and thousand grain weight.


INTRODUCTION
Barley (Hordeum vulgare ssp. vulgare L.) is one of the most important crops worldwide because of its multipurpose usage in human diet and as animal feed. Although it is known to adapt to a wide range of environments, in Europe new cultivars have been bred under favorable conditions, which led to the narrowing of genetic diversity in agronomical properties, including resistance to environmental stresses like shortage of water. The gene pool for tolerance to water scarcity in modern elite cultivars is very limited. Broadening of the genetic diversity provides the basis for plant improvement. Hence, success in breeding new varieties with improved tolerance to water shortage or heat can be achieved through the use of wild relatives, landraces or varieties growing in dry areas as donors of the resistance (Ceccarelli, 1994;Ceccarelli and Grando, 1999;Grando et al., 2001;Cattivelli et al., 2008).
Understanding adaptation of plants to their environment is a key issue that is addressed in many fields of study and may also contribute to breeding crop plants adapted to sub-optimal conditions. Water scarcity (also referred to as drought stress) is a condition associated with insufficient soil moisture available to provide satisfactory crop production. Development of drought tolerant cultivars becomes increasingly important in changing climate. Many different approaches have been used to study the nature of plant reactions to defined levels of drought stress, including physiological processes investigations (Buschmann et al., 2000;Sanchez et al., 2002;Jones et al., 2003) and phenotyping on conveyor systems in glasshouses with controlled irrigation systems (Tuberosa, 2012;Honsdorf et al., 2014).
Accurate quantification of traits observed in plants grown under conditions of limited watering is crucial for identifying loci of interest in the genome. While in construction of genetic maps high-throughput genotyping platforms are used routinely, plant phenotypes are still assessed mainly by conventional procedures, which are time-consuming, labor-intensive, low-throughput, and usually destructive. Studying the drought stress response is particularly challenging as its impact on plant performance is a dynamic process that occurs across space and time. Endpoint measurements are insufficient to asses and analyze dynamic responses. Thus, phenotyping has become the major operational bottleneck limiting the power of genetic analysis (Hartmann et al., 2011;Cabrera-Bosquet et al., 2012).
In recent years, automation, imaging, and software solutions have paved the way for many high-throughput phenotyping studies (Munns et al., 2010;Busemeyer et al., 2013;Chen et al., 2014;Honsdorf et al., 2014;Paulus et al., 2014). Regardless of the platform setting, the goal is to evaluate phenotypic properties of plants using automated processes in a non-invasive way. Automated systems have been successfully applied to investigate numerous components of plant growth, and may be used to help tackle basic research questions when combined with genetic strategies (Famoso et al., 2010;Muraya et al., 2017;Neumann et al., 2017).
The aim of the present study was to screen changes in growth, architecture and physiology traits measured by high-throughput imaging in spring barley RILs population under drought and control (well-watered) conditions, to identify traits relevant for the description of drought response and to assess an association between the traits measured during the time-course experiment and yield-related traits measured after harvest. The purpose was also to propose a method of selection of traits produced by image analysis in order to obtain a non-redundant feature set characterizing barley drought response. We assume that quantitative analyses of plant structure traits during development under water-limited environments will permit the identification of lines with enhanced resilience to water deficit.

Plant Material
Spring barley (Hordeum vulgare L.) recombinant inbred lines (RILs) population (hereafter referred to as MCam) derived from a cross between European and Syrian spring barley cultivars -Maresi and Cam/B1/CI08887//CI05761 -was used in our studies (Mikołajczak et al., 2016). Maresi is a German semidwarf cultivar with the pedigree Cebeco-6801/GB-1605//HA-46459-68, and Cam/B1/CI08887//CI05761 (hereafter referred to as Cam/B1/CI) is a Syrian breeding line adapted to dry environments. The Syrian genotype was supplied to Dr A. Górny (Institute of Plant Genetics PAS, Poznań) by Drs S. Grando and S. Ceccarelli from ICARDA in Aleppo, and European cultivar was obtained from the collection of IPG PAS Poznań. RILs were derived by the single-seed descent technique (Goulden, 1939) until F 8 generation. 95 RILs and the parental genotypes Maresi and Cam/B1/CI were examined in the experiment.

Plant Cultivation and Phenotyping
After 4 weeks of pre-cultivation, plants were transplanted onto a high-throughput platform. They grew under controlled greenhouse conditions and were phenotyped on a daily basis using the fully automated system consisting of conveyor belts, a weighing and watering station, and three imaging sensors. Each genotype was represented by 16 plants located in four carriers, each with four pots. Genotypes were randomized throughout the plant growth area. Two environmental conditions were simulated: control -optimal watering (for two carriers), and stress -limited watering (for another two carriers). The growth conditions in the greenhouse were set to 20 • C during the day and 16 • C at night with relative humidity >65%. The daylight period lasted 16 h starting at 6 AM. Using automated, targetweight based watering, control plants remained well-watered at a field capacity of 70%, and those in stress conditions were kept at 20% field capacity.
Drought stress was applied from the fifth day after placing the plants on the platform (35th day after sowing, DAS 35) until day 18th (DAS 48). After the stress period, plants were rewatered to optimal field capacity and kept well-watered again for another 2 weeks. Fertilization was carried out twice with Combo Hakaphos blau (200 cm 3 absolute volume per plant). In order to avoid any position effects, the carriers were shuffled daily for one lane, every 3rd day within each lane for 11 positions. Imaging started 4 days after transplanting onto the platform (DAS 31) and was performed daily for 33 days until DAS 63. Observations were not done on DAS 56, but, for simplicity, the data were analyzed as a continuous series of 32 time points. One top and three side view images were taken covering the visible range of the light spectrum (VIS), static fluorescence signals (FLUOR), and the near infrared spectrum (NIR) (Junker et al., 2015). In this manner, 190,493 images were acquired for all genotypes and two treatments during the whole image phenotyping period. The image processing pipeline (Integrated Analysis Platform, IAP, Klukas et al., 2014) provided 376 traits that were subjected to analysis (Supplementary Table 1).
After imaging in the phenotyping platform, the plants were moved from the platform to a non-automated greenhouse, where they were cultivated for another 40 days under optimal watering conditions until ripening. After harvesting, 12 plant architectureand yield-related traits were measured for each genotype on two plants grown in control conditions and two plants grown in drought ( Table 1).

Statistical Analysis
Data (per-plant observations) obtained from the image analysis pipeline were first submitted to outlier removal by application of the Grubbs test (Grubbs, 1950). Then, data were averaged over plants within carriers giving as input for statistical analysis, for each trait, a table for 97 genotypes × 2 conditions × 2 biological replications (carriers) × 32 days. Pearson correlation coefficients between time courses of observations done in replicates (carriers of the same genotype) were computed and only the traits for which both the correlation between two carriers of non-stressed plants and the correlation between two carriers of stressed plants were both >0.2 were selected for further analysis. A linear mixed model with fixed effects of treatment and random effects of genotype × treatment interaction (with unstructured covariance matrix for treatment term) was fitted for data obtained on each day, with REML estimation of variance components, estimation of genetic correlation between conditions, and Wald test of mean drought effect. Relative drought effects (RDE) for genotypes were computed as 100 * (mean value under drought -mean value under control)/mean value under control. Grouping of image-based traits with respect to time profiles of computed parameters (drought effects, genetic correlations, correlations with after-harvest traits) was performed by searching for the optimum clustering into 1, 2, . . . , 10 groups using the minimum withingroup sum of squares (WSS) criterion and selecting the proper number of groups by inspection of the scree plot of WSS values.
Phenotypic data obtained after harvesting were analyzed by a mixed linear model analogous to the one used for imagederived data. Pearson correlation coefficients were computed between RDE for each after-harvest trait (computed as above) and each image trait on each day (using the sets of RDE for all genotypes). Prediction of RDE for after-harvest traits by RDE for image-based traits was done by fitting a partial least squares (PLS) regression model with two-groups crossvalidation.
Data processing until obtaining of mean values over carriers was done in R 3.3.1 (R Core Team, 2018). Subsequent statistical

Analysis of Post-harvest Traits
For all 12 traits observed after harvesting, the distributions of mean values for lines under drought was shifted to lower values (Supplementary Figure 1A), which resulted in significant mean drought effects (Table 1; P < 0.001). The variance components for RILs were the largest for traits describing main spikes (LSm, NSm, NGm, GWm) and larger under control than under stress conditions for all traits except TGW ( Table 1; see also differences in range of distributions shown in Supplementary Figure 1A). The genetic correlation between conditions was the largest for length of spikes and for number of spikelets (LSm, LSl, NSm, NSl), and the smallest for GWl and TGW. Phenotypic correlation between conditions was very low for number of tillers (NTP, NTT); genetic correlation could not be estimated for these traits due to the presence of negative estimates of variance components. The biplot in Supplementary Figure 1B shows a smaller role of traits describing the numbers of tillers than of other traits in discriminating non-stressed and stressed plants. Figure 1C) revealed high positive correlations between both traits describing number of tillers (NTP, NTT) and between traits describing spike properties, with correlations larger in control conditions than under stress. Correlation coefficients revealed a weak association between plant height and other traits both in non-stress and stress conditions, although in drought, this association appeared to be slightly stronger. Total number of tillers in non-stress conditions was negatively correlated with spike traits, i.e., length, number of spikelets and grains in main and lateral spikes, as well as thousand kernel weight, but in drought stress, these correlations became weaker.

Analysis of correlations between traits (Supplementary
For the purpose of further analysis, relative drought effects (RDE) were computed for all lines and all after-harvesting traits; they were on average negative ( Figure 1A). Strong correlation of effects was found for two traits describing numbers of tillers and for subsets of traits describing main and lateral spikes ( Figure 1B). On the basis of traits with largest inter-line variation, i.e., four main spike traits (LSm, NSm, NGm, GWm), we identified two lines that were very distant in optimal conditions and quite similar in drought, therefore, contrasting with respect to reaction to drought: LA (MCam001) with large trait values under control conditions and a large relative loss under drought, and LB (MCam079) with low values under control and small relative loss under drought (Supplementary Figures 1A,B). These two lines are extreme in terms of the scores of the first principal component obtained for RDE of main spike traits ( Figure 1C and Supplementary Table 2). Lines LA and LB were contrasting in a similar manner also with respect to traits concerning lateral spikes, TGW (note a positive drought effect for LB, Figure 1A) and, to a smaller extent, number of tillers. Line LB was somewhat taller than LA, but the drought effects for plant height were similar for both lines. Parental lines Maresi and Cam/B1/CI were not as extreme as LA and LB, which indicates existence of transgression effects with respect to reaction to drought.

Image-Based Phenotyping of the Dynamics in Drought Response
Out of 376 characteristics provided by the image analysis pipeline, seven contained no observations, and 83 provided data that did not satisfy the criterion requesting the correlation between carriers to be bigger than 0.2 (Supplementary Table 1). The classification of all image-extracted traits with respect to their type and origin is shown in Table 2. One can note that a rather small proportion of traits derived from near infrared imaging (41.18%) was accepted for further analysis. All traits but one (no. 21, "leaf length") exhibited a significant mean drought effect at least at one time point (ANOVA F-test, p < 0.05, Bonferroni correction over 286 analyzed traits). First significant effects were observed on the fourth day of drought, and the number of traits with a significant effect grew until first day after drought (  size. The largest and consistent drought effects were observed in clusters 1 and 3. Cluster no. 1 contained four highly similar profiles of negative effects for traits measured in visible or fluorescent light describing projected plant area (traits 38, 44, 45, 277). Cluster no. 3 contained six profiles for traits describing proportion of yellow, red or brown pixels (traits 139, 142, 144, 324) and skewness of colors (traits 328, 334), measured in visible light by side and top cameras, with positive drought effect increasing from about 12 to 18 day, meaning increased intensity of these colors under drought (plants less green). Groups 4, 5, and 6 contained profiles for, respectively, 26, 116, and 130 traits, with positive (group 4; all but one were fluorescent traits) and negative (group 6; 90% were geometric traits) drought effects. In terms of drought effects progression, for most of the traits, the maximum effects were observed at the end of drought period, but there was a group of traits in cluster 4 with a large drought effect appearing about 5 days later (traits 124, 132, 166, 167, 284, 286, 297, 301, 302, 308, 309, 317, 344, 345). Clusters 2 and 4 were composed of traits that indicated a lack of recovery of plants by the end of the observation period.  The differences between genetic variance in drought and control, measured by the fold change in logarithmic scale (log 2 FC), were visible, especially after 18th day (end of drought period). The median values of FC (over traits) first slightly decreased, but then increased, especially for texture traits ( Figure 2B), meaning bigger variance among lines subjected to drought for many image-derived features. Clustering of time profiles of genetic correlation between conditions (GC), the indicator of the genotype by environment (GE) interaction, revealed seven groups of traits ( Figure 3B; Supplementary Figure 2B). A temporary decrease of GC after the end of the drought period could be seen in cluster 2 containing five traits (153,160,181,287,288; color uniformity, texture, skewness of saturation). In cluster 4 (28 traits; 46% texture traits, 75% top camera), a very strong decrease of correlation to values close to -1 (indicating a strong cross-over interaction) with a subsequent increase to values close to zero was observed. For traits in group 5 (34 traits; 67% color, 2 geometric traits: compactness 245, area/skeleton length 272), the correlation started decreasing at the end of drought period. For three traits in group 1 (43, 44, 45; projected plant area), the correlation was small on all days; for the rest of the traits, in groups 6 (80% geometric traits, 82% top camera) and 7 (64% geometric, 62% side camera), the correlations were close to 1 on all days, indicating a lack of GE interaction.
Simultaneous classification of traits with respect to time profiles of both drought effects and genetic correlation provided 21 combined clusters ( Table 3). Among clusters with a moderate or a large number of traits, one can note cluster (4, 4) of 14 traits characterized by positive drought effects, no recovery and genetic correlation decreasing after drought period, approaching -1 during recovery, and vanishing at the end of observations. All traits in this cluster were obtained from fluorescent images and describe color (eight traits, yellow or red fluorescence intensitypositive effect means more red or yellow; color or brightness uniformity -positive effect means less uniform) or texture of plants (six traits, positive effect means coarser surface). Profiles of mean values and drought effects for an exemplary trait from this group, trait no. 124 (yellow fluorescence intensity), for all lines, are shown in Figure 4A; characteristic is a late appearance of positive drought effects (increased yellowness in drought),  Cluster numbers correspond to the ones shown in Figure 3.
no recovery for most of the lines, and differentiation of the effects among lines leading to decreased genetic correlation between conditions. We also note cluster (6, 6) of 36 traits with negative drought effects, but with no evident decrease of genetic correlation (all traits in this group but one are geometric ones, measured in fluorescent or visible light by top view camera); profiles for a trait from this group, no. 184 (diameter of the smallest circle drawn around the plant), shown in Figure 4B, show that the drought effects were relatively homogeneous for lines, which led to a rather large genetic correlation. The most numerous classes, clusters (5, 7) enriched in texture traits and (6,7) enriched in geometric traits, contain traits with moderate drought effects and no interaction of genotypes and environment. The behavior of lines LA and LB with respect to traits shown in Figures 4A-C can be described as follows. The (relatively) less resistant line LA showed a larger "yellowing" effect [trait 124, cluster (4,4)] during drought (days 13-14) than the more resistant line LB; after drought, the effects for LA were much smaller, and LA, unlike LB, showed recovery by the last day of observations. The smallest circle diameter [trait 184, cluster (6,6)] was higher for LA than for LB both in control and drought conditions, and the negative drought effects (smaller diameter in drought), starting during the drought period, were initially larger for LB, then similar for both lines, and, after drought, indicated a reduced recovery for LB. The projected plant area [trait 39 in Figure 4C, cluster (6,7)] was larger for LA than for LB; the negative drought effects for LA was bigger than for LB on days 14-20, but then was smaller, with, again, no recovery for LB. So, for two of the presented traits (124, 39), the reaction of LA during drought could be stronger than that of LB, but the difference between LA and LB was more visible during the re-watering period than during drought.
The above characteristics are further illustrated by selected side view VIS images of lines LA and LB shown in Figure 5. For example, a good recovery for LA and its lack for LB with respect to projected plant area is well-visible. The pictures show differences in developmental stages between lines (with approximate BBCH ratings; see section "Discussion"). LB has reached heading during the drought period; this process affected its behavior with respect to many traits, and could cause apparently smaller drought effects in terms of yellowing (trait 124) or plant area loss (trait 39). LA has started heading at the end of drought period, when heading in LB was completed; this explains smaller drought effects for LA immediately after drought. During re-watering, LA continued heading and showed recovery, whereas LB, being at the flowering stage and going toward maturity, did not show signs of recovery.

Correlations Between Traits in Time and Post-harvest Traits
After computing -for each day of experiment -correlation coefficients between relative drought effects for traits observed in time and for traits observed after harvesting, the numbers of days were counted on which the correlation was significant, and these counts were summed up over time traits (Table 4).
The total numbers of "days with correlation" were the highest for traits describing properties of main spike, GWl and TGW. The  lowest numbers were obtained for plant height and tiller counts. Taking into account this information, and the information on correlations of RDE for traits observed after harvesting (section "Analysis of Post-harvest Traits"), the following detailed analyses were chosen for interpretation of relationships between timecourse and after-harvest phenotyping: a) Analysis of correlations between traits measured in time and the (correlated complex of) traits describing main spikes -traits no. 5-8; the correlation profiles were averaged over four after-harvest traits taken into account, b) Analysis of correlations between traits measured in time and TGW.
For case (a), clustering of profiles of correlation coefficients ( Figure 6A) revealed seven relatively homogeneous clusters, representing various time profiles of correlation; in some clusters, a trend in time is evident despite the fact that some correlations are not statistically significant (but note that smoothed profiles were drawn, with correlations on particular days possibly significant). Correlations decreasing or increasing during or immediately after drought, and then taking opposite values, were observed in, correspondingly, clusters 1 (27 fluorescent traits) and 2 (32 traits; 53% in color traits, 25% RGB), and 7 (58 traits; 81% geometric traits); or approximately constant until a few days after the end of drought period, and then going up or down, in, correspondingly, clusters 3 (32 traits; no overrepresentation)  and 6 (57 traits; 75% top view). Clusters 4 (36 traits; 97% geometric, 100% side view) and 5 (44 traits) were made of traits with rather low and constant correlation with RDE for traits describing main spikes. For case (b), grouping of correlation profiles provided five clusters of traits ( Figure 6B), with clusters 4 (63 traits; 65% fluorescent; 41% texture) and 5 (102 traits; 77% geometric) containing features characterized by, correspondingly, decreasing and increasing correlation in drought period, and going to opposite values afterward. Cluster 2 contained 33 traits (90% geometric), for which correlations with TGW increased during the drought period and decreased afterward. By merging clusters characterized by similar trends in time: (1,2,3), (6,7), (4,5) for correlations with main spike traits, and (4), (2,5), (1,3) for correlations with TGW, we obtained the less fine divisions of traits (Table 5). For most of the traits, the profiles of correlations with spike traits and with TGW had similar characteristics (both increasing, both decreasing or both constant -the diagonal in Table 5). The joint cluster of traits which did not correlate with main spike traits but revealed a correlation with TGW decreasing to negative values contained 25 traits (92% geometric, 80% side view). In addition to correlations of drought effects, we also computed correlations between absolute values of traits measured in time and absolute values

Selection of Traits Based on Correlations of RDE
The joint classification of 286 traits observed in time with respect to the shape of profiles of mean drought effects and of genetic correlations between conditions (section "Image-Based Phenotyping of the Dynamics in Drought Response"), and shape of the profiles of correlations with after-harvest traits (section "Correlations Between Traits in Time and Post-harvest Traits"), provided: 59 groups of traits in case of correlation with spike traits, and 46 groups in case of correlation with TGW. The final selection of traits observed in time during drought stress, with various characterizations with respect to timecourse of mean relative drought effects, genetic correlation (GE interaction), and correlation with after-harvest traits, was made by finding in all joint categories defined above the time traits characterized by the number of correlation days maximal in the given category and bigger than the upper quartile of the distribution of number of correlation days for all traits. This operation provided the list of 21 image-derived traits for which RDE were highly correlated with RDE for spike traits, and 11 traits with the same property but with respect to TGW ( Table 6 and Supplementary Table 4). The set of time traits selected on the basis of correlation with both main spike traits and TGW consists of traits characterizing: plant convex hull circularity (16), plant's convex hull geometry (75), texture (104), proportion of brown pixels in relation to green pixels (144), and hull compactness (256). Characteristics of these traits are summarized in Supplementary Table 5.

Prediction of RDE for Post-harvest Traits Based on Image Phenotyping During Drought
Utilizing image-based traits selected in previous section, we performed partial least-squares regression model fitting to predict RDE for after-harvest traits, separately for each day of observations. The profiles of the measure of model fit, predictive residual error sum of squares (PRESS), for prediction of RDE for four main spike traits using RDE of 21 image traits are shown in Figure 7A. The prediction quality was relatively good when traits measured on days 15-20 were used. Toward the end of imaging period, prediction by 21 selected traits was better than by using all traits. For prediction of TGW, the profile of PRESS corresponding to 11 selected traits shows a steady improvement, whereas the profile for all traits -a temporary one on days 18-19 ( Figure 7B). On days 31-32, quality of prediction for both spike traits and TGW using the selected subset of traits is as good or better than prediction using all 286 traits.
The image-based traits with a dominant role in prediction can be identified in PLS biplots shown in Figure 8. Genotypes more resistant to drought in terms of main spike characteristics (namely, losing less main spike length and grains) were characterized: The above relationship can be also confirmed by comparing the behavior of lines LA (relatively susceptible) and LB (relatively resistant) with respect to above-mentioned image traits (Supplementary Figure 3). The standard deviation of the saturation values of the plant pixels. The lower this value, the more uniform is the saturation of the plant colors.

(Continued)
Frontiers in Plant Science | www.frontiersin.org Standard Deviation -first order texture property (independent of pixel neighbors). Calculated on grayscale image derived from channel Blue (RGB).

Yes
Correlations computed for drought effects.

DISCUSSION
This paper discusses the employment of high-throughput, noninvasive imaging platform to the characterization of phenotypic reaction of barley lines to limited irrigation. We analyzed the data set provided by the image analysis pipeline containing observations of more than 350 traits with the aim of defining a non-redundant set of image-derived phenotypic characteristics useful for a further inference on the fitness of genotypes. The criteria for selection of traits were based on genetic considerations and on relating the picture of RILs observed during drought stress to their characteristics with respect to yield-related traits observed after harvesting.

Analysis of Post-harvest Traits
In our experiment, variance components for yield-related traits of RILs were lower under drought than well-watered conditions. This result confirms observations from practical breeding that differences between genotypes (lines/cultivars) in crop performance features are more evident in favorable conditions than under stress. That is widely used in breeding high-yielding cultivars which are selected most often under optimal conditions (Laingi and Fischer, 1977;Atlin et al., 2017). Genetic correlations between conditions estimated in our experiment showed that the strategy of selection for any environment may work better for traits characterizing spike architecture than for yield or TGW, and will not work for plant architecture (number of tillers).
Our inference was based not on absolute trait values, but on relative drought effects characterizing the reaction of genotypes to water shortage. We utilized these effects for four characteristics of main spikes to identify genotypes with a relative good (like LB) and poor (like LA) resistance to water shortage. We attribute the presence of such genotypes in the studied population to the properties of its parents: Maresi, a good-quality European variety, and Cam/B1/CI, a line of Asian origin. The characterization by relative effects is independent of the fact that LB (and, indeed, Cam/B1/CI) behaved much worse than LA (and Maresi) in optimal conditions.

Phenotyping in Time
The analysis of image-derived data acquired over 32 days of the experiment run on the high-throughput phenotyping platform provided information on the behavior of genotypes during drought period and during re-watering. Low replicability of observations obtained from infrared imaging was revealed; this is in contrast with findings of Chen et al. (2014), who reported good reproducibility of NIR-based traits. Relative drought effects (mean -averaged over all genotypes; positive or negativedepending on the trait nature) were observed for almost all features extracted by the image analysis pipeline; they started to be significant on the fourth day of limited watering. The earliest and largest negative effects were observed for plant area measured in visible or fluorescent light by side-view camera. The largest positive effects, starting in the middle of drought period (similarly as in Chen et al., 2014), were noted for proportions of yellow, red or brown pixels measured in visible light by side-and top-view cameras. So, the first observable effect of reduced watering was the slowed down, and then strongly reduced, plant body development; it was followed by plants turning less green as a result of chlorophyll degradation (Ghandchi et al., 2016). Reducing the chlorophyll content causes a decrease in photosynthesis intensity and thus limits the growth and development of plants (Buschmann and Lichtenthaler, 1998;Buschmann et al., 2000;Jansen et al., 2009).
In addition to relative drought effects, we studied the time profiles of genetic variances and of genetic correlation coefficients between conditions, which informed us about the variability and co-variability (GE interaction) of the traits under different conditions. The genetic variance under drought was bigger than under control conditions in the re-watering period (especially for many texture traits). This is in contrast with what was observed for after-harvest phenotypic traits and suggests that drought resistant genotypes might be selected under limited watering using properly chosen image-derived traits. This possibility is further confirmed by the existence of many traits (especially those measuring color or texture) for which the GE interaction increased toward the end of the observation period.
The simultaneous classification of traits with respect to time profiles of drought effects and of genetic correlation between drought and control showed that numerous traits describing color or texture of plants were characterized by positive drought effects and decreased correlation. It indicates that under drought the leaves discolored (turned brown or yellow), while under optimal watering they were still green, but the effect sizes were different for different genotypes and the correlation between drought and control decreased. On the other hand, many traits describing geometrical properties of plants, measured in fluorescent or visible light by top view camera, showed negative drought effects, but with no evident decrease of genetic correlation between treatments. This again suggests the use of color and texture traits in preference to geometric ones in selection.

Relationships Between Traits in Time and Post-harvest Traits
The numbers of days on which the significant correlation coefficients between relative drought effects for traits observed in time and traits observed after harvesting were the highest for spike characteristics and TGW, while the lowest for final number of tillers and plant height. That was due to the duration of drought in the context of plant development. Phenotyping started when plants were at the stage of tillering (26-29 in the BBCH scale). During the period of water stress and 2 weeks after drought, they have undergone the phases of tillering, and shoot and spike development, which may explain observed associations. On the other hand, re-watering contributed to grain filling and the emergence of new shoots which did not form grains. Plants continued growing also after the period of image phenotyping; therefore, there is no clear relationship between image-derived traits and the final number of tillers and plant height. This result suggests that a different drought application scenario would be necessary to study the reaction of genotypes to drought with respect to plant architecture.
About 60 image traits were selected (mainly color and texture) with RDE characterized by correlation with RDE for post-harvest spike traits and TGW increasing in time (to positive values). As we noted, this was due to the fact that the lines classified by us as "relatively resistant, " e.g., LB, did not recover from the drought by the end of the observation period (positive color/texture effect), whereas the "relatively susceptible" lines, e.g., LA, did recover (negative color/texture effect). For example, we have shown that during drought, the leaves of LB turned yellow to a lesser extent than those of LA, but after re-watering, toward the end of observations, the "yellowing" effect of the LB deepened and remained, while for LA, gradually decreased and become even negative at day 32. It indicates that LA has the ability to produce new shoots and green leaves after re-watering, which reduced the overall "yellowing" effect.
On the other hand, about 100 traits (mostly geometric) had RDE characterized by correlation with RDE for post-harvest spike traits and TGW decreasing over time (to negative values). This was also caused by the lack of regeneration of "relatively resistant" genotypes in terms of plant architecture by the end of the observation period. For example, drought caused a reduction in the projected area of plants (trait 20) and skeleton length (trait 261); the effects were similar for LA and LB, or a bit smaller for LB, until approximately day 20, when recovery started for line LA but not for LB.
We can safely assume that the differences explained above between lines LA and LB, more pronounced during rewatering than during drought, originated from the properties of the parents: Cam/B1/CI is a very early line as compared to Maresi, with about 10 days earlier heading and earlier maturity (Mikołajczak et al., 2016). A corresponding difference between lines LB and LA was visible in image phenotyping (see section "The Trait Selection Procedure" for more discussion of phenology). The lack of recovery of LB could be caused by its inability to regenerate after heading, in a later developmental stage leading to maturity, whereas less yellow appearance of LA at the later stages could be caused by spikes and new leaves emerging after drought period. The inference that we could make on the basis of relative drought effects and identification of "relatively" resistant and susceptible genotypes shows important differences in the behavior of genotypes. While LA showed a "stay green" and recovery behavior, LB seemed to promote senescence, potentially with rapid mobilization of resources into grain filling during drought (thus avoiding drought effects by a shorter life cycle). The shorter vegetative phase of LB, with fewer tillers, may cause lower yield in the non-stress condition.
In addition to correlations of relative drought effects, we presented results showing that image-based traits expressing good correlations with absolute values of most of the afterharvest traits could be found. This finding creates a possibility of using the approach presented in this paper for characterization of genotypes better suited for cultivation in optimal or suboptimal conditions. Such analysis is in progress and will be reported elsewhere in the context of utilization of molecular polymorphism data for localization of dynamic quantitative trait loci.

The Trait Selection Procedure
The classification and selection of traits was based on the assumption that useful image-derived features should provide information on three components of RIL's characterization: environment (drought) effects, genetic correlation between conditions related directly to genotype-by-environment interaction, and the link between behavior during the stress and the final yielding performance. First two parameters are widely used for static characterization of RIL populations observed in traditional experiments with only after-harvesting phenotyping; they were also used for dynamic description of a smaller plant population by Chen et al. (2014) who imaged 18 barley cultivars under well-watered and limited irrigation conditions. Having observations of post-harvest yield-related traits, we were able to extend this idea to correlations between relative drought effects of two sets of traits. By a multiple use of a simple clustering procedure, we divided the set of all traits into disjoint groups with various characteristics, and, finally, selected representatives of these groups valuable for inference. We think that the described procedure can be used extensively in image phenotyping for an unsupervised selection of representative traits obtained from various existing image analysis pipelines. Pipelines used by image phenotyping platforms are continuously modified, and new pipelines are being developed. Their potential to provide new phenotypic characteristics is almost infinite. In such situation, the presented selection procedure can be useful. Linking this procedure to emerging taxonomies of image-based plant traits (Fahlgren et al., 2015;Das Choudhury et al., 2019) will be beneficial for better annotation of data sets.
In our studies, the applied statistical approaches allowed choosing 21 traits out of 286 observed in time for which drought effects were highly correlated with drought effects for spike traits, and 11 traits with the same property but with respect to TGW. The traits common to these two sets characterized plant convex hull geometry related to the growth habit, and texture, proportion of brown pixels and chlorophyll intensityfeatures related to the earliness. In the studied population, growth habit was mainly determined by the semi-dwarfing gene from Maresi, while earliness by gene(s) from Cam/B1/CI (Mikołajczak et al., 2016Ogrodowicz et al., 2017). Segregation of genes determining these two properties largely influenced variation of all the observed traits, both measured in time and after harvest.
In the context of observed phenological differences between studied parental forms and RILs, the problem arises of the influence of these differences on the results of trait selection and on the conclusions. Due to the absence of image-derived traits annotated directly as corresponding to reaching consecutive developmental stages, we attempted to target this problem by visual scanning of (top view RGB) images to score, for each plant, the day of reaching BBCH49 (the stage when flag leaf sheath opens; all data not presented), and using these data as a covariable in the linear mixed model in which drought effects and variance components for RILs were estimated. The differences between significance of mean drought effects in these two models were quite minor (Supplementary Figures 4A,B); the number of significant drought effects in the ANCOVA model was slightly smaller until day 22, which means that for some traits the difference between control and drought conditions was partially explained by differences in phenology. Bigger differences between two models were observed for fold change of genetic variance between drought and control conditions (Supplementary Figures 4C,D). The proposed trait selection procedure performed using results of the ANCOVA model provided 26 traits, out of which 16 were the same as in the procedure based on ANOVA model reported in Results (Table 6), including all five traits selected for both correlation with spike traits and with TGW. Thus, some effect of the fact that RILs were, at a given time point, at different developmental stages, was observed. This points to a need for the development of image recognition algorithms able to score phenology in an automatic way. This additional analysis also further illustrates the flexibility of the proposed trait selection procedure with respect to the statistical model of the experiment. In reference to the discussion on earliness of parental genotypes and RILs, we note that, on average, lines LB and LA reached BBCH49 on imaging days, correspondingly, 6.8 and 25.0 in control conditions, and on days 5.5 and 26.4 under drought, which is consistent with Figure 5. For parental lines Cam/B1/CI and Maresi, the corresponding values were 4.9 and 26.6 in control conditions, and 3.5 and "more than 32" (after the imaging period) under drought, which, in terms of direction of differences, agrees with the results of Mikołajczak et al. (2016).
The final step of our analysis, prediction of RDE for afterharvest traits by RDE of selected image-based traits, was performed to illustrate a possible application of the trait selection procedure. We do not claim that the PLS regression is the best method to perform such prediction. However, we showed that a good fit of a prediction model can be obtained for a small set of selected traits.

CONCLUSION
Image phenotyping allows observation of changes in plants during growth. The results of our experiment show these changes in various aspects: changes in color, texture, and geometrical properties of plants under optimal and water stress conditions, and association of the traits evaluated in time with post-harvest characteristics. In most cases, these changes were expected and they are in line with practical observations. Our results permitted to distinguish features related to the convex hull geometry, texture and proportion of brown pixels for which effects of drought were highly correlated with drought effects for spike traits and TGW. Estimated genetic correlations between drought and well-watered conditions showed that strategy of selection for a broad range of environments may be more effective for traits characterizing spike architecture than for grain yield or thousand grain weight, and can be ineffective in the case of number of tillers. RILs relatively susceptible and tolerant to drought were selected and although the yellowing and area loss effects were somewhat smaller for the latter during drought, after re-watering such effects for the tolerant lines deepened and remained, while for susceptible ones gradually decreased which indicates that in the studied population, tolerant lines escaped drought through accelerated plant development. Further results concerning the described experiment and RIL population will combine the presented barley dynamic phenomic data with highresolution linkage mapping to illustrate the evolution of the genotype-phenotype relationship in time, under water scarcity. Such approach is relevant to identify genes underlying the response of barley plants to drought.

DATA AVAILABILITY STATEMENT
The data set containing images and observations of imagederived and post-harvest traits in the MIAPPE-compliant ISA-Tab format is available online (Kuczyńska, 2020).

AUTHOR CONTRIBUTIONS
AK, MS, and PK conceived and planned the study. AK, KK, KM, and PO performed the experiment. HĆ-K and PK processed the image-derived and post-harvest data and did statistical analysis. MS, TAd, AK, KM, and PO reviewed the literature and provided biological interpretation. AJ, KW-F, and TAl planned and supervised the experiment on the phenotyping platform. SM processed the imagederived data. All co-authors critically reviewed the manuscript and provided constructive feedback on presentation and interpretation of results.