Spatial and Texture Analysis of Root System distribution with Earth mover’s Distance (STARSEED)

Abstract

Purpose

Root system architectures are complex and challenging to characterize effectively for agronomic and ecological discovery.

Methods

We propose a new method, Spatial and Texture Analysis of Root SystEm distribution with Earth mover’s Distance (STARSEED), for comparing root system distributions that incorporates spatial information through a novel application of the Earth Mover’s Distance (EMD).

Results

We illustrate that the approach captures the response of sesame root systems for different genotypes and soil moisture levels. STARSEED provides quantitative and visual insights into changes that occur in root architectures across experimental treatments.

Conclusion

STARSEED can be generalized to other plants and provides insight into root system architecture development and response to varying growth conditions not captured by existing root architecture metrics and models. The code and data for our experiments are publicly available: https://github.com/GatorSense/STARSEED.

Introduction

Studying plant roots is one of the keys to achieving the second Green Revolution [1]. The study of plant roots requires, among other things, effective characterization and development of comparative methods for root development, architecture, and spatial distribution, including early in the root system’s development. Current methods, such as the widely used $\hbox {WinRHIZO}^{\hbox {TM}}$ and $\hbox {WinRHIZO}^{\hbox {TM}}$ Tron software suite (Regent Instruments Inc., Quebec, Canada) use 2D images of roots to measure individual parameters related to morphology and topology. These standard software packages have also been developed to further make accessible the analysis of individual root traits [2, 3]. These tools capture information relevant to characterize root spatial distribution (RSD). However, standard software packages typically provide parameters that are only capturing one or a few specific aspects of a root system such as total root length, surface area, branching angle and root order. These packages are less effective at holistically quantifying intact root systems and how they change in response to edaphic conditions. To address this shortcoming, machine learning can be used in tandem with these biological features to automate and improve understanding of root systems [4].

A subarea of machine learning, deep learning, has been widely used for automating root detection and segmentation from 2D images of root systems [5,6,7,8,9]. Deep learning methods can require high computational power and large amounts of annotated data however, as well as lack explainability [10, 11]. Yet these techniques have shown promise for localizing and identifying unique 2D texture features useful for root phenotyping [6]. Texture, meaning the spatial arrangement of the pixel values in a raster grid [12], is a powerful cue that can be used to identify patterns tied to RSD characterization. A straightforward approach to compute texture from an image is to use counts (i.e., histograms) of local pixel intensities [13]. This representation of texture can be defined as counting the percentage of root pixels in an area of the image. By counting the number of root pixels in each cell of a raster grid, this allows for quantification of localized pixel intensity depending on the cell size of the grid. In addition to this baseline approach, one can compute more complex texture features such as fractal dimension and lacunarity.

Root systems can be described as approximate fractal objects over a finite range of scale. A fractal object is defined by two fundamental characteristics: self-similarity (i.e. shape variations of the object on one scale are repeated at another) [14] and the possession of a fractal dimension that is expressed as a non-integer dimension, as opposed to the more familiar integer Euclidean dimensions [15]. This fractal dimension value can be calculated for root systems using a specific approach called the box-counting method [16, 17] and has been used successfully to capture root system architecture [16, 18, 19].

A related feature to fractal dimension is lacunarity. Lacunarity captures the distribution and size of “gaps” in an object [15, 20]. An object, in our case the image of a root system, will be highly lacunar if it contains many gaps distributed over a greater range. Lacunarity will thus be smaller in value when root systems are dense and increase with the presence of gaps (i.e., areas of the image where roots are not located or areas of the growth media matrix where roots are not located) and a coarse spatial arrangement of roots [15, 21]. Since some visually distinct images can have the same fractal dimension, lacunarity provides a feature that can aid in discriminating between these distinct textures [21,22,23].

The aforementioned texture features are promising candidates in building a new tool to accurately and precisely describe RSD. Therefore, a novel algorithm is needed that would allow for the calculation of these texture features from 2D root images and spatially explicit comparisons between root systems based on the feature values. For these aims, Earth Mover’s Distance (EMD) appears to be a promising candidate approach [24]. EMD, also known as the Wasserstein-1 distance [25], is used to determine quantitative differences between two distributions [24]. EMD has several advantages such as allowing for partial matching (i.e., comparisons can be made between representations of different sizes such as comparing smaller and larger root systems) and matches our human visual perception when the chosen ground distances (i.e., distance between feature vectors) is meaningful [24]. Essentially, we want to use EMD to find the minimal distance or amount of “work” to transform one root spatial distribution to another.

In this study, we developed the Spatial and Texture Analysis of Root SystEm distribution with Earth mover’s Distance (STARSEED) approach, to better characterize RSD and overall root soil exploration. Our approach used EMD to quantify the heterogeneity among root structures. The proposed STARSEED approach allowed us to establish meaningful biological connections between the treatments such as cultivar and moisture level, and the observed RSD as well as provide detailed qualitative and quantitative analysis. The aims and contributions of STARSEED are the following: (1) characterize the spatial arrangements of roots in local (i.e., smaller regions of the image) and global (i.e., whole image) contexts, (2) extract useful information to describe the distribution of roots in the image, and (3) give precise insight for biological interpretation.

Materials and methods

Greenhouse setup and data collection

Data presented and used in this work is coming from a previous study that has been published [26]. A summary of the greenhouse setup and data collection of the materials and methods are presented here. Sesame was grown in 64 custom-made rhizoboxes. Each rhizobox was filled with 1550g of inert calcine clay (Turface Athletics, Buffalo Grove, IL, USA), henceforth referred to as soil. A single seed from one of four non-dehiscent sesame cultivars, all provided by Sesaco Inc. (Sesaco32, Sesaco35, Sesaco38 and Sesaco40, referred to as S32, S35, S38 and S40 thereafter), was planted per box. Four soil water content treatments were implemented: 60, 80, 100 and $120\%$ of the soil water holding capacity, corresponding to 688, 918, 1147 and 1376 mL of water per rhizobox, respectively. The soil and water were mixed thoroughly together before filling the rhizoboxes to promote homogeneity of the soil water content throughout. The top of each rhizobox was covered with Press’N Seal Cling Wrap (Glad, Oakland, CA, USA) to minimize water evaporation. No water was added afterwards throughout the experiment with the assumption that soil water content, in the absence of evaporation, stayed relatively constant throughout the 16–21 days of each run experiment. A small hole was pierced in the plastic wrap upon seedling emergence to allow for leaf and stem growth. The two factors were arranged in a complete randomized design with 4 replications in a greenhouse where the daily temperature was maintained between 25 and 35 degrees Celsius.

Four runs of 64 rhizoboxes were completed. Run 1 and 2 were prepared with soil and water only, while run 3 and 4 were fertilized with 1.51mg of 15-5-15 + Ca + Mg Peters Excel mix fertilizer (ICL Specialty Fertilizers, Summerville SC) and 0.29g of ammonium sulfate were dissolved in the water applied to each rhizobox prior to mixing with the soil. Baking soda was added as needed to neutralize the fertilizer solution to a pH of 6. Plants were grown for a duration of 21 days after planting (DAP) for the runs without fertilizer and 16 DAP for the runs with added fertilizer due to their fast growth.

On the last day of the experiment, all rhizoboxes were scanned with a Plustek OpticSlim 1180 A3 flatbed scanner (Plustek Inc., Santa Fe, CA, USA) to generate the raw images with the size $6800\times 4676$ at a resolution of 400 Dots Per Inch (DPI). Rhizoboxes were all set on the scanner the same way using a physical guide so that all rhizoboxes had the same position in the scanned images. For each run, images from seeds that did not germinate or seedlings that died during the experiment were not considered in the analysis; 11 images were thus removed from Run 1, 10 from Run 2, 11 from Run 3 and 5 from Run 4. The total number of analyzed images was 107 Runs 1 and 2, and 113 for Run 3 and 4. Representative RGB image samples for each cultivar and water level are shown in Additional file 3: Figure S3 and Additional file 4: Figure S4. The roots from all analyzed images were hand-traced with $\hbox {WinRHIZOTron}^{\hbox {TM}}$ (Regent Instruments Inc., Quebec, Canada) and total root length (TRL) was obtained.

Spatial and Texture Analysis of Root SystEm distribution with Earth mover’s Distance (STARSEED)

The overall pipeline of our proposed method is illustrated in Fig. 1 and explained in more detail in Algorithm 1. STARSEED consisted of five major steps. The first step was to take the input images and perform preprocessing to isolate the roots pixels from the background to focus on the architecture. For step 2, we divided the image into equal-sized bins by generating local regions of the image to improve computational efficiency and incorporate spatial information. Step 3 involved extracting texture features to describe the root pixels within each bin and generate a compact representation of the image (i.e., signature). The fourth step used EMD to measure the magnitude and direction of changes between each root signature. Lastly, we projected the pairwise distance matrix to perform qualitative and quantitative analysis of the root images and their corresponding treatment. Details and rationale behind each step of this algorithm are provided in the following subsections.

Fig. 1

Overview of the proposed STARSEED method step by step

Image Preprocessing—Step 1 The initial step of the proposed approach was to perform preprocessing of the images to isolate the root pixels from the background. Raw images were manually traced using $\hbox {WinRHIZO}^{\hbox {TM}}$ Tron software to generate the labels for root pixels in each image. The black and white binary images were reconstructed based on these labels, in which pixels that corresponded to roots were assigned a value of 1 in the mask and the non-root pixels were assigned a value of 0. The binary images were then cropped and calibrated, and downsampled by a factor of eight through average pooling and smoothed by a Gaussian filter ($\sigma =1$) to mitigate noise before extracting texture features.

Sub-region Generation and Feature/Signature Extraction—Steps 2 & 3 In order to apply EMD to characterize and compare RSDs using STARSEED, we constructed a signature representation, a composition of texture features computed spatially, for each image. Since texture is undefined at a single point [28], a local neighborhood needed to be identified to compute the texture features comprising the signature. In our approach, we divided each image into bins of equal sizes using a grid, each bin serving as the local neighborhood. Within each bin, texture features were computed to locally characterize the RSD. The size of the grid/bins provides trade-offs between computational efficiency and localization. For example, a larger grid yielding large bins corresponds to less computation cost, but each texture feature is computed over a larger area, resulting in a loss of localization. The computational time for the method scaled linearly with the number of bins used (e.g., more bins corresponded to increased computational time) and this information is captured in Additional file 1: Figure S1 and Additional file 2: Figure S2. For the proposed method, the grid size that corresponded to the maximum CH score was used to further analyze STARSEED as discussed in the following Experiments and Results section.

Fig. 2

For step 2, the image is divided into B equally-sized bins after preprocessing. Next, the feature/signature extraction for step 3 is completed by representing each bin as a $\mathbb {R}^{3}$ vector that consist of the bin center location (x and y) and the texture feature (f) describing the pixels in the bin. The set of these B vectors is the final signature representation. The signature can also be represented as a heatmap that shows the feature value for the pixels in a region. In this example, the percentage of root pixels is used as the texture feature

Each bin of the image was represented through a cluster representative. The cluster representative consisted of spatial coordinates (i.e., horizontal and vertical position of the bin center) and a weight to create a $\mathbb {R}^3$ vector to represent a local spatial bin of the image. The weight given to each bin was computed by the texture feature extracted from the pixels in the bin. An example of steps 2 and 3 is shown in Fig. 2. We calculated three feature values in each bin: percentage of root pixels, fractal dimension, and lacunarity. To calculate the percentage of root pixels, we computed the total number of root pixels present in each bin divided by the area (i.e., total number of pixels) of that same bin. The percentage of root pixels provides direct insight into the distribution of roots in the image. If a region of the image has denser roots, the percentage of root pixels in the bins present in that region will be higher. The final representation of the image was a signature, which was the set of clusters from the image. By constructing our signature in this manner, we incorporated spatial (locally at the bin level and globally over the whole image) and texture information, to represent the root distributions.

Novel Earth Mover’s Distance Application—Step 4 EMD can be used to effectively evaluate certain characteristics between different distributions, including images [24]. Images are comprised of pixels and these pixels can be clustered or assigned to meaningful groups based on shared characteristics such as spatial location. The set of these clusters are used to form a signature, a more compact representation of the image to increase computational efficiency [24]. Given an image with C clusters, the signature representation, P, is given by Eq. 1 where $\textbf{p}_i$ is the cluster representative and $w_{\textbf{p}_i}$ is the weight of the cluster:

$$\begin{aligned} P = \{(\textbf{p}_1,w_{\textbf{p}_1}),...,(\textbf{p}_C,w_{\textbf{p}_C})\}. \end{aligned}$$

(1)

Typically, the cluster representative is a feature vector and the weight is the percentage of pixels in a cluster expressing that feature [24]. While the selection of the cluster representative is application dependent, defining the cluster representative as the information/descriptor (i.e., features) and importance (i.e., weight) provides a clear interpretation for EMD within disparate areas of an image. In our proposed STARSEED method, the signature representations were generated as mentioned in Step 3.

Once the signature representations are constructed, we used EMD to generate a pairwise distance matrix $\varvec{D}$ between each pair of RSD signatures, choosing one RSD signature as the reference within each pair. After determining the magnitude of changes between two clusters, the minimal flow matrix, $\textbf{F}$, is used to compute EMD [24]. Given two image signatures, P and Q, with S and T representatives respectively, EMD is formulated in Eq. 2 where $d_{ij}$ is the ground distance between the centroid of regions $\textbf{p}_i$ and $\textbf{q}_j$ and $f_{ij}$ is the optimal flow between $\textbf{p}_i$ and $\textbf{q}_j$:

$$\begin{aligned} EMD(P,Q) = \frac{\sum _{i=1}^{S}\sum _{j=1}^{T}d_{ij}f_{ij}}{\sum _{i=1}^{S}\sum _{j=1}^{T}f_{ij}} \end{aligned}$$

(2)

EMD captures the dissimilarity between spatial distributions: larger values indicate more dissimilarity or “work” to move the defined “earth” feature to the cluster representative. EMD allows one to measure the global change between two spatial distributions (i.e., distance measure) as well as local changes between the two sources of information through the flow matrix, $\varvec{F}$. If we can translate a 2D root image into a spatial distribution of the values for a given texture feature calculated from this image, EMD would then be able to measure the changes in this distribution from one root system image to the next. Consequently, EMD would provide a two-level comparison between these distributions: one at the global scale (i.e., the whole distribution or the whole image) that we can term “holistic” through the EMD measure, and one at the local scale through the flow matrix that shows the distinct magnitude and direction of changes between each image.

In our experiments, we first computed EMD between each pair of images, separately between Runs 1–2 and Runs 3–4. Regions containing no roots were removed from our calculations to mitigate the impact of background on our results. Therefore, we were comparing root signatures of different sizes. In order to not favor smaller signatures, the normalization factor was added to the EMD calculation as shown in Eq. 2. Any distance can be used to define the ground distance, $d_{ij}$ [24]. To consider the root spatial distribution, we compared the center location of each spatial bin by selecting the Euclidean distance as our ground distance to compute $d_{ij}$.

Fig. 3

Example of qualitative results produced by projecting EMD matrix through MDS for the percentage of root pixels as the feature and 2000 spatial bins for images from Runs 3 and 4. MDS served as a visualization approach to show the relationships in the pairwise distance matrix, D. The RSD images are arranged in the 2D space such that images with similar RSD are near one another (i.e., similar MDS coordinates). Additional analysis can be completed by incorporating the moisture or cultivar information. The CH index (step 6) is computed to assess the MDS projection. The different colored frames represent the different moisture levels. The CH indices for 60%, 80%, 100%, and 120% are 8.67, 9.59, 0.23, and 55.68 respectively. The overall CH index for this projection is 15.57. The 120% water level is the most different from the other architectures and the 120% water level is the most compact and separated class visually as well as quantitatively through the CH index

Projection Using MDS—Step 5 Once the distances between each image signature representation were computed to generate the pairwise distance matrix, we can use a dimensionality reduction approach to project the matrix into two dimensions. We projected the data into two dimensions to qualitatively and quantitatively evaluate the pairwise distance matrix. We used multi-dimensional scaling (MDS) [27] as this method had been shown to work well with EMD to identify patterns among images that share some characteristics [24, 29].

To illustrate step 5 in the algorithm, Fig. 3 provides a visualization of the EMD matrix projection in a 2D space through MDS for Run 3 and Run 4 with the percentage of root pixels as the feature calculated for the highest number of spatial bins (2000 in each image). The output from the MDS method is an $\mathbb {R} ^{N \times 2}$ embedding matrix, E, that is a projection of the pairwise distance matrix in the Euclidean space. Once E is computed, the embedding is then used in Step 6 to calculate Calinski-Harabasz (CH) index [30] to quantitatively assess relationships among the root images when considering the treatment values (cultivar and water levels).

Assessment of Relationships for RSD - Steps 6 Lastly, we calculated the CH [30] index to measure the intra-cluster similarity (i.e., RSD representatives of the same cultivar or water level) and inter-cluster dissimilarity (i.e., RSD representatives of different cultivar or water level) across treatment levels in each pair of runs. By using the EMD approach, we were able to score the differences between each treatment level, highlighting: (a) the magnitude and (b) the location of the RSD dissimilarities. Ideally, samples that belong to the same cluster should be “close” (i.e., smaller EMD distances or intra-cluster variance) while samples from different clusters should be “well-separated” (i.e., larger EMD distances or inter-cluster variance). When the clusters are dense and well-separated, the CH index is larger. Finally, a non-parametric Kruskal-Wallis test was performed on the CH scores within each treatment for each pair of runs across features, and then separately for each feature across cultivars or across water levels. When the test was rejected ($p<.05$), Dunnett’s test with Bonferroni’s adjustment was used to compare treatment levels with one another.

Experiments and results

Our dataset consisted of 220 images across four different runs. To select the best bin size of for the experiments, we investigated the impact of coarser (i.e., small number of spatial bins—large grid) and finer (i.e., large number of spatial bins—small grid) scales. The number of spatial bins per image varied from 100 to 2000 in steps of 100 for a total of 20 values. Then, we performed quantitative and qualitative analysis for local and global RSD supported by the STARSEED method. For our local analysis, we calculated the CH index using different features to evaluate the effectiveness of these extracted features at capturing RSD differences between different growing conditions and cultivars. For our global analysis, we aggregated each centered root image to create a representative for each treatment to elucidate the RSD spatial differences between these representative images. We discuss the results of the local and global analysis of STARSEED in the following sections.

Local analysis: relationships to treatment level representatives

The average scores and standard deviations of the CH indices across all number of spatial bins for Runs 1 and 2 are shown in Table 1 while the scores for Runs 3 and 4 are shown in Table 2. The small standard deviations indicate that our method was robust across the sizes of spatial bins (i.e., the scale) used for our analysis. The statistical analysis of the CH scores showed that % root pixels tended to be the feature with the most discriminating power as evidenced by the highest CH score out of the three features considered, except for both fractal dimension and lacunarity on cultivar in Runs 1 and 2. For the rest of this study, we displayed the EMD results for the % root pixels feature, except when presenting results for the cultivar differences in Runs 1 and 2 where we will use the fractal dimension feature. The number of spatial bins that corresponded to the highest CH score is shown in Table 3. In Runs 1 and 2, the CH indices for water levels (1.50 overall average) were not as high ($p<.05$) as for cultivars (3.96 overall average), indicating more overall variability of RSD without added nutrients between cultivars than between moisture levels (Table 1). When looking specifically at each treatment level’s score, both the S32 cultivar and the 80% water level had the largest CH indices within cultivar and water level treatments respectively, indicating these treatments had the most distinct RSD. For Runs 3 and 4 with added nutrients, the trends were inverse in that the average CH indices were higher across all three features and spatial bins for water level (14.70 overall average) than for cultivar (1.62 overall average). The highest average CH score was observed for the 120% soil water level, indicating a distinct RSD for this treatment compared to the three other soil water treatments (Table 2).

Table 1 Average CH Indices and their associated standard deviation for each feature by cultivar and water level for Runs 1 and 2 across all sizes of spatial bins

Table 2 Average CH Indices and their associated standard deviation for each feature by cultivar and water level for Runs 3 and 4 across all sizes of spatial bins

Table 3 Number of spatial bins that scored the maximum CH index for each feature, cultivar, water level in each pair of runs

Global analysis: EMD between treatments

To identify global trends in RSD for each treatment, we aggregated each centered root image to create a representative for each treatment. For each representative root image, we computed the feature representation as shown in Fig. 1 from the average feature values within each spatial bin for images with the same cultivar or water level. After we obtained these representatives, the EMD was computed between each cultivar or water level. By using the EMD approach, we highlighted the differences between each treatment level by displaying (a) the magnitude and (b) the location of the RSD dissimilarities, as shown in Figs. 4 through 8.

Fig. 4

RSD, EMD score, and value changes between S32 (a) and S35, S38 and S40 (b–d respectively) for Run 1 and 2. The color scale for a–d represents the fractal dimension value of each spatial bin. e shows the overall EMD score. f–h Represent the 20% largest value changes between S32 and the three other cultivars; red spatial bins are unique to S32, green spatial bins are unique to S35, S38 or S40, and yellow spatial bins are common to both cultivars

The EMD calculation captured the variability of the RSDs between cultivar and soil water levels within each pairs of runs, as illustrated by Fig. 4. The white arrows (Figs. 4, 5, 6, 7, 8f, g, h) represent the 20% largest value changes between the reference treatment level (Figs. 4–8a) and the other levels (Figs. 4, 5, 6, 7, 8b–d). The EMD scores between the reference treatment level and the other levels are shown in Figs. 4, 5, 6, 7, 8e. Smaller EMD scores indicated higher similarity between the two RSDs.

Differences between cultivars Figure 4 shows the processed RSD images for each cultivar and three EMD results when these cultivars are compared with S32 as the reference cultivar. The EMD score (Fig. 4e) was smallest between S32 and S35, indicating that the two RSDs were most similar. This is confirmed by comparing Fig. 4f, g and h: Fig. 4f comparing S32 to S35 has the shortest and fewest arrows. Arrow abundance and length is a direct representation of the magnitude of the value changes between the two RSD. Lengthier and/or more numerous arrows are therefore indicative of more dissimilar RSDs compared with EMD results showing fewer and/or smaller arrows. Figure 4 shows that the changes in RSDs between S32 and S35 are relatively small, and for the most part are concentrated in the deeper section of the root system. Specifically, S35 seemed to have a more laterally spread out root distribution towards the deepest portion of the image, and a more concentrated root mass along the tap root right above that. Out of the three other cultivars, S32 was most different from S38, and Fig. 4g shows that the majority of the 20% highest value changes in RSD between these cultivars were recorded toward the middle region of the image. S38 also appears to spread out more laterally in the mid-section of the image compared to S32.

Fig. 5

RSD, EMD score and largest value changes between S40 (a) and S32, S35 and S38 (b–d respectively) for Run 1 and 2. The color scale for Figure 5b–d represents the fractal dimension value of each spatial bin. e shows the EMD score overall. f–h represent the top 20% value changes between S40 and the three other cultivars; red spatial bins are unique to S40, green spatial bins are unique to S32, S35 or S38, and yellow spatial bins are common to both cultivars

The differences between S32 and S40 were located in the upper part of the root system, as shown by the location of the top 20% value change between the two RSD in Fig. 4h. S32 tended to have a more spread out root distribution towards the top of the soil surface compared to S40. According to Table 1, the RSD for S40 was different from all three other cultivars’ RSD. The EMD score (Fig. 5) yet showed that the S40 RSD was further from the S32 RSD than the S35 and S38 RSD (Fig. 5e).

Fig. 6

RSD, EMD score and largest value changes between S40 (a) and S32, S35 and S38 (b–d respectively) for Run 3 and 4. The color scale for b–d represents the percentage of root pixels in each spatial bin. e shows the EMD score overall. f–h represent the top 20% value changes between S40 and the three other cultivars; red spatial bins are unique to S40, green spatial bins are unique to S32, S35 or S38, and yellow spatial bins are common to both cultivars

All cultivars appeared to expand laterally with the addition of nutrients as shown in Fig. 6. S40 and S35 were the cultivars with the most distinct RSD as indicated by the highest CH score in Table 2. Figure 6 shows specifically that for S40, the differences between this RSD and the three others remained mostly marginal: the EMD score remained low (< 3.5, Fig. 6e) and the 20% largest value changes between the S40 RSD and the others only represented a few number of changes (Fig. 6f–h). Interestingly, the RSD for S40 was most dissimilar to the RSD for S32 both with or without added nutrients.

Differences between moisture levels For Runs 1 and 2, when considering soil water levels, the RSD were overall very similar as shown in Table 1. The RSD that was most distinct from the others was that of the 80% moisture treatment (Fig. 7). The 80% moisture treatment was overall most different from 100%, and was the most similar to 120% and the 60% treatment RSD (Fig. 7e). Specifically, the 80% treatment had a visibly denser central region than the three other moisture levels (Fig. 7a through Fig. 7d). This is reflected more clearly in Fig. 7f through h: the 80% treatment showed a more concentrated root distribution closer to the top of the root system compared to the other treatments, and so the arrows tend to point toward the upper center of the root system.

Fig. 7

RSD, EMD score and largest value changes between 80% (a) and 60%, 100% and 120% (b–d respectively) for Run 1 and 2. The color scale for Figure 7b-d represents the percentage of root pixels in each spatial bin. e Shows the EMD score overall. f–h represent the top 20% value changes between 80% and the three other soil water levels; red spatial bins are unique to 80%, green spatial bins are unique to 60%, 100% or 120%, and yellow spatial bins are common to both water levels

The RSD for the 120% moisture level appeared extremely different from the other moisture levels (Table 2, Figure 8). The roots did not grow very deep inside the rhizoboxes, but tended to spread out laterally to a greater extent than in the other soil water conditions. Our STARSEED approach captured this distinct difference between the root images, as most changes between the RSDs are shown through the arrows pointing down from the lateral roots of the 120% moisture level representative RSD. The root system for 120% also tended to be denser, as evidenced by the amount of high color value bins that indicate a higher percentage of root pixels in Fig. 8a.

Fig. 8

RSD, EMD score and largest value changes between 120% (a) and 60%, 80% and 100% (b–d respectively) for Run 3 and 4. The color scale for Fig. 8b-d represents the percentage of root pixels in each spatial bin. e shows the EMD score overall. f–h represent the the top 20% value changes between 120% and the three other soil water levels; red spatial bins are unique to 120%, green spatial bins are unique to 60%, 80% or 100%, and yellow spatial bins are common to both water levels

Discussion

Avoidance vs tolerance in waterlogged soil conditions

One interesting observation is that the three features that were considered (% root pixels, fractal dimension and lacunarity) had very similar discriminating power between cultivars and between moisture levels (Tables 1 and 2). Both fractal dimension and lacunarity have been used successfully in RSD and root architecture analysis [18, 31]. However, the main differences between these previous studies and the current one is that in this work these features were not calculated across the entire image but for each spatial bin. By using these texture features to describe the individual regions of an image, we improved the RSD characterization and comparisons for each treatment. While we calculated the fractal dimension and lacunarity values for the whole images, we observed limited differences among treatments for the CH index (e.g., S32 $= 0.06$ and S38 $= 0.03$ for Runs 3 and 4 with fractal dimension). The Kruskal-Wallis/Dunnett’s test with Bonferroni’s adjustment showed that % root pixels appeared to be a better feature than fractal dimension and lacunarity, as % root pixels led to the highest CH score for images from Runs 3 and 4 across treatments. For Runs 1 and 2, fractal dimension was the highest ranked feature for cultivar, while % root pixels was better for water levels. Therefore, for future similar work, % root pixels, an easily calculated parameter, may be preferential to fractal dimension and lacunarity for distinguishing between sesame RSDs.

When looking at these results more specifically, we saw that without fertilizer, cultivar appears to be the dominant driver of RSD, and that soil water only plays a minimal role (Tables 1 and 2). This can be seen in the highest CH index score on average across all cultivars compared with the averages for water level (Table 1). S32 was the cultivar with the most distinct RSD, and in Fig. 4, we could hypothesize what made this RSD different from the others. Specifically, the arrows representing the 20% highest value changes show that the main differences between the S32 and S35 RSDs were located mostly in the lower (deeper) section of the root system (Fig. 4f), and that on the contrary the main differences between the S32 and S40 RSDs were situated in the top (shallower) section of their root systems.

The situation is reversed when fertilizer is added, and soil water level becomes the factor with the biggest impact on RSD, as evidenced by the much higher CH index value for the soil water level treatment. This high score is actually driven by a single treatment, 120% soil water level, which had a very different RSD compared to the other soil water treatment (Fig. 8). The roots appeared to not grow as deep and had an increased lateral growth compared to the RSD for the three other water levels. This can be seen through the location of the white arrows, all pointing at changes between the RSD at 120% and the three other RSDs. This very distinct RSD in flooded conditions when fertility is adequate in Runs 3 and 4 contrasts strongly with that of Run 1 and 2. Sesame is known to be highly sensitive to waterlogged soil [32], though this sensitivity does not always translate into early total root length differences [26]. Instead, we saw here that the plant seemed to activate a variety of morpho-anatomical responses to cope with flooding stress, echoing the literature [33]. Here, we can highlight two different morpho-anatomical strategies employed by the crop depending on the soil fertility status. Without adequate nutrient availability, the crop seemed to opt for a tolerance strategy to flooding, developing roots into the bottom part of the rhizobox where the soil was saturated with water (Fig. 7). We can assume that this apparent tolerance strategy is accompanied by other changes not captured in this study such as the formation of aerenchyma and changes in the enzymatic activities [33]. However, when there are enough nutrients in the soil, the plants seemed to adopt an avoidance strategy, and did not grow roots down into the waterlogged soil, but tended to proliferate laterally (Fig. 8). We note that these observations are constant across the cultivars; as a result, fertility may condition the response of early sesame RSDs to flooding stress.

Validity of the EMD method

Although of interest in estimating early root vigor and biomass, only considering TRL provides limited insight about RSD as TRL does not consider the spatial arrangements of the roots [26]. Increasing our knowledge and understanding of root development in response to environmental stresses is of capital importance for global agricultural production, and has been defined as the pillar of the Second Green Revolution [1]. Many methods are now available to observe RSD and root architecture, yet some of the more advanced techniques such as CT imaging remain expensive and impractical, underlining the need for improved trans-disciplinary phenotyping approaches that can capture and quantify the complexity of RSD and root system architecture [34,35,36].

Our proposed STARSEED provides a visual, explicit, and quantitative characterization of RSD, allowing for spatially precise comparisons between RSD. A valuable aspect of the analysis is the introduction of the EMD class score, which summarizes all the differences between two RSDs down to a single number. The EMD method, combined with the CH index, can be used to very quickly know which treatment led to the most distinct RSD, and the degree to which these RSDs relate to one another globally. The STARSEED method could also be used in combination with other RSD and root system architecture measurements such as branching, angle, hierarchy, and fine-root distribution to achieve a comprehensive mechanistic understanding of root system.

The actual spatial EMD result (i.e., Figs. 4 through 8) can then used to further elucidate specific differences between RSDs. To the best of the authors’ knowledge, this is the first time EMD has been used to characterize RSD. The method has been thoroughly validated and used in varied fields of science, including fluid mechanics [37], linguistics [38], or, more related to the present study, image classification and comparison [39]. The CH index used in STARSEED serves as an evaluation step for the method to not only quantify what we observe, but also to ensure the approach captures distinct features to effectively group RSD based on shared spatial characteristics, which can then be linked to the treatment structure of the experiment.

One of the most critical strengths of STARSEED is that the proposed method allows for the generation of an “average” RSD visual representation called the representative RSD image for that treatment, as shown in Figs. 4 through 8a through d. By partitioning this newly generated image into local bins and calculating a feature value for each of these spatial bins, we can generate a visually explicit average RSD that can be compared to other representatives RSD to study differences between genotypes or environmental conditions (or any other treatment given a different experiment). Recent studies have attempted to either develop new techniques to directly measure RSD and/or root system architecture [3], find a way to accurately represent average root systems corresponding to specific growing conditions [40], and many more studies have developed or refined root development models [41,42,43,44]. One study in particular created similar 2D heat maps of “root frequency” observed on the four transparent surfaces of rhizoboxes but did not perform a spatially explicit analysis to the level of STARSEED [45]. Furthermore, the STARSEED approach is able to characterize and compare 2D RSD images of root systems at various scales. RSD is typically measured through soil coring at various depths, using a grid on a profile wall, or with mini-rhizotron systems [46,47,48]. These techniques lack an efficient way to integrate local data at the whole root system’s scale across all depths. STARSEED provides a way to get information about the entire RSD (with the outputs presented in Figs. 4 through 8 a, b, c and d) and at the same time highlight the largest RSD differences between two RSD signatures (Figs. 4 through 8f–h). This feature is quite unique and could be of great value to breeding programs where RSD need to be compared. Some of the shortcomings of STARSEED include the need to hand-label root images with WinRhizoTron. This is an extremely time consuming task that represents a serious limit to the rapidity and ease with which our approach can be used. However, recent work in machine learning is paving the way for fully automated and robust root tracing algorithms [49], which would dramatically increase the applicability of STARSEED. In addition, STARSEED remains an indirect measure of similarity between root systems as the proposed method is based solely on image analysis. STARSEED only indicates relative similarity or differences between two root systems captured with the same imaging system.

Going beyond the 2D images of rhizoboxes

Given the rhizobox set up of this study and the use of flatbed scanners for root images acquisition, the observations are confined to 2D when root systems always exist in 3D spaces. Imaging and characterizing a whole root system in 3D in situ or even in vitro is usually low throughput and expensive, if not nearly impossible. Most scientists thus resort to using models in lieu of direct observations, as previously mentioned. One of the main difficulty of using these 3D root development models is correct parameterization, which is crucial to the accuracy and validity of a model’s prediction [50]. Recent work has shown that 2D measurements of root systems could be used to adequately inform the parameters for 3D models in wheat [51]. It is thus a reasonable hypothesis to say that the current STARSEED approach could be further refined and subsequently used to generate and interpret such 3D models. There is potential to go even further and directly use STARSEED, and more broadly EMD, on 3D data [52], although acquiring such data still remains challenging.

Combining the average EMD visual maps with new RSD and root system architecture characterization techniques [53] have the potential to boost our understanding of root system plasticity. These tools would allow for both global (i.e., EMD scores) and local (i.e., magnitude and direction of changes) comparisons between genotypes and environmental conditions. We can also quantify the distribution of arrows as this information is captured from the flow matrix. We could record the distribution of the magnitude and direction of the arrows (similar to histogram-based features such as histogram of oriented gradients [54] or edge histogram descriptors [55]) for further analysis. Another interesting direction would be to develop a way to connect the observations made with STARSEED to more mechanistic measurements.

Conclusion

In this paper, we presented STARSEED, an approach to compare root system distributions from 2D images. Qualitative and quantitative analysis demonstrate the effectiveness of the proposed method. STARSEED successfully incorporated spatial and texture information to describe root architectures in both global and local contexts. The method is explainable and provides clear connections to biological aspects of each root image. STARSEED also allows for aggregating individual architectures to assess and average responses for each environmental condition and genotype. Future work includes clustering based on root architecture (i.e., using the pairwise EMD matrix for relational clustering), comparative study of different measures (e.g., convex hull) for feature extraction, applying our general framework over time to characterize how RSA develops in various conditions, along with scaling up to 3D comparisons.