High‐throughput profiling and analysis of plant responses over time to abiotic stress

Abstract Sorghum (Sorghum bicolor (L.) Moench) is a rapidly growing, high‐biomass crop prized for abiotic stress tolerance. However, measuring genotype‐by‐environment (G x E) interactions remains a progress bottleneck. We subjected a panel of 30 genetically diverse sorghum genotypes to a spectrum of nitrogen deprivation and measured responses using high‐throughput phenotyping technology followed by ionomic profiling. Responses were quantified using shape (16 measurable outputs), color (hue and intensity), and ionome (18 elements). We measured the speed at which specific genotypes respond to environmental conditions, in terms of both biomass and color changes, and identified individual genotypes that perform most favorably. With this analysis, we present a novel approach to quantifying color‐based stress indicators over time. Additionally, ionomic profiling was conducted as an independent, low‐cost, and high‐throughput option for characterizing G x E, identifying the elements most affected by either genotype or treatment and suggesting signaling that occurs in response to the environment. This entire dataset and associated scripts are made available through an open‐access, user‐friendly, web‐based interface. In summary, this work provides analysis tools for visualizing and quantifying plant abiotic stress responses over time. These methods can be deployed as a time‐efficient method of dissecting the genetic mechanisms used by sorghum to respond to the environment to accelerate crop improvement.

domesticated 8,000 -10,000 years ago. Thousands of genotypes displaying a wide range of phenotypes have been collected and described (Deu, Rattunde, & Chantereau, 2006;Lasky et al., 2015;Paterson et al., 2009). Sorghum bicolor, the primary species in cultivation today, has many desirable qualities including the ability to thrive in arid soils with minimal inputs, and many end-uses (Morris et al., 2013;Vermerris & Saballos, 2013). For example, grain varieties are typically used for food and animal feed production, sweet sorghum genotypes accumulate non-structural, soluble sugar for use as syrup or fuel production, and bioenergy sorghum produces large quantities of structural, lignocellulosic biomass that may be valuable for fuel production (Murray, 2013;Rooney, 2014). Sorghum genotypes can be differentiated and categorized by type according to these enduses.
Rising interest in sorghum over the last forty years has led to efforts to preserve and curate its diversity. To maximize utility, these germplasm collections must now be characterized for performance across diverse environments (Araus & Cairns, 2014;Fiorani & Schurr, 2013;Furbank & Tester, 2011). Deficits in our understanding of genotype-by-environment interactions (G x E = P, where G = genotype, E = environment, and P = phenotype) are limiting current breeding efforts (Zamir, 2013). Controlled-environment studies are quantitatively robust but are often viewed with skepticism regarding their translatability to field settings. Further, they can often accommodate only a limited number of genotypes at a time. In contrast, field-level studies allow for large numbers of genotypes to be evaluated simultaneously. However, these studies provide limited resolution to resolve the effect of environment on phenotype and often require multiyear replication. This conundrum has motivated enthusiasm for both controlled-environment and field-level high-throughput phenotyping platforms. However, the use of large-scale phenotyping and statistical modeling to predict field-based outcomes is challenging (Deans et al., 2015;Lipka et al., 2015;Zivy et al., 2015).
Here, we sought to define a set of measurable, environmentally dependent, phenotypic outputs to aid crop improvement. We utilized automated phenotyping techniques under controlled-environmental conditions to characterize G x E interactions on a diverse panel of sorghum genotypes in response to abiotic stress. Specifically, we describe and quantify statistically robust differences among the genotypes to nutrient-poor conditions using three phenotypic characteristics: biomass, color, and ion accumulation. Using image analysis to characterize leaf color and biomass over time in conjunction with ionomics, we report measurable, genetically encoded, phenotypic traits that are affected by nitrogen treatment. This work presents a foundation for understanding the range of sorghum early responses to abiotic stress and provides tools for analyzing other available datasets.

| Plant growth conditions
Round pots (10 cm diameter) fitted with drainage trays were prefilled with Profile â Field & Fairway TM calcined clay mixture (Hummert International, Earth City, Missouri), the goal being to minimize soil contaminates (microbes, nutrients, etc.) and maximize drainage.
Before the beginning of the experiment, the thirty genotypes of Sorghum bicolor (L.) Moench (Table S1)  To get useful information from a given image, the plant must be segmented out of the picture using various mask generation methods to remove the background so all that remains is plant material (see Figure 1). A pipeline was developed to complete this task for the sideview and top-view cameras separately and they were simply repeated for every respective image in a high-throughput computation cluster. For this dataset of approximately 90,000 images with the computation split over 40 cores, computation time was roughly four hours. Upon completion, data files are created that contain parameterizations of various shape features and color information from several color spaces for every image analyzed.

| Outlier detection and removal criteria
Each treatment group began with nine reps per genotype for the 100/100 and 50/10 treatment groups and six reps per genotype for the 10/10 treatment group. Outliers were detected and removed by implementing Cook's distance on a linear model (Cook, 1977) that only included the interaction effect of treatment, genotype, and time. That is, for each observation (every image, for every plant, every day), an influence measure is obtained as the difference of the model with and without the observation. After getting a measure for all observations in the dataset, outliers were defined as having an influence greater than four times that of the mean influence and were subsequently removed from the remaining analysis. In total, 5.8% of the data, 1598 images, was removed using this method.

| PCA
Three types of PCAs are generated: one for the shape features, color features, and ionomics. All shape parameterizations that are generated from PlantCV are included in the dimensional reduction. Principal components of color, as defined by the hue channel in two degree increments, are examined using all one hundred and eighty bins in the dimensional reduction. Ionomic PCA was generated using every element that passed internal standards of quality.

| GLMM-ANOVA
Using area as the response variable, a general linear mixed model was created to identify significant sources of variance adjusting for all other sources, otherwise known as type III sum of squares. Designating genotype as G, treatment as E, and time as T, there are six fixed effects: G, E, GxE, GxT, ExT, GxExT. The mixed effect is a random slope and intercept of the repeated measures over time. Wald chi-square statistic was implemented and is a leave-one-out model fitting procedure which allows for adjustment of all other sources.

| Heatmaps
Every cell is a comparison of treatments using a one-way ANOVA wherein the p-value is obtained from a F-statistic generated from the sum of squares of the treatment source of variation. After getting all the raw p-values, a Benjamini-Hochberg FDR multiple comparisons correction is made to aid in eliminating false positives. The p-value distribution was very left-skewed so a log-transform is used to normalize them. Agglomerative, hierarchical clustering was used on the corrected p-values. Each genotype had an associated vector of p-values and a Canberra distance is calculated for all pairwise vectors which are then grouped by Ward's minimum variance method.

| Color processing
PlantCV returns several color space histograms for every image that is run through the pipeline (RGB, HSV, LAB, and NIR). Every channel from each color space is a vector representing values (or bins) from 0 to 255 which are black to full color, respectively. All image channel histograms were normalized by dividing each of the bins by the total number of pixels in the image mask ultimately returning the percentage of pixels in the mask that take on the value of that bin. The hue F I G U R E 1 Experimental Overview. (a) Watering regime used for nitrogen deprivation. The x-axis shows the age of the plants throughout the experiment and the y-axis indicates the estimated volume of water plus nutrients (ml), calculated based on the weight change of the pot before and after watering. Each dot represents the average amount of water delivered each day with vertical lines indicating error (99% confidence interval). Watering regime was increased due to plant age (shades of blue). The experimental treatments are listed above the plots. Volume of water and source of nitrogen are indicated and were scaled based on the 100% (100/ 100) treatment group (1 mM ammonium/14.5 mM nitrate for 100% treatment group). (b) Image analysis example (genotype NTJ2 from 100/100 treatment group on day 16 is shown). Top row: Example original RGB image taken from phenotyping system and plant isolation mask generated using PlantCV. Bottom row: two examples of attributes analyzed (area and color). Scale bar = 15 cm channel is a 360-degree parameterization of the visible light spectrum and contains the number of pixels found at each degree. The colors of most importance are between 0 and 120 degrees which correspond to the gradient of reds to oranges to yellows to greens.
Colors beyond this range, like cyan and magenta, have values of all zeros and are not shown. Means and 95% confidence intervals were calculated on a per-degree basis over the replicates. Area under the curve calculations were made using the trapezoidal rule within the two ranges of 0 to 60 degrees and 61 to 120 degrees which are designated as yellow and green peaks, respectively.

| Ionomic profiling and analysis
The most recent mature leaf was sampled from each plant on day 26 of each experiment, placed in a coin envelope, and dried in a 45°C oven for a minimum of 48 hrs. Large samples were crushed by hand and subsampled to 75 mg. Subsamples or whole leaves of smaller samples were weighed into borosilicate glass test tubes and digested in 2.5 ml nitric acid (AR select, Macron) containing 20 ppb indium as a sample preparation internal standard. Digestion was carried out by soaking overnight at room temperature and then heating to 95°C for 4 hrs. After cooling, samples were diluted to 10 ml using ultra-pure water (UPW, Millipore Milli-Q). Samples were diluted an additional 5 9 with UPW containing yttrium as an instrument internal standard using an ESI prepFAST autodilution system (Elemental Scientific). A PerkinElmer NexION 350D with helium mode enabled for improved removal of spectral interferences was used to measure concentrations of B, Na, Mg, Al, P, S K, Ca, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Rb, Mo, and Cd. Instrument reported concentrations are corrected for the yttrium and indium internal standards and a matrixmatched control (pooled leaf digestate) as described (Ziegler et al., 2013). The control was run every 10 samples to correct for elementspecific instrument drift. Concentrations were converted to parts per million (mg analyte/kg sample) by dividing instrument reported concentrations by the sample weight.
Outliers were identified by analyzing the variance of the replicate measurements for each line in a treatment group and excluding a measurement from further analysis if the median absolute deviation (MAD) was greater than 6.2 (Davies & Gather, 1993). A fully random effect model is created for every element and partial correlations are calculated for treatment, genotype, and the interaction using type III sum of squares.

| RESULTS
3.1 | Phenotypic effects of nitrogen treatment on a sorghum diversity panel Next to water, nutrient supply (most notably nitrogen availability) is often cited as the most important environmental factor constraining plant productivity (Chapin, Bloom, Field, & Waring, 1987;Liu et al., 2015). The initial goal of our experimental design was to enable the early detection and quantification of stress responses in plants. Figure 1 illustrates the overall experimental design we used to test the phenotypic effects of nitrogen treatment on sorghum over the course of a three-week-long experiment using high-throughput phenotyping. Three nitrogen treatments were designed to analyze the effects of source (i.e., ammonium vs nitrate) and quantity of nitrogen on plant development over time (Figure 1a, methods). For this study, sorghum was chosen for its genetic diversity and wide range of abilities to thrive under semiarid, nutrient-limited conditions. In order to test the role that genotype plays in response to nitrogen treatment, a panel of 30 sorghum lines was assembled (Table S1). This panel includes sorghum accessions from all five cultivated races (bicolor, caudatum, durra, guinea, and kafir), representing a variety of geographic origins and morphologies Kimber, Dahlberg, & Kresovich, 2013). The genotypes also display a range of photoperiod sensitivities and are categorized into three general pro-  Figure 2a). Our results indicated that "area" was the plant shape feature that displayed the largest treatment effect. We consider area measurements from plant images as a proxy for biomass measurements as these traits have been shown to be correlated for a number of plant species, including sorghum (Fahlgren et al., 2015;Neilson et al., 2015). Additionally, the effect of low nitrogen on plant color is well established and RGB image-based methods have been described to estimate chlorophyll content of leaves (Cendrero-Mateo et al., 2016;Hu et al., 2010;Junker & Ensminger, 2016;Mishra et al., 2016;Shibghatallah et al., 2013;Wang, Wang, Shi, & Omasa, 2014). In contrast to shape, PCA of color (hue and intensity) attributes at the end of the experiment only separated the high nitrogen treatment group away from the two lower nitrogen treatment groups ( Figure 2b). These data indicate that the different nitrate concentrations in the two lower nitrogen treatment groups significantly affect shape but not color. To further explore the effect that our experimental treatments had on the measured shape characteristics and color for each individual genotype, an interactive version of the generated data is available here: (http://plantcv.danforthcenter. org/pages/data-sets/sorghum_abiotic_stress.html).
Many factors contribute to the ability of plants to utilize nutrients and presumably, much of this is genetically explained. Correspondingly, genotype was a highly significant variable (pvalue = .003 when measuring area) within this dataset. To investigate how much nitrogen treatment response is explained by major genotypic groupings, we calculated the contribution of type, photoperiod, or race on treatment effect. Of these, photoperiod was the only grouping that significantly contributed to area (Fig. S2).

| Size and growth rate during nitrogen stress conditions
Nitrogen stress tolerance is a plant's ability to thrive in low nitrogen conditions. To identify sorghum varieties tolerant to growth in nutrient-limited conditions, we considered plant size at the end of the experiment within the most severe nitrogen deprivation treatment group for all genotypes (Figure 3a)   | 5 Therefore, we hypothesized that lines would be late-responding either (i) because they were proficient at using any level of available nitrogen or (ii) because they grew slowly regardless of quantity of nitrogen supplied. We found that the early-responding lines were larger, on average, than the late-responding lines within the 100/ 100 treatment group (Figure 4b, bottom panel), suggesting that these lines are more competent at using available nitrogen. A subset of these genotypes are displayed in Figure 4c to illustrate our observations. The genotype Atlas is an example of a very early-responding line, and it was one of the largest plants in the 100/100 treatment group, but also one of the worst-performing lines in the 10/10 treatment group (Figures 3a,b and 4c). In contrast, China 17 performed relatively well under nitrogen-limited conditions (10/10), but when nitrogen was abundant (100/100), the biomass accumulation was relatively poor (Figures 3a,b, and 4c). A similar phenotype was observed for PI_510757. In addition to varying the amount of nitrogen available, we also tested whether any lines harbor a preference for nitrogen source. Nitrogen is typically available in two ionic forms Notably, Atlas translated an increased ammonium level into larger plant size. In contrast, San Chi San showed no change in average plant size between the two lower nitrogen treatments (Figure 4c).
Among the 30 tested genotypes, 16 displayed little difference between the 50/10 and 10/10 groups in terms of plant size toward the end of the experiment (Fig. S3). This highlights the importance of considering both quantity and source when investigating nitrogen responses.

| Combined size and color analysis over time
In addition to affecting shape attributes, nitrogen starvation generally results in reduced chlorophyll content and increased chlorophyll catabolism. Other groups have used image analysis to estimate chlorophyll content and nitrogen use in rice (Wang et al., 2014). The RGB images contain plant hue channel information, and this was found to be a separable characteristic within the nitrogen deprivation treatment groups (Figure 2b). We assessed color-based responses to nitrogen treatment in the early-and late-responding genotypes as defined in Figure 4a (Figure 5). To facilitate this analysis, we used the generated histograms of images of the individual plants from each day of the experiment and averaged those from the early and late categories within each treatment group (Figure 5a, day 13). We found that the histograms of the plant images contained two primary peaks: yellow and green. For both early-and late-responding lines, the yellow peak was larger than the green for the plants in the 10/10 treatment group as compared to the 100/ 100 treatment group. Early-responding lines within the 100/100 treatment group displayed the largest green channel values. Lateresponding lines grown under nitrogen-limiting conditions displayed the largest yellow channel values. In order to further visualize color-based treatment effects, we subtracted the 10/10 histograms from the 100/100 histograms and plotted this difference (Fig. S4). This revealed that although the late-responding lines were more yellow, the magnitude difference from the treatment was similar for early and late lines in the yellow channel. In contrast, the early-responding lines tended to have a larger green channel difference between the 10/10 and the 100/100 treatment groups, with early-responding lines showing a larger difference in the green channel.
To assess color-based treatment effects over time, we took the area under the histograms (e.g., Figure 5a) for all time points and plotted them against plant age (Figure 5b). Given the peaks within the histograms mentioned above, we focused on these regions and defined yellow (degrees 0-60) and green (degrees 61-120) to facilitate quantitative analysis. As expected, plants within the 10/10 treatment group were generally more yellow (and consequently less green) over the course of the stress treatment. We detect a peak difference between yellow and green occurring on day 13, and then the effect diminishes. A similar peak and overall pattern is seen in the 100/100 treatment group, with plants greening after day 13.
Focusing on either the green or the yellow hue, there was no discernable difference in color over time between late-and earlyresponding lines within the 10/10 treatment group (left panel, and green (degrees 61-120) hues over time for 100/100 (left) and 10/10 (right) treatment groups. Plotted is the area under the curves presented in A (y-axis) over the duration of the experiment (x-axis) for early-and late-responding genotypes. Gray areas indicate standard error dotted versus dashed lines, p > .05). However, within the 100/100 treatment group, early-responding lines were consistently greener, while late-responding genotypes became increasingly yellow until day 13. Late-and early-responding lines behaved differently under the 100/100 nitrogen treatment conditions, becoming significantly different quickly (day 10, p < .05) and remaining so for the duration of the experiment, with the most significant difference occurring on day 13 (p < 1 9 10 À15 ).
Combining the above plant size-and color-based data, we conclude that the "early-responding phenotype" indicates that these plants are able to take better advantage of available nutrients.
Importantly, both size and color phenotypes indicate that the earlyresponding genotypes do not display the fitness advantage in low nitrogen conditions. Together, these data demonstrate that colorbased image analysis is consistent with and complementary to the more-established biomass measures of fitness and performance.
Phosphorus (mg)/kg tissue In addition to the image-based analysis used above to reveal measurable size-and color-based outcomes in response to nitrogen treatment, we also performed ionomic analysis to gain better insight into the physiological changes that occur in response to nitrogen (Fig. S5). It has been established that both genetic and environmental factors and their interactions play a significant role in determining the plant ionome (Asaro et al., 2016;Baxter & Dilkes, 2012;Baxter et al., 2008;Chao et al., 2012;Shakoor et al., 2016;Thomas et al., 2016). Thus, this analysis was used to explore alterations that might not be revealed by shape or color analysis but would still contribute to the effect of nutrient availability. Each element was modeled as a function of both genotype and treatment, and genotype was a significant factor for most elements with Mo, Cd, and Co being the most affected by genotype (Figure 6,  Figure 6a). In Arabidopsis, the presence of nitrate has been shown to inhibit phosphorus uptake (Kant, Peng, & Rothstein, 2011;Lin, Huang, & Chiou, 2013). Consistent with this, dry weight-based concentrations of phosphorus were inversely proportional to administered nitrate treatment, with the 100/100 treatment group accumulating less on average than either 50/10 or 10/10 ( (Coskun, Britto, & Kronzucker, 2016;Jackson & Reynolds, 1996). Additionally, root architecture is affected by nitrogen source and nutrient availability. It has been shown for a number of species, including maize and barley, that ammonium causes a reduction in lateral root branching that can be reversed with the addition of phosphorus (Drew, 1975;Giles et al., 2017;Ma, Zhang, Rengel, & Shen, 2013;Thomas et al., 2016).
Compounding this equation, ammonium also causes acidification of the soil, which affects the uptake of other nutrients and likely alters the root microbiome, further complicating most analysis.
Under the tested experimental conditions, some genotypes were more affected by nitrogen source in terms of end-biomass than others, for example, the difference between San Chi San and Atlas  (Kimber et al., 2013). The largest such institution, the US National Sorghum Collection (GRIN database), provides agronomic characteristic information for 40-60% of the collection (e.g., growth and morphology characteristics, insect and disease resistance, chemical properties, production quality, photoperiod in temperate climates). Thus, much work is yet to be done to fully characterize and maximize the potential of this hearty, productive crop species.
Nitrogen-use efficiency is traditionally defined by the difference in biomass or grain production between plants grown in resourcesufficient versus resource-limited conditions at the end of the growing season. Stated differently, this measure asks the question: How efficient is a plant at translating a provided resource (nitrogen) into plant biomass. Equally important is the ability to efficiently use a limited resource. Factors that play into these distinct definitions of resource-use efficiency include ability to survive periods of extreme stress and rapid utilization of resources as they become available. In this manuscript, we make progress toward deconstructing the building blocks that make up nitrogen response phenotypes. These analyses reveal diverse quantitative indicators of abiotic stress and genotypic differences in stress mitigation that can be used to further crop improvement. Having made progress toward deconstructing these building blocks, we are now in a position to discover the underlying genetic explanations for genotypic variability in resourceuse efficiency and tolerance to resource-limited growth conditions. This work forms a foundation for future research to overlay additional abiotic and biotic stress conditions to achieve a holistic view of sorghum G x E phenotypes. The overall goal of this research was to support such efforts and expedite the process of meaningful crop improvement.

| CONCLUSION
Plant stress tolerance is important for food security, and sorghum has potential as a high-yielding, stress-tolerant crop. "Resourceuse efficiency" is often measured in one of two ways: (i) a comparison between yield production under resource-sufficient and resource-limited conditions, (ii) the ability to survive within resource-limited environments. Here we describe and apply highthroughput phenotyping methods and element profiling to sorghum grown under variable nutrient levels. We quantify nitrogenuse efficiency in genetically diverse sorghum accessions based on color fluctuations and growth rate over time and elemental profile. Through this analysis, we report a time-efficient, robust approach to identifying resource-use-efficient and abiotic stresstolerant plants.