Uncovering natural variation in root system architecture and growth dynamics using a robotics-assisted phenomics platform

  1. Therese LaRue
  2. Heike Lindner
  3. Ankit Srinivas
  4. Moises Exposito-Alonso
  5. Guillaume Lobet
  6. José R Dinneny  Is a corresponding author
  1. Department of Biology, Stanford University, United States
  2. Department of Plant Biology, Carnegie Institution for Science, United States
  3. Institute of Plant Sciences, University of Bern, Switzerland
  4. UCLouvain, Faculty of Bioengineering, Belgium
  5. Forschungszentrum Jülich, Agrosphere Institute, Germany

Abstract

The plant kingdom contains a stunning array of complex morphologies easily observed above-ground, but more challenging to visualize below-ground. Understanding the magnitude of diversity in root distribution within the soil, termed root system architecture (RSA), is fundamental in determining how this trait contributes to species adaptation in local environments. Roots are the interface between the soil environment and the shoot system and therefore play a key role in anchorage, resource uptake, and stress resilience. Previously, we presented the GLO-Roots (Growth and Luminescence Observatory for Roots) system to study the RSA of soil-grown Arabidopsis thaliana plants from germination to maturity (Rellán-Álvarez et al., 2015). In this study, we present the automation of GLO-Roots using robotics and the development of image analysis pipelines in order to examine the temporal dynamic regulation of RSA and the broader natural variation of RSA in Arabidopsis, over time. These datasets describe the developmental dynamics of two independent panels of accessions and reveal highly complex and polygenic RSA traits that show significant correlation with climate variables of the accessions’ respective origins.

Editor's evaluation

The authors present an automated system for phenotyping root system architecture based on bioluminescent roots resulting from a constitutively expressed luciferase transgene (GLO-Roots). They have developed a robotics-assisted phenotyping platform and an automated image analysis pipeline for high throughput analysis. An impressive array of 93 luciferase expressing Arabidopsis thaliana accessions provides a major resource for understanding the genetic basis for root system architecture variation under physiologically relevant conditions. The work will be of great interest to plant biologists and all those studying genetic variation in plants.

https://doi.org/10.7554/eLife.76968.sa0

Introduction

The diversity of shoot and root system forms within the plant kingdom mirrors the many functions required for plant survival in diverse habitats. However, these forms are not static, but dynamically change during the lifespan of a plant as it acclimates to a changing environment, and over evolutionary time scales as a consequence of natural selection and genetic drift. The visible above-ground shoot system is optimized for photosynthesis and reproduction while the below-ground root system is responsible for anchoring the plant and allowing it to reach essential nutrients and water, which is mediated, in part, through interactions with microbes in the soil (Lundberg et al., 2012; Bulgarelli et al., 2012; Castrillo et al., 2017; Harbort et al., 2020). While shoot forms and structures have undergone centuries of scrutiny under the human eye, underground root systems have received much less attention due to impeded accessibility.

Root system architecture (RSA) refers to the spatial arrangement of roots in soil, which encompasses the geometric nature of root connectivity and distribution (Lynch, 1995; Rellán-Álvarez et al., 2016). This arrangement of roots depends on relationships among the growth rates, branching frequency, and gravity set point of different root types. In allorhizic root systems, common in Eudicotyledons like the model plant A. thaliana, the embryonically formed primary root is the first to emerge from the seed at germination and ultimately branches to give rise to secondary roots, which, in turn, branch and produce tertiary and higher order roots (Osmont et al., 2007; Fitter, 2017). Each step of forming a branch, from the priming of a lateral root primordia to emergence from the parent root, is a tightly controlled process regulated by endogenous and exogenous signals (Moreno-Risueno et al., 2010; Péret et al., 2009; Van Norman et al., 2013; Motte et al., 2019). Modulation and crosstalk of plant hormones change local growth patterns. The ratios of the phytohormone auxin and its antagonist, cytokinin, control growth throughout the plant and, importantly, regulate primary and lateral root initiation (Tian et al., 2014; Morris et al., 2017). The hormone abscisic acid (ABA), on the other hand, inhibits lateral root primordia from emerging (De Smet et al., 2003). ABA levels can increase in response to external stress and cause an extended quiescent phase after lateral root emergence (Duan et al., 2013; Xiong et al., 2006). Hence, the precise balance of these hormones influences root development and shapes overall RSA.

In addition to their genetically determined growth program, plants also react to spatial and environmental cues to tune their form for certain environments. Roots must strategically navigate through the soil environment for efficient acquisition of water and nutrients (Rellán-Álvarez et al., 2016; Dinneny, 2019). Understanding the mechanisms by which roots sense and respond to their environment remains a significant challenge since existing methods either compromise on physiological relevance or on throughput, which can impact the ability to deploy functional genomic tools to identify regulatory genes. Many RSA studies are performed on young root systems using easily accessible, in vitro gelmedia-based systems and leverage the powerful genetic and molecular resources of A. thaliana (Rosas et al., 2013; Ogura et al., 2019; Waidmann et al., 2019; Julkowska et al., 2014). Whereas bigger and more mature root systems of crop species have been studied in soil environments, however, compromising on throughput (Morris et al., 2017; Jiang et al., 2010).

Previously, we developed GLO-Roots (Growth and Luminescence Observatory for Roots), an imaging platform enabling visualization of soil-grown Arabidopsis root systems from germination to maturation by combining custom growth vessels, luminescence-based reporters, an imaging system, and an image analysis suite (Rellán-Álvarez et al., 2015). Here, we present the automation of GLO-Roots by creating GLO-Bot, a robotic platform that enables large-scale capture of Arabidopsis root growth over time in a soil-like environment. Along with the automated imaging system, we developed an improved root analysis pipeline for quantification of various RSA traits. Together, these technological advancements enabled time-lapse imaging of RSA development of 93 Arabidopsis accessions from 14 to 28 days after sowing, creating a unique dataset that allows for wide-ranging analyses of highly complex, adult root systems. Quantitative genetic analyses and genome-wide association studies (GWAS) yielded insight into the heritability and genetic architecture of distinct RSA traits, and significant associations with climates across the species’ geographic distribution highlighted the power of our approach and its relevance for understanding how root architectural traits contribute to local adaptation.

Results

Automation enables high spatio-temporal resolution measurements of root system dynamics

Establishing the robotics infrastructure for automated root imaging

To further develop the potential of GLO-Roots for time-lapse imaging of root development, from early stages to mature RSAs, we developed GLO-Bot, a Cartesian gantry robotic system consisting of an exterior frame and arm on a linear rail system built from T-Slotted aluminum extrusion, which houses the Growth and Luminescence Observatory 1 (GLO1) imaging system (Rellán-Álvarez et al., 2015; Figure 1A and B). In addition to GLO1, the frame contains a growth area for seven black polyethylene growth boxes (Figure 1B), each of which holds 12 rhizotrons (also known as rhizoboxes) arranged in a two by six grid (Figure 1C, Figure 1—figure supplement 1). In total the system can hold 84 rhizotrons, which typically contain one plant.

Figure 1 with 1 supplement see all
Growth and Luminescence Observatory for Roots (GLO-Roots) automation: GLO-Bot.

(A) Schematic of the GLO-Bot Cartesian gantry system includes: (1) The Growth and Luminescence Observatory 1 (GLO1) imaging system published in Rellán-Álvarez et al., 2015, which houses two cameras and a rotating stage for root imaging, (2) a station for general watering or treatment with diluted luciferin solution prior to imaging, (3) a barcode scanner to identify the rhizotron and load a specific watering and imaging protocol, (4) a robotic arm, which moves in the x-, y-, and z-directions and has a hook at the end to pick up rhizotrons, and (5) an area for plant growth, which can be seen in the photograph of GLO-Bot (B). (C) Automation updates required modification of the GLO-Roots growth vessel design to include a black acrylic plate and hooks for rhizotron handling as well as dividers within the growth boxes and guides along the bin top, which allow the rhizotron to hang and shield the roots from light. Copper tape along the edge of the bin top enables positioning. Gray-scale bars denote 30 cm.

The primary goal of developing the robotics platform was to automate plant imaging since the luciferin watering regime of each rhizotron and the long exposure time required to capture luminescence signal were the most time-consuming steps of the GLO-Roots system. An important consideration was the layout for the robot and range of motion for the arm (Figure 1A). The perpendicular position of the two cameras relative to the rhizotron being imaged in GLO1 lent itself well to the final linear rail system and arm with x-, y-, and z-motors. To minimize positional errors associated with finding a particular rhizotron, the robotic arm checks the coordinates of each growth container prior to imaging (Video 1). The arm moves to each plant position and lifts the rhizotron via a metal hook, which gives the arm a robust surface to pick up the rhizotron (Figure 1C, Figure 1—figure supplement 1A). Both rhizotron and growth box designs were adapted for robotic use in a way that prevents light exposure of the root system during growth (Shi et al., 2018) while also facilitating easy removal and repositioning of rhizotrons in the growth container during imaging. Each rhizotron is shielded from light exposure by creating individual chambers within the growing box (Figure 1C, Figure 1—figure supplement 1B). To facilitate removal and replacement of each rhizotron for imaging, we changed the design of the growth box from a standing to a hanging position by adding a piece of acrylic at the top of each rhizotron, which, with the help of small guide pegs, allows the rhizotron to hang in its proper position and shields light from entering the growth container (Figure 1C).

Video 1
GLO-Bot running.

Programed movements of the robotic arm bring the rhizotron to a barcode scanner, which allows the system to recall the specific treatment and imaging protocol to use (Video 1). Subsequently, the arm brings the rhizotron to a watering station where a luciferin solution is applied prior to imaging, which provides the substrate for the constitutively expressed luciferase (Rellán-Álvarez et al., 2015; Figure 1A; Figure 2). Slow back and forth movements of the rhizotron under a pipette tip connected to a peristaltic pump provide the solution ( Video 1). As the arm moves, it bypasses a specified gap in the center of the rhizotron to avoid contact with the shoot. The entire watering process takes approximately 10 min, which allows the water to slowly seep through and spread throughout the rhizotron. The automation was set up such that while one plant is being imaged, the next plant is being watered. To minimize arm movement time, watered rhizotrons are placed on a ledge while the arm removes the imaged rhizotron from the imager and returns it to its growth box (Video 1). If the system is running at maximum capacity, GLO-Bot can process 96 rhizotrons in 24 hr, by growing an eighth box on a nearby shelf and swapping out one of the seven boxes in the system during the day. However, the following experiments were all imaged at approximately half capacity in order to respond to possible technical errors.

Workflow of image analysis pipeline that enables robust trait measurement.

(A) Sample workflow for processing time-series images from multiple days after sowing (das) starts by (1) combining raw images to merge front, back, upper, and lower images, (2) a blank image is then subtracted to remove inherent noise, (3) images are registered using the ImageJ plugin Template Matching and Slice Alignment to account for x- and y-movement caused by slight position changes during rhizotron insertion into the imager, (4) registered images are then de-noised and overlaid, which helps overcome luciferase signal loss in older root tissues; however, this step also compounds any background noise. (5) Image subtraction between each day removes this noise and isolates new growth. Processed images are run through GLO-RIAv2 (6) and the output is formatted, such that in silico vectorized roots can be reconstructed (7). Additional noise is cleaned up (8) using an iterative distance-clustering method. Traits are then extracted from cleaned roots (9) and the outputs can be used for downstream analysis (10). (B) Comparison of 85 images measured with previous root analysis methods (GLO-RIAv1, gray) and the new analysis methods (GLO-RIAv2, coral) vs. manually traced ground truth measurements from SmartRoot (x-axis) reveals that the new analysis method increases accuracy, as demonstrated by the r2 values, higher slopes, and lower y-intercepts.

Development of an image analysis pipeline for dynamic trait quantification

Image alignment, signal optimization, and background noise removal

The increased throughput of GLO-Bot enabled time-lapse imaging of root systems and created new opportunities for the measurement of dynamic root growth traits, but also revealed new challenges for root image analysis. Luciferase signal intensity varies between transgenic LUC-expressing lines and decreases in older root tissues, making it even more difficult to distinguish between roots and noise during the course of the experiment (Figure 2). While this signal variation is problematic for analysis of raw images, we leveraged our time-lapse data to capture a strong signal for the entire root system by sequentially adding images together. However, in addition to summing the root signal, this step also sums background noise (Figures 2 and 3). Therefore, we utilized a combination of pre- and post-measurement processing to amplify true signal and reduce noise (Figure 2). To account for some of the inherent noise introduced by the imaging system itself, we run the background removal macro remove_background.ijm, to subtract a blank image, taken in GLO1 without a rhizotron present, from every root image (Figure 2A-2). True root signal can be amplified, or inferred from previous days in the time-lapse, if all images of a root system have the root in the exact same position relative to the exterior image dimensions. Therefore, each series of root images is manually checked to ensure root position remains stationary over time. Misalignments were fixed using the ImageJ plugin Template Matching and Slice Alignment (Template_Matching.jar, Tseng et al., 2011). This plugin works by user selection of a region present in all images followed by automated detection and translation (x-, y-movement) of subsequent images to align image features, a method that is also known as image registration. This step occasionally required subsetting the images to yield precise alignment. During the alignment process, all root images were screened and roots were checked for: (1) appropriate number of images; (2) erroneous image dates; (3) excessive background noise; (4) any other anomalies from image testing (described in the Materials and methods ‘Image analysis’ section) or occasional camera malfunction. An incorrect number of images or erroneous image dates were due to technical issues with the automated plant handling and imaging pipeline. Images with excessive background noise, or extraneous images, were removed. Next, aligned and cleaned-up images were run through the macro denoise_overlay.ijm, designed to strengthen signal intensity and decrease noise through a series of image-wide addition and multiplication steps (see detailed steps in the Materials and methods ‘Image analysis’ section). The output from this macro is a set of images in which each day of imaging is added to the previous days, thus accumulating the strongest signal throughout the root system (Figures 2A4).

Figure 3 with 3 supplements see all
Root system architecture (RSA) of six Arabidopsis accessions over time.

(A) Six accessions from diverse locations were imaged continuously from 9 to 31 days after sowing (DAS). *Col-0 location reflects where accession was originated, not collected. (B) Each accession had a unique RSA pattern (five root systems per accession at 21 DAS each shown in a different color overlaid on top of each other). (C) Principal component analysis (PCA) using raw trait values for eight time points and nine traits shows time (red is 9 DAS and yellow is 31 DAS), as the largest source of variation (PC1) and angle and depth-to-width ratios further distinguishing the root systems (PC2). When colored by accession, this analysis shows separation between the accessions. Loading plot illustrates the impact of each trait on the overall variation. (D) Four root system raw traits for six accessions growing over time from 9 to 31 DAS. Trends were calculated using LOESS smoothing (n=12–15 plants at each time point). 95% confidence interval shown by gray shading. Non-smoothed raw data and additional traits are plotted in Figure 3—figure supplement 2. Two vertical dashed lines highlight the interval (14–28 DAS) chosen for future experiments. Depth and width measurements are constrained by the rhizotron size. Colors correspond to accessions in panel B.

Figure 4 with 3 supplements see all
Root system architecture (RSA) of the Growth and Luminescence Observatory for Roots (GLO-Roots) diversity panel.

(A) Locations of all accessions used in this study. Inset highlights the concentration of accessions from the Swedish population. (B) RSAs of the GLO-Roots diversity panel used in the genome-wide association studies (GWAS) at 20 days after sowing (DAS). Six root systems overlaid on top of each other with root system color indicating each replicate. Root systems are arranged in order of median average angle of the root system from deepest in the upper left to shallowest in the bottom right.

With luminescent root signal maximized throughout all of the images in a time-lapse series, the macro invert.ijm is used to invert the gray-scale values of the images, thus making the root system itself black and the background white for downstream whole root system analyses. In addition to determining architectural parameters of complete root systems, we also established methods to quantify the actively growing regions of root systems. For this, we use the macro tip_tracking.ijm to isolate the new growth between images by a series of dilating and subtracting of two consecutive images (see detailed steps in the Materials and methods ‘Image analysis’ section; Figures 2A5).

Figure 5 with 1 supplement see all
Trait relationships among the Growth and Luminescence Observatory for Roots (GLO-Roots) diversity panel.

(A) Principal component analysis (PCA) of nine fitted traits (i.e. fitted values) for the 93 accessions of the GLO-Roots diversity panel reveals traits influenced by time as the predominant source of variation within the dataset (PC1), while angle measurements and the depth-to-width ratio further distinguish accessions (PC2). (B) A correlogram of the fitted trait values at each time point, increasing left to right and top to bottom, further shows the negative correlations between the depth-to-width ratio and the angle traits and displays the relationships between all nine traits at each time point. Correlations where p>0.01 have been crossed out with an ‘X’. While some traits remain stable over time, other traits change over time. The extent to which a trait changes varies between accessions. (C) In some accessions, the fitted angle of new growth (average angle per day) remains constant each day, while in others the angle of new growth increases or decreases over time. The accessions are ordered based on variability of average angle per day with similar angles to the left and strikingly different angles on the right. (D) Ale-Stenar-64–24 and San Martin demonstrate contrasting root architectures, with opposite depth-to-width ratios and opposite average angles. (E) TDr-18 and Ting-1, on the other hand, are examples of accessions with constant and variable (respectively) average angle per day measurements.

Downstream analyses of root growth dynamics requires the distinction between true root signal and background pixels

All images were run through a modified version of GLO-RIA (Growth and Luminescence Observatory Root Image Analysis; Rellán-Álvarez et al., 2015), GLO-RIAv2 (Figures 2A6). GLO-RIAv2 removes features deprecated from the 2017 ImageJ update (Schneider et al., 2012; Schindelin et al., 2012) and includes root angle output ranging from 0° to 180°.

Figure 6 with 1 supplement see all
Genome-wide association studies (GWAS) and climate correlations of root traits of the Growth and Luminescence Observatory for Roots (GLO-Roots) diversity panel.

(A) Manhattan plot for average angle per day at 28 days after sowing (DAS). Black points indicate SNP positions that pass the Bonferroni threshold (solid line) at least once throughout the time series, as shown in the inset plots. Gray points are those in linkage disequilibrium (LD) with the significant SNPs. (B) Correlation between the average angle per day and the mean diurnal temperature range indicates the potential importance of root angle for surviving highly variable climates. (C) A world map showing mean diurnal range (black represents larger fluctuations) with the distribution of accessions (points). Points in (B) and (C) colored by cluster identity based on changes in average angle per day over time and computed by between-group average linkage hierarchical clustering with distance = 1.75, as shown in the heatmap (D). Heatmap depicts average angle per day for each day with blue indicating steeper angles and yellow indicating shallower angles. Changes in average angle per day divided the accessions into three clusters: the yellow cluster with consistently shallow root growth, the blue cluster with consistently deep root growth, and the red cluster with intermediate or changing root angles.

Figure 6—source data 1

Manhattan plots and SNP development of all nine root system architecture (RSA) traits.

(A) Manhattan plots of genome-wide association (GWA) run on each of the nine RSA traits at each time point with corresponding Q-Q plots (B). SNPs highlighted in black pass the threshold, gray indicates SNPs in linkage disequilibrium (LD). Genes in this region are reported in Supplementary file 4. (C) SNP development of the 29 unique SNPs that pass the Bonferroni threshold at least once throughout the time series.

https://cdn.elifesciences.org/articles/76968/elife-76968-fig6-data1-v2.pdf

Downstream analyses of the GLO-RIAv2 outputs were performed in R (R Development Core Team, 2019). Most root system structures visible on the GLO-Bot images are disconnected in many locations, due to variability in the reporter’s intensity and the presence of physical obstructions, such as soil particles. We call each detected bioluminescent piece of root a ‘root segment’. Data for each root segment include (1) the x-, y-coordinate of the upper left corner of the bounding box for each segment, (2) the segment length, and (3) the angle with respect to gravity, such that root segments can be referenced with 0° pointing straight down and 90° pointing straight out to the side. These data give us the ability to reconstruct the root system in silico as a series of vectors (Figure 2).

Next, we remove non-root particles by iteratively clustering the starting points of the measured segments. Clusters are tested for distance between points, the number of segments, and the linearity of the points within the cluster (detailed steps in the Materials and methods ‘Image analysis’ section, Figure 2A-8). Segments that do not meet empirically defined criteria are removed, leaving behind only ‘true root’ segments. Using the ‘true root’, we calculate traits associated with new growth including: number, length, and angle of new root segments, as well as traits at the whole root system level: width, depth, center of mass, convexhull, depth-to-width ratio, total length, and average angle (Figure 2A-9). When computing average angle measurements, vector segments are weighted by length, meaning that small segments contribute less to the overall trait value compared to longer root segments. In addition to these seven traits, new root length (length per day) and average angle (average angle per day) are also summarized on a per-day basis. Additional traits can be derived from these metrics and extracted from the raw images for future analysis. It is important to note that each computed trait must be considered in the context of the experimental system. For example: (1) width and depth measurements are ultimately constrained by rhizotron size making these traits relevant up to a certain time point depending on the accession, (2) correspondingly, the depth-to-width ratio will always converge close to 2:1 as the root system fills the rhizotron, and (3) root length is an approximation since it is a sum of the visible root segment lengths and therefore does not account for the distance between disconnected segments. Despite these caveats, the nine traits measured (listed above and summarized in Supplementary file 1) encompass the basic parameters commonly used to describe the root system as well as new traits that are made possible through our time-lapse imaging and automated plant handling system.

To test the accuracy of our trait extraction methods, we used SmartRoot (Lobet et al., 2011) to manually trace the root systems for five growing plants and ground truth our parameters. We observed a strong correlation between the traits extracted via our new image processing pipeline and those that were measured manually (Figure 2B). Manual SmartRoot tracing of a single GLO-Roots image can take multiple days, thus making analysis of large-scale time course experiments unfeasible. While the initial GLO-Roots pipeline works well for end-point image analysis, the above-mentioned updates allow for more accurate root trait measurements of time course images (Figure 2B) and gives us ways to explore new temporal traits.

Time-lapse imaging of six Arabidopsis accessions reveals distinct root growth trajectories

The robotics and image analysis pipeline we developed allowed us to increase the throughput of the GLO-Roots system, which we then applied to address two important questions: How do root system traits vary over developmental time, and how do dynamic root system traits vary across a species? To address the first question, we wanted to characterize root system traits with high temporal resolution and with greater replication than possible with a larger survey of genotypes. We selected a set of six Arabidopsis accessions based on their geographic dissimilarity and availability of homozygous LUC2-expressing lines with strong reporter expression (Figure 3A; Rellán-Álvarez et al., 2015). Each accession was grown in 10 replicates, where 5 replicates were imaged daily from 9 days after sowing (DAS) until 31 DAS and the other 5 replicates from 21 to 31 DAS. In depth comparisons between replicates can be done using the RShiny App (https://tslarue.shinyapps.io/rsa-6access-app/). The six accessions that were imaged displayed diverse RSAs (Figure 3B, Video 2) and principal component analysis (PCA) demonstrated that we could distinguish the accessions based on nine of the traits extracted using our image analysis pipeline (Figure 3C, Supplementary file 1). Expectedly, time accounted for the largest source of variation in PC1 (56.99%, Figure 3C), indicating that the size and developmental stage of a plant is the most important variable affecting RSA. We could see that depth-to-width ratios and average angle measurements further separated the accessions along PC2 (21.67%, Figure 3C). In total, the first four PCs capture 94.02% of the variation with the greatest distinction between accessions apparent in PC2 (Figure 3—figure supplement 1). Looking at the traits over time demonstrated that each of the accessions had its own unique growth trajectory with average angle and average angle per day exhibiting the most diversity (Figure 3D, Figure 3—figure supplement 2). Early during plant growth all accessions have fairly similar average angles, presumably because they have few lateral roots, but over time their average angles for the whole root system shift and distinguish the accessions. By looking at these data, we established that most traits stabilized by the end of the 28 DAS growth period (Figure 3D, Figure 3—figure supplement 2). Importantly, having imaged the plants both continuously and only at 21 DAS, we were able to confirm that these traits were not influenced by continuous imaging and the associated luciferin watering (Figure 3—figure supplement 3). Time-lapse imaging of these six accessions proved the robustness of GLO-Bot running for 15 hr a day over several weeks, the ability of our image analysis pipeline to determine various root system traits over time, and that even a small subset of accessions showed substantial variation in RSA traits. In addition, the application of our newly developed imaging protocols to these data allowed us to optimize experimental parameters, including luciferin watering speed and frequency, imaging frequency, and imaging time period. Those parameters were foundational to examine a larger, more variable population.

Video 2
Root system architecture (RSA) growth of six Arabidopsis accessions over time.

Natural variation in dynamic RSA traits is highly complex and polygenic

The GLO-Roots diversity panel

To explore the diversity of root architectures in the Arabidopsis species and define the genetic basis of these traits, we sought to apply the GLO-Roots method to a larger panel of curated accessions. We therefore transformed an initial set of 192 A. thaliana accessions with the ProUBQ10:LUC2o reporter for use with the GLO-Roots system (Rellán-Álvarez et al., 2015; Figure 4—figure supplement 1A, B). These accessions comprise the extended GLO-Roots diversity panel, which originated from a densely sampled Swedish population (Long et al., 2013) and were supplemented with accessions from extreme environments with high and low annual average precipitation (Fick and Hijmans, 2017), as well as accessions with high and low sodium accumulation levels in the shoot (Baxter et al., 2010). Ultimately, the GLO-Roots diversity panel, a set of 93 accessions, were selected for the GWA population (Figure 4ASupplementary file 2). This set balances genetic similarity and environmental variation by combining 72 Swedish accessions with 20 accessions from extreme environments. These lines were selected based on location, reporter expression strength, and the number of independent transformation events recovered to avoid insertion effects. Four Col-0 control plants grew in different growth boxes to examine position effects (Figure 4—figure supplement 2). This population not only represents a foundational resource for describing phenotypic trait diversity in the Arabidopsis species, but can later be used to identify genes underlying this trait diversity over time, using the tool of GWA mapping.

Although GLO-Bot could image 96 rhizotrons in 24 hr, we decided to run it at half capacity in order to have time to respond and adjust to technical errors as they arise. To capture the breadth of RSA dynamics within the diversity panel, we imaged 48 rhizotrons on 1 day (Figure 4—figure supplement 2, part 1) and the other 48 rhizotrons on the next day (Figure 4—figure supplement 2, part 2) from 14 to 28 DAS (Figure 4B, Video 3, Figure 4—figure supplement 1C). We grew all 93 accessions at once and replicated the experiment six times. Next, using our image analysis pipeline, we focused on the nine root traits that distinctly separated RSAs in the subset of six accessions, and which encompass RSA at a broad scale (Supplementary file 1).

Video 3
Root system architecture (RSA) growth of Growth and Luminescence Observatory for Roots (GLO-Roots) diversity panel over time.

Growth of root systems in soil is likely to incur an experimental cost due to the greater variability in experimental conditions, compared to gel-based media, and the resultant increase in variability of the root traits measured. Furthermore, the size of the diversity panel we used required that biological replication be performed on separate dates, due to limits in the capacity of the automation pipeline. To account for these additional sources of variation, we developed an analysis pipeline that allowed us to summarize trait data for our replicated time course experiments. For each measured trait, we estimated ‘fitted values’ to account for replicate noise (Figure 4—figure supplement 3). Fitted values were calculated using a generalized linear mixed model for each trait using the R package MCMCglmm (additional information can be found in the Materials and methods ‘Trait processing’ section; Wilson et al., 2010; Mrode, 2014; Hadfield, 2010). The model provided an intercept and slope for each genotype which were used to calculate the expected values at 0, 48, 96, 144, 192, 240, 288, and 336 hr after the start of imaging (equivalent to 14, 16, 18, 20, 22, 24, and 28 DAS). These fitted values were used for further analyses.

Root trait diversity and phenotypic relationships among the GLO-Roots diversity panel

Similar to the analysis of the first set of six accessions (Figure 3C), PCA of the GLO-Roots diversity panel revealed that the first principal component captures trait variation over time (Figure 5A, Figure 5—figure supplement 1A), as well as traits that are strongly associated with time and, thus, developmental stage and plant size (Figure 5A, Figure 5—figure supplement 1B) such as the root system’s area (convexhull area), length (total length), and depth. The second principal component demonstrates that both the average angle of the root system and the average angle of new growth account for a large amount of root system variation (Figure 5A, Figure 5—figure supplement 1C, D). Average root angle and depth-to-width ratio are negatively correlated (Figure 5A, B), as demonstrated by the accession San-Martin, which has the smallest depth-to-width ratio and the highest average root angle (i.e. most shallow) as well as Ale-Stenar-64-24, which has one of the highest depth-to-width ratios and lowest average root angles (Figure 5D). Additional data exploration and accession comparisons can be done using the RShiny App (https://tslarue.shinyapps.io/rsa-app/), which plots the vectorized root systems for all of the accessions. The raw data is available through Zenodo at https://doi.org/10.5281/zenodo.5709009. Accessions and replicate numbers with excessive background noise or other anomalies are identified in Supplementary file 3.

Further analysis of the traits via pairwise comparison at all time points not only allowed us to visualize existing correlations within traits over time, but also among traits (Figure 5B). Looking at average angle per day, we see that consecutive days show strong correlation across accessions (e.g. day 14 compared to day 16: Spearman’s rank correlation 0.99), whereas later time points are less well correlated with early time points (day 14 compared to day 28: Spearman’s rank correlation 0.56, Figure 5B). This suggests that the average angle of lateral root tips is a trait that changes over time with some accessions showing greater change than others (Figure 5C). Indeed, examining this trait across our 93 accessions reveals great diversity in how stable lateral root tip angle is across development (Figure 5C). While many accessions have a consistent angle of new growth, such as TDr-18 and San2 (Figure 5C, Figure 5E), others change their angle of new growth over time, such as TAA14 and Ting-1 (Figure 5C, Figure 5E).

RSA traits are heritable and highly polygenic

We next wanted to investigate the genetic basis of RSA traits. Some traits, such as the angle measurements, had high broad-sense heritability while others, such as total length, had low broad-sense heritability (Supplementary file 5). While some of the low heritabilities could be due to the limitations of the experimental method of growth or trait measurement, these data indicate measurable genetic control over some traits, which led us to use GWA to identify potential causal loci. Using the fitted values, we conducted GWA analyses for all nine traits at each time point using the AraGWAS database SNP matrix and the program GEMMA (Zhou and Stephens, 2012; Togninalli et al., 2018). In total, 29 SNPs were significant for at least one time point using Bonferroni multiple-test p-value correction (Supplementary file 4). These alleles also had a minimum allele frequency greater than 0.05. These SNPs highlight regions of the genome with elevated significance and the temporal nature of our data allows us to track the changes in significance over time (Figure 6A, Figure 6—source data 1). For example, when looking at average angle over time, we see a greater strength of association, as indicated by lower p-values over time, which likely corresponds to the stabilization of the trait (Figure 6A, Figure 6—source data 1B). In contrast, we see the strength of the association decrease for depth-to-width ratios over time (Figure 6—source data 1B), which likely reflects that growth of root systems becomes constrained at later time points by the physical limits of the rhizotron. We gain confidence in the chromosomal regions identified by tracking the significance over time and observing the slight enrichment of experimental p-values in our Q-Q plot distributions (Figure 6—source data 1).

Out of the 29 identified significant SNPs, four synonymous SNPs can be found within the coding sequence of two genes, eight were in intergenic regions, whereas the majority of SNPs were found in upstream- and downstream regions of genes, potentially affecting expression patterns (Supplementary file 4). It should be noted that the most significant SNP is not always the causal SNP. For this reason, Supplementary file 4 lists all of the genes within a 20 kb region surrounding the SNP (10 kb upstream and 10 kb downstream). The identification of many significant SNP associations suggests that RSA development is a polygenic trait.

Looking at the two traits with the highest heritability (Supplementary file 5), average angle (H2=0.45) and average angle per day (H2=0.31), we found 12 and 8 significant SNPs, respectively (Supplementary file 4). Interestingly, both traits revealed one common SNP upstream of a hypothetical protein (AT3G42390) with 30% protein identity to a phospholipase D alpha 3 (At5g25370)(Supplementary file 4). Furthermore, the only intragenic SNP for average angle per day was found in an intron of another hypothetical protein with 35% identity to a non-specific phospholipase C2. Phospholipases are lipid-hydrolyzing enzymes that are known to play a role in signaling during plant development, stress responses, and responses to environmental cues (Takáč et al., 2019). We identified SNPs in two genes encoding proteins that show similarity to phospholipases suggesting a potential role of phospholipids in adjusting root angles to local environments.

Root depth was associated with a SNP upstream of a gene encoding Cyclin A2;3 (AT1G15570; Supplementary file 4). This gene was previously found to play a role in auxin-dependent mitotic-to-endocycle transition that is involved in the transition from cell proliferation to cell differentiation in the Arabidopsis root meristem (Ishida et al., 2010). Thus, alterations in gene expression could influence root depth in certain Arabidopsis accessions.

The GWAS on the depth-to-width ratio yielded an intronic SNP in a P-loop nucleoside triphosphate hydrolases superfamily protein with Calponin Homology domain-containing protein (Supplementary file 4). This gene is expressed in the lateral root cap (http://bar.utoronto.ca/eplant/; Waese et al., 2017) and changes in expression may influence root gravitropism leading to a more shallow root system. Future work will establish the molecular basis for the SNP-phenotype relationships identified here.

RSA traits of the GLO-Roots diversity panel significantly correlate with climatic variables

To gain insight into the relevance of the different root architectures of the GLO-Roots diversity panel in the natural environment, we correlated our nine root traits with published bioclimatic variables (Fick and Hijmans, 2017). We see the strongest correlations between these traits and climatic variables related to temperature (Figure 6—figure supplement 1). Specifically, average angle per day shows a notable correlation to mean diurnal temperature range (Pearson r=0.373, p=2.56 × 10–4; Figure 6B). This result indicates that accessions with shallow RSAs tend to grow in climates that exhibit larger changes in the daily minimum and maximum temperatures. Likewise, convexhull shows a positive correlation to mean diurnal range, albeit less significant (Figure 6—figure supplement 1). These data suggest that plants with shallow roots may actually grow better than those with deep roots in a more variable environment. Clustering of accessions by change in average angle per day over time divides the accessions into three groups: consistently shallow, consistently deep, and intermediate or changing root angle (Figure 6C, D). Interestingly the consistently shallow and consistently deep accessions segregate from each other in a global map, which is consistent with our finding that climate is likely to be an important variable determining the distribution of accessions with these traits (Figure 6C).

Discussion

Implementation of automation expands the number of individuals that can be phenotyped and the phenotypes that can be captured (Gehan et al., 2017). By automating the GLO-Roots system, we created GLO-Bot, a robotic platform that enables unprecedented insight into Arabidopsis root growth over time in soil and at developmental stages rarely observed. Careful design of the system to facilitate automated handling allowed us to maintain physiological relevance while increasing throughput and allowing for novel trait measurements.

Along with the automated imaging system, we developed an improved image analysis pipeline, GLO-RIAv2, for quantification of RSA. The previous image analysis program was semi-automated, and often needed intervention to define regions of interest. While this is feasible for end-point imaging with a small number of samples, it was not scalable with GLO-Bot. For automated image analysis, however, one of the biggest challenges was the combination of maximizing root luminescence signal intensity while decreasing background noise so that GLO-RIAv2 would only analyze true root signal. We saw a strong correlation between our traits extracted via our new image processing pipeline and those that were measured manually, demonstrating that our updates allowed for accurate examination of growth dynamics in both new growth and at a whole root system level. Recent progress in machine learning is demonstrating the opportunities within the field for segmenting out root and noise and could further improve accuracy (Wang et al., 2019; Smith et al., 2018). Whereas RGB imaging together with software tools like RootPainter has facilitated automated root system analysis in species with thicker roots (Smith et al., 2020; Bauer et al., 2022), bioluminescence still has advantages with species with finer roots, like Arabidopsis, since it improves the contrast between the root and soil signal. Using transgenic lines homozygous for the luciferase transgene would lead to stronger signal and therefore better root detection. New advances in recapitulating a full bioluminescence pathway in planta (Khakhar et al., 2020) could provide an exciting alternative for increasing root luminescence signal and therefore simplifying image segmentation.

The automation of GLO-Roots allowed for imaging larger sample numbers, which enabled the visualization of root growth of different Arabidopsis accessions over time. At maximum capacity, GLO-Bot could image 96 rhizotrons within 24 hr. However, all experiments were run at approximately half capacity in order to facilitate the management of handling errors by the system. The biggest source of errors occurred during the return of a rhizotron back to its growth box, which we addressed by adapting the design of the holding bins. Further advances in robotics and in the precision with which components of the system are manufactured will likely improve the throughput further.

Our automated plant handling and image analysis pipeline allowed us to explore the importance of time and genetic diversity in driving changes in RSA. A set of six accessions that were homozygous for the luciferase reporter allowed us to image root systems every day with a high degree of replication. Such experiments also tested the robustness of the system and facilitated optimization of imaging parameters and watering times to extract root traits that reflect true differences in RSA. Nine root traits were automatically extracted through GLO-RIAv2 and showed striking differences even in this small subset of accessions. Our analysis of these data revealed that differences between accessions often only emerged over time, such as average root angle, which highlights the importance of quantifying root traits over time and at these later developmental stages.

The GLO-Roots diversity panel provided access to a greater amount of genetic variation, but required changes in experimental design and data analysis methods. The GWAS ultimately led to the identification of significant SNPs associated with variation in six of nine root traits measured. However, given the polygenic nature of these traits, our ability to detect individual causal SNPs or peaks may be low. Additionally, our relatively small population size (93 accessions) challenges our ability to discover SNPs significantly associated with the traits (Gibson, 2012). Although a surprisingly rapid linkage decay in A. thaliana (Nordborg et al., 2005) allows within-gene trait mapping (Atwell et al., 2010; Exposito-Alonso et al., 2018), it will be ideal to follow up the SNPs identified above and other close variants through further analyses, such as accession-specific gene expression and complementation, knock-out and reporter gene assays (Ogura and Busch, 2015). Besides the analyses conducted in this study, the generated datasets present valuable resources for further studies related to the genetic analysis of the principal components to explore the phenotype space, multi-trait analyses, modeling of phenotypic variation across time, and structure-function modeling between RSA and water use.

The correlation of our nine extracted root traits to published bioclimatic data from the origin sites of the respective accessions reveals a significant correlation between average angle to mean diurnal temperature range. Thus, accessions with shallow RSAs tend to grow in climates that exhibit larger changes in the daily minimum and maximum temperatures reminiscent of desert-like climates. This fits with the known pattern that many plant species growing in the desert show rather shallow root systems with advantages including reduced energy input, increased ability to capture moisture from precipitation, and high nutrient availability in upper soil layers, which may improve survival in a more varied climate (Pierret et al., 2016; Schenk, 2008; Ogura et al., 2019).

Surprisingly, we found low or some non-significant correlations between the nine root traits and total annual precipitation, a factor we had initially used to select the accessions characterized in this study. This finding could be due to a variety of factors underlying the design of our study including the limited number of accessions analyzed, the specific selective pressures acting on Swedish Arabidopsis populations, where the majority of our characterized accessions originate (Long et al., 2013), and the well-watered conditions that we used to grow our plants. It may also be that annual precipitation is less important as a selective pressure acting on Arabidopsis root architecture compared to temperature. Indeed vapor pressure deficit (VPD), which is the driving force determining the rate of water loss during transpiration (Fricke, 2017), is exponentially related to temperature and modern crop plants appear to be particularly susceptible to elevated VPD (Lobell et al., 2014). To reliably link environmental parameters to RSA responses, it will be necessary to survey our accessions in conditions similar to their native environment or in conditions replicating environmental stresses. This will ultimately expand our understanding of how root architectural traits relate to different climate parameters and how the root system will react to a warming climate, both in agricultural and in natural ecosystems.

Materials and methods

Growth system

Rhizotron design

Request a detailed protocol

As described in Rellán-Álvarez et al., 2015, rhizotrons were built with two sheets of 1/8” clear abrasion-resistant polycarbonate plastic (Port Plastics, Portland, OR, Product: MarkolonAR) and were water jet cut (AquaJet LLC, Salem, OR) into 15 cm × 30 cm rectangles with 14 small gaps running down each side (clear_sheets.dxf). Two spacers (spacers.eps) were laser cut (Stanford Product Realization Lab) from 1/8” black cast acrylic (TAP Plastics, Mountain View, CA, and Calsak Plastics, Hayward, CA). Two rubber U-channels were cut into 29 cm long strips from Neoprene Rubber Trim, 5/16” Wide × 23/32” High Inside (McMaster-Carr, Elmhurst, IL, #8507K33). Additional adaptations for automation included: beveling the top of each polycarbonate sheet to a thin 45°̊ angle using a belt sander. Half of the polycarbonate sheets had 2.5 mm holes drilled in 5/8” down from the top edge and 3/16” in from each side. The small holes were used to screw a 1/16” aluminum water jet cut hook, outer dimensions 18 mm × 57 mm with a 13 mm diameter circle cutout (rhizotron-machine-hooks-13mm.pdf), to both sides of the sheet using 2×4 M3 stainless steel socket cap screws (McMaster-Carr, Elmhurst, IL, #91292A109). Not only did the hooks provide a way to pick up the rhizotron, they also were threaded through a rhizotron top (rhizo-bin-top-enlarged.pdf) cut (Pagoda Arts, San Francisco, CA) from 1/8” black acrylic, which served as a light shield and a bearing surface for the hanging rhizotron. The other half of the polycarbonate sheets had 8 mm holes drilled in centered at the same location and were used as counterparts to the sheets with tops and hooks.

Boxes and holders

Request a detailed protocol

Similar to Rellán-Álvarez et al., 2015, rhizotrons were placed in black polyethylene boxes with acrylic holders during plant growth and imaging. Modifications were as follows: 12”W × 18”L × 12”H black polyethylene boxes (Plastic-Mart, Fort Worth, TX, Part number R121812A) were outfitted with a top holder and inner divider system to grow 12 rhizotrons at a time. The inside dividers were composed of 1/8” black acrylic cut into five sheets (bin-b.pdf, bin-c.pdf, and bin-d.pdf), which interlock perpendicularly with a large middle sheet (bin-d.pdf), as well as two smaller pieces (bin-b.pdf, bin-c.pdf) at each end. This divider system creates chambers, which isolate the rhizotrons from each other, thus preventing light affecting neighboring rhizotrons when one rhizotron is removed. Each bin has a bin top (rhiz-bin-top-enlarged.pdf), also cut from 1/8” black acrylic, has 12 cutouts arranged in two rows of six, where the rhizotrons sit. The cutouts have two vertical tabs, which were glued in using TAP Acrylic Cement (TAP Plastics, Mountain View, CA) and provide placement guides for proper rhizotron insertion. The bin tops had two posts in diagonally opposite corners, which were connected with sticky copper tape (McMaster-Carr #76555A641) placed on the underside of each bin top to create a conductive path that the robotic arm used to find each bin and calculate the location of each rhizotron. Every post was made from an M5 stainless steel button head screw (McMaster-Carr #92095A484), 35 mm long – 13 mm OD 18-8 Stainless Steel Unthreaded Spacer for M5 screws (McMaster-Carr, #92871A276), 25 mm long – 10 mm OD 18-8 Stainless Steel Unthreaded Spacer for M5 screws (McMaster-Carr, #92871A093), and an M5 stainless steel lock nut (McMaster-Carr #93625A200). The screw went through the larger spacer, through the bin top, through the smaller spacer, and everything was held together by the nut. Each post extended down and into a 10 mm hole precisely positioned at the edge of each black box.

All build files for the rhizotrons and holders can be found at https://doi.org/10.5281/zenodo.6694558.

Biological components

Six accessions

Request a detailed protocol

Bay, Col-0, Est, Sha, Tottarp, and Tsu-0 were transformed using the floral dip method (Zhang et al., 2006) with the ProUBQ10::LUC2o reporter, which consists of the UBQ10 promoter region and first intron amplified from Col-0 genomic DNA, the plant-codon-optimized LUC2 coding sequence, and a plasma membrane-localized mCherry coding sequence driven by the 35S promoter (additional details in Rellán-Álvarez et al., 2015). Positive transformants were selected under a fluorescent dissecting scope (Leica M165 FC) using the mCherry marker visible in the mature seed. Bay, Col-0, Est, Sha, and Tsu-0 seeds were carried through to the T3 generation. Tottarp was used in the T2 generation.

Diversity panel

Request a detailed protocol

A. thaliana accessions were ordered from the Arabidopsis Biological Resource Center (ABRC: CS78885; Supplementary file 2). For vernalization, plants were grown in pots filled with Pro-Mix PGX soil (Premier Tech, Canada) and placed in a growth cabinet with long day conditions (16 hr light, 8 hr dark) at 10°C and 50% humidity using fluorescent bulbs with a light intensity of about 200 μE m−1 s−1. Upon flowering, plants were transformed as described above and then moved to a 22°C greenhouse with long day conditions and 150–250 μE m−1 s−1 light intensity. Transformation of the 192 initial accessions yielded 187 positive T1 lines, 171 of which could be confirmed in the T2 generation. Ten positive transformants from each of the T2 lines were screened on plates for strong luminescent signal once primary root length reached approximately 1.5 cm, usually around 2 weeks, but up to 26 days. Seedlings with strong root signal were given a line identification letter, then transferred to pots, and grown in the 10°C growth cabinet until inflorescences emerged (27–109 days after transfer to soil depending on the accession), at which point all plants were photographed and moved to a 22°C greenhouse until seeds could be harvested. Ultimately, 148 lines yielded at least three independent lines with strong root luminescence. Before starting an experiment, seeds were again screened for positive transformants using the mCherry seed coat marker.

Growth method

Rhizotron preparation

Request a detailed protocol

Rhizotron preparation was done as described in Rellán-Álvarez et al., 2015, with minor modifications: a polycarbonate plastic sheet with 8 mm holes and a sheet with the black rhizotron top and metal hooks were each laid on the table. Since the acrylic rhizotron top extends from the top in every direction, those sheets were laid down on the building surface with the acrylic piece and hooks lying just-off the edge of the surface and the screws pointing up, which ensured that the sheet lay flat. To prevent excess buildup of soil, the spacers were inserted into the sheets with 8 mm holes, and ‘flags’ (spacers modified with small pieces of material at the top to cover the bottom part of the hook and screw, Flags-final.eps) were inserted into the sheets with the attached tops. Using an electric paint sprayer (Wagner Spraytech Control Spray Double Duty HVLP Sprayer Model #0518050), a mist of water was applied to the sheets. Peat based Pro-Mix PGX soil (Premier Tech, Canada) was gently sifted over the sheets using a 2 mm sieve (US Standard Sieve Series N ̊ 10) and excess soil was gently shaken off. The mist, sift, shake procedure was done for a second time, which creates a two-layer thick surface of soil thin enough to see a small amount of light through when held up to a lightbulb or window. To prevent soil from falling out, a folded and trimmed piece of paper towel was moistened then placed at the bottom of each rhizotron. The spacers and flags were removed and the flags were replaced with clean spacers. The two sheets are carefully placed together and rubber U-channels were slipped on to each side. A small handful of sifted soil was placed into the groove in the rhizotron top and gently pushed down. Completed rhizotrons were hung in the boxes and, due to the increased size of the boxes, 3 L of water were added to the bottom of each box.

Sowing and growth

Request a detailed protocol

Two transfer pipettes (each ~2 mL) of Peters Professional 20-20-20 General Purpose fertilizer were added to each rhizotron immediately after building. Fertilized rhizotrons in filled boxes were covered and left to soak overnight. Seeds were stratified for 2 days at 4°C in distilled water. Three seeds were sown in the center of each rhizotron using a p200 pipette (Eppendorf) and a unique barcode was placed on the rubber U-channel. Rhizotrons were sprayed down with water and a clear acrylic sheet was placed on top of the box then sealed with packing tape to create a humid environment. Three days after sowing, clear sheets were unsealed and rhizotrons were watered with two transfer pipettes of water. The following day, the clear sheets were removed and rhizotrons were watered again. Rhizotrons were watered with two pipettes of water until 9 days after sowing. Plants were thinned to one central plant 5 days after sowing. Plants were grown under long day conditions (16 hr light, 8 hr dark) at 22/18°C (day/night) using LED lights (Valoya, C-series N12 spectrum) with a light intensity of about 130 μE m−1 s−1. For the six accessions, plants were grown until 9 or 21 DAS and then were imaged every day until 31 DAS. Five replicates of each accession were grown for each treatment and plants were watered with luciferin every 6 days. For the diversity panel, the population was grown up, as a whole, six times, thus having temporal population replicates rather than internal replicates in order to reduce variation effects of growth condition. Since we imaged the set of 93 accessions in two parts (parts 1 and 2), we mostly split the accessions alphabetically and added two Col-0 controls per imaging session. The rhizotrons of each imaging session were then swapped once between boxes to avoid rhizotron building effects. For the second imaging session (part 2), we made sure to include a Col-0 control to the bin that grew in the shelf outside of GLO-Bot (part 2, Bin D) to examine position effects. After the first randomization, accessions were always planted in the same locations relative to each other (Figure 4—figure supplement 2) so the population was treated identically in each replicate. Plants were grown on shelves in the growth chamber next to the imager until 14 DAS. Subsequently, they were transferred into the imager and imaged every other day from 14 to 28 DAS with luciferin additions every 6 days.

Plant imaging

Request a detailed protocol

On each imaging day (timing depending on experiment), GLO-Bot loaded each rhizotron into the imager then closed the door and triggered μManager to capture 5 min exposures on each side of the rhizotron. Inside the imaging system, there are two cameras on top of one another, with one taking an image of the top part of the rhizotron and the other, the bottom part of the rhizotron. After one side was imaged, the rhizotrons were rotated using a Lambda 10-3 Optical Filter Changer (Sutter Instrument, Novato, CA) and the other side of the rhizotron was imaged using the top and bottom camera, yielding four images per rhizotron. If it was the first imaging day or a designated luciferin day (every 6 days), GLO-Bot added 50 mL of 300 μM D-luciferin (Biosynth International Inc, Itasca, IL) to the top of each rhizotron immediately before loading the rhizotron into the imager. Before every replicate, a blank image was taken to account for background noise within the imager.

Image analysis

Request a detailed protocol

Initial root image preparation to combine and align the four raw images was done as described in Rellán-Álvarez et al., 2015. Initial background noise removal was done using remove.background.ijm. By manually opening the blank image for the respective replicate and then running this macro, the open blank image will be subtracted from all of the files in the folder being processed by the macro. The ImageJ plugin Template Matching and Slice Alignment was run as instructed in the manual to register all root images. Images are then run through the macro denoise_overlay.ijm to remove additional noise and add images together, which maximizes signal intensity throughout the root by: (1) subtracting 1.5 from all pixels to remove pixels with low values; (2) multiplying all values by 3, which greatly amplifies pixels with high values; (3) subtracting 8 from all values, which again removes pixels with low values (i.e. those that did not get amplified); (4) running the ImageJ command Despeckle, which replaces every pixel with the median value in its 3×3 neighborhood (Schindelin et al., 2012; Schneider et al., 2012); (5) setting and using a threshold to generate a binary image; (6) running the ImageJ commands Despeckle, Erode, and Dilate to further remove small pieces for noise; (7) using the ImageJ command Analyze Particles to get rid of small particles; and (8) merging the cleaned image with that of the previous day using the ImageJ command Add Create. This macro can be adjusted using the subtract-multiply-subtract sequence to change how signal intensity is amplified by manually testing the values on a handful of images.

Next, macros invert.ijm and tip_tracking.ijm are run to prepare images for GLO-RIAv2 (https://doi.org/10.5281/zenodo.5574925). This macro works by: (1) opening two consecutive images; (2) using the ImageJ command Dilate twice on the image from the earlier time point; (3) using the ImageJ command Subtract Create to subtract the earlier day from the later day. By dilating the earlier root system, we are able to ensure complete subtraction from the later time point and, therefore, only analyze the new root growth. Again, the thresholding in the new growth isolation can be refined if needed. Images were then run through GLO-RIAv2 and downstream analysis was performed in R (R Development Core Team, 2019) for additional cleanup and calculating traits. Initial formatting and combination of GLO-RIA output files is done using 1-format.Rmd. This file uses the raw output, the ‘key’, which contains information about the plants, and the experimental start date to calculate relative imaging times (very important for comparing multiple experiments), extract identifying information for each rhizotron, and, for the local files, calculate the root angles with respect to gravity using trigonometric calculations based on whether the raw output is greater than or less than 90°. Although we re-calculate the angles from 0 to 90°, having the 0 to 180° angles tells us which way the vector is pointing and which trigonometric function should be used to calculate the end point for the root segment. It should be emphasized that angles are with respect to gravity and not with respect to the parent root. Formatted files are run through 2-clean.R to remove stray particles. This process works by (1) computing the distances between all root x-y points in an image; (2) using hclust(method = “single”) to cluster those distances; (3) split the cluster tree into two groups using cutree(); (4) calculating the distance between the groups and determining if it is greater or less than the predetermined proximity maximum; (5) if the distance is smaller than that distance the cluster is kept, if it is larger than the distance, cluster with the larger minimum distance to the center-x is kept and the other cluster is further examined; (6) in the cluster being examined, the cluster is again tested for whether the distance between the clusters is far too large (nearly twice that of the initial proximity maximum) or whether there are less than four points in the cluster, which would indicate the cluster is likely noise. If either of these are true, the cluster is discarded. If not, the cluster is further examined to test if the points are linear, which would indicate the cluster could be a root segment. The above sequence is done for 1 day, then the next day is added in and the process begins again, and this continues until only the ‘true root’ segments remain. Parameter adjustments in this file primarily depend on the strength of the root signal and can be optimized through trial and error. The 2-clean.R script outputs a text file describing which clusters were kept and removed, a pdf for each rhizotron showing all of the removal steps, two csv files, clean_removed_points.csv and clean_true_root.csv, which list all of the points that were either kept or removed, as well as two initial trait files: clean_ROIs.csv and clean_traits.csv. Having these ‘true root’ segments allows for calculation of traits at the new growth level as well as the whole root level. Some of these calculations are done during the root cleaning process, while others are re-calculated later in 3-summarize_traits.Rmd.

Images were manually curated to eliminate those with very high levels of noise or from occasional camera malfunction. Additionally, phosphorescent (‘glow-in-the-dark’) stars in the upper right corner of some images were manually cropped out. Initially, these stars were placed in the images to aid image alignment, but instead generated too much noise. For the GLO-Roots diversity panel, a total of 34 images were removed (Supplementary file 6). Those images were replaced by an image of the previous day to ensure the correct number of images for GLO-RIAv2 analysis.

Trait processing

Request a detailed protocol

Raw trait measurements were converted into single values for each genotype in the software R (R Development Core Team, 2019) using the script mcmcglmm_best.R. This script is run for each trait, and the user inputs the raw phenotypic data for each trait. Several models were tested, including a second-order polynomial for time, different error correction factors, and models with and without kinship-informed random factor. Each model was examined for parameter convergences using the deviance information criteria (DIC), a Bayesian model version of Akaike information criterion (AIC). The selected model was chosen based on speed, simplicity, and generalizability. It should be noted that some traits may benefit from a higher order function. The final model had the form:

y=β0+(β1+ugenot)t+ugenoi+urep+e

where the trait (y) and time (t) were fitted as fixed effects, genotype (ugenoi), replicate (urep), and the time and genotype interaction (ugenot) were treated as random factors. The effects capture the average deviation for each genotype from the mean of the population. The script runs the data through a null model, which is a model run with no parameters and provides the baseline comparison for future models using the same input data, and then runs the previously described mixed model to estimate fixed effects by running 1,000,000 iterations in a Monte Carlo Markov chain (MCMC) and 200,000 burn-in using the MCMCglmm R package (Hadfield, 2010). Since MCMCglmm utilizes an MCMC walk, the model has inherent stochasticity and can produce slightly different values each time it is run. For reproducibility, we set a ‘seed number’ for the script then ran the model three times for each trait then selected the model with the best DIC.

Once the model is fit, we extract the genotype-specific parameters which capture differential starting values of a trait and differential trajectories over time. With these intercept and slope parameters, we can calculate the value of the trait for each genotype at any desirable time point, to understand the genotype’s trait behavior. We do this for the time points 0, 48, 96, 144, 192, 240, 288, and 336 hr to match our imaging frequency. These fitted genetic values are termed our ‘fitted values’: the coefficients (intercept and slope) as well as the expected value at each time point for each genotype. Broad-sense heritability was calculated using standard approaches in MCMCglmm (Hadfield, 2010), which essentially quantify variation on the trait trajectories associated to genotype relative to the total (Hadfield, 2017).

In addition, once we extracted genotype-specific values of each trait (either the intercept or slope coefficient or the predicted value at a time point), we conducted GWAs using GEMMA and the AraGWAS database SNP matrix (Zhou and Stephens, 2012; Togninalli et al., 2018). GEMMA was run as a linear mixed model (lmm) and with a kinship matrix to correct for accession relatedness. Only SNPs with a minimum allele frequency greater than 0.05 were examined. Bonferroni threshold, which was calculated using the number of linkage blocks within our population as computed by PLINK, was used as the significant threshold (0.05/113831 LD blocks). Manhattan plots were generated in R.

Data availability

GLORIAv2 is available through Zenodo, DOI: https://doi.org/10.5281/zenodo.5709009 https://doi.org/10.5281/zenodo.5574925. Image analysis pipelines and scripts are available through Zenodo, DOI: https://doi.org/10.5281/zenodo.5708430. RShiny App for exploring root system architecture of the six accessions is available through Zenodo, DOI: https://doi.org/10.5281/zenodo.6757675. RShiny App for exploring root system architecture of the diversity panel is available through Zenodo, DOI: https://doi.org/10.5281/zenodo.5708422. Imaging data and images are available through Zenodo, DOI: https://doi.org/10.5281/zenodo.5709009. General code for software operating robotics available: GitHub: https://github.com/rhizolab/rhizo-server. Rhizotron laser cutting files are available through Zenodo, DOI: https://doi.org/10.5281/zenodo.6694558. Previously published datasets used: WORLCLIM2: Fick SE, Hijmans RJ, 2017, https://worldclim.org/, https://doi.org/10.1002/joc.5086.

The following data sets were generated

References

Decision letter

  1. Caroline Gutjahr
    Reviewing Editor; Technical University of Munich, Germany
  2. Jürgen Kleine-Vehn
    Senior Editor; University of Freiburg, Germany
  3. Niklas Schandry
    Reviewer; Ludwig Maximilians University Muenchen, Germany
  4. Larry York
    Reviewer; Oak Ridge National Laboratory, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Uncovering natural variation in root system architecture and growth dynamics using a robotics-assisted phenomics platform" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jürgen Kleine-Vehn as the Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Niklas Schandry (Reviewer #1); Larry York (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

As you see from the reviewers' comments they (and myself) are excited about your manuscript, although it requires some revision and clarifications to become acceptable for eLife. Please find the reviewer comments below. We have discussed the experiment with 6 accessions (mentioned by reviewer #3) and we agree that it should stay in the manuscript but that it also causes some distraction from the main point. We suggest to more strongly focus the text on these results to better integrate them with the next chapter describing results from the full panel. Otherwise, we consider all reviewer comments to be important to improve the manuscript.

Reviewer #1 (Recommendations for the authors):

Comments on GWAS:

– I noticed that sometimes, the data cleaning procedure appears to retain pixels which are not part of the root system (e.g., Ting-1 in GWA-5; Tomegap2 in GWA-1 and GWA-5 or Dra2-1 in GWA3).

– The authors use “depth” to refer to the vertical distance covered by the overall root system, and “length” to refer to the sum of the length of all vectors that represent the root system. In the extreme case of a single main root growing straight down depth would be equal to cumulative length and depth should never be larger than cumulative length, yet there are cases where depth is larger than the total length of the root (Kz9, GWA3 and Kor3 GWA2). The length measurement in several accessions suffers from unconnected vectors, while the overall shape of the root system is probably correct (depth, width, hull). See for example accession Ge-0 GWA3, GWA5 or GWA6, where the root appears to be made up of many very short vectors with gaps. This raises the question if the length variable is suitable for quantitative genetic analysis. I also wonder how this issue affects the accuracy of the vector angles with respect to the true root system.

– Breeding values were computed for the quantitative genetic analysis. It should be clarified whether “breeding value” refers to the estimated model coefficients (p. 8), or to the fitted values (p. 54). All traits were overall using the same formula. It is not really clear why the authors opted for this approach instead of using different models for the different traits. Some traits, such as width (see Figure S5) do not display a linear relationship with time.

Average angle per day appears to extremely dynamic and I am not sure if the breeding values for this trait are meaningful with regards to the original data, raising the question if the fitted values for this trait are a good input for GWAS. Is it possible that the overall observation that p-values of some SNPs decrease over time simply reflects the linear relationship with time specified in the model, as fitted values increase over time?

What is the advantage of using fitted values at different timepoints over using the model coefficients for each accession as the input?

Since the model assumes a linear relationship, which appears to not capture the dynamics of all traits, would it not be reasonable to simply compute the trait mean per day per accession and use these values as input for GWAS, instead of using fitted values?

Reviewer #2 (Recommendations for the authors):

1. While the GLO-Bot system is using potting mix, it will still be very different from soil grown roots. With GLO-roots, the roots are grown in a thin sheet of soil (essentially restricted to 2 dimensions) and the roots will become constrained by the width of the system as well. This will most likely reduce differences between shallow rooting accessions and intermediate accessions. It will be important to discuss the constraints of the system and discuss whether any time point will be closer to a 3D soil growth environment. One way how this might be elucidated would be to look at broad sense heritability at different time points. Such an approach might also indicate an optimal timepoint for the GWAS and the correlations with environmental variables.

2. Page 8: "For each of these traits, we used a generalized linear mixed model in the R package MCMCglmm to account for replicate and block effects noises and estimate the "breeding value" of each accession for a given trait (Wilson et al. 2010; Mrode 2014); (Hadfield 2010) (Supplemental Figure 5). The denoised trait values for each genotype at 0, 48, 96, 144, 192, 240, 288, and 336 hours after the start of imaging (equivalent to 14, 16, 18, 20, 22, 24, and 28 DAS) were used for further analyses."

It remains a bit unclear what the process was (how were the block effects defined, etc.) and what the rationale was to use "breeding value". Also, most readers will not know what a "breeding value" is. Are the "denoised" trait values breeding values? Why do the authors think that the raw trait values are not helpful for further analysis (after all, many of the BSH seems reasonable)?

3. Page 9ff: The candidate gene follow up is very incomplete. For instance, the authors show in their Manhattan plots that some associations are linked to other SNPs in that region. It therefore would be good to also include genes in proximity of the linked SNPs in the table. Also, the SNPs identified with a GWAS are not necessarily the causal SNPs (in particular, this is an issue if not the full genome sequence was used for GWAS). The causal SNPs could be linked or partially linked to the GWAS SNPs (e.g. due to allelic heterogeneity). It would therefore be helpful if for the genes discussed in the text, the genomic region could be shown with haplotypes or just using the 1001 genome project SNPs in accessions with alternative GWAS SNPs.

4. Page 5: Root Angle output is unclear (also in the Suppl. Table). This is important as it is a major trait in the manuscript. Is this the average of all roots? Do all roots need to be detected for that? What are the two points that are used to measure the angles? What is a root segment? How many are considered?

5. Page 6: How many root systems were compared with SmartRoot? Is each point in Figure 2B an entire root system? Is the comparison similarly accurate for very different RSAs? What is the advantage of GloRoot vs Smartroot (i.e. how much time does it take to trace a GLO-bot RSA manually)?

6. The trait data is not deposited with the manuscript or on a public repository. This would be very important to allow the community to build on the work of the authors.

Reviewer #3 (Recommendations for the authors):

The comments below may be helpful to the authors.

Since there weren’t line numbers, I give page numbers below and tried to ‘quote’ key words that should allow the authors to find what my comments pertain to.

1: ‘above ground’ is not hyphenated here but is elsewhere.

1: 'climate variables of the accessions respective origins' or something as to not be confused with growth conditions, etc.

2: 'between' should be among.

2: "Most RSA studies…" – can the most here really be substantiated? If not, "many" may suffice. I am not sure that most are in gel or in Arabidopsis.

3: rhizotrons are often called rhizoboxes. If rhizotron is preferred, perhaps "(or rhizobox)" can be added at first mention to orient some readers.

5: 'detection and translation of subsequent images' – I believe this describes what is commonly called registration so the word may be useful to add here for clarity.

5: 'invert the images' may be expressed as 'invert the colors of the images' to be precise.

5: 'root segment' needs defined.

6: 'series of vectors' – seems like RSML and some would ask if RSML is supported.

6: 'half of which' could be confused for accession or replicate (I know it's replicate) but consider saying "where five replicates" instead.

6: for the 6 accessions, you can add a bit more quantification by using a calculation of heritability or repeatability. You could have also used the additional replication for power analysis to determine the need for replications numbers, for example boot strapping from 2 – 10 reps to see if replicates improve measures or heritability.

7: in the last paragraph, what type of randomization used should be specified (maybe in methods, I don't like this format where some methods are in results but not all, but I know it's because results are first). On page 8 it is mentioned that an R package is used for block. I think this is leading me to conclude the results simply have to methods mixed in.

8: So is it right that a subset of days are used? Why?

8: PCA over time implies time dependent traits are just sized related.

9: ‘fitted values’ could be clarified.

10: ‘published bioclimatic variables’ could include for clarity ‘…from the origin sites of the respective accessions.’

11: ‘progress in machine-learning’ – I remember I think one of the original rationales for this system was to image on both sides for a more complete root system? Was that in the original paper? Is that done here? With regards to bioluminescence – is it still clear that this method is advantageous to RGB imaging or other forms of imaging? For example, the Smith et al. 2020 paper shows many examples of identifying roots from complex backgrounds. The RootPainter software has been used for a rhizobox study:

https://www.biorxiv.org/content/10.1101/2021.08.24.457510v1

Can the possible advantages be discussed, especially if there are specific examples from images that show color imaging alone is not sufficient?

11: sample sizes could be sample numbers I think.

Figure 1 – scale needed in all panels.

48: Not clear why Materials and methods are in supplemental?

51: "Accessions were always planted in the same locations relative to each other so the population was treated identically in each replicate." Strictly speaking – this is a bad choice, somehow you decided to intentionally NOT randomize which goes against statistical training and a strict reviewer may say it compromises the experimental design where all the measurements could be artifacts of position. Is there a good reason to not randomize? It seems odd in what otherwise seems like such careful design.

51: 'four raw images' – root and shoot from each side is four? Maybe make clear.

51: what is the 'correct blank images' – is this process entirely automated? Sounds like a user walks the processing through several fairly automated steps but not quite automated from image to data it seems.

https://doi.org/10.7554/eLife.76968.sa1

Author response

As you see from the reviewers' comments they (and myself) are excited about your manuscript, although it requires some revision and clarifications to become acceptable for eLife. Please find the reviewer comments below. We have discussed the experiment with 6 accessions (mentioned by reviewer #3) and we agree that it should stay in the manuscript but that it also causes some distraction from the main point. We suggest to more strongly focus the text on these results to better integrate them with the next chapter describing results from the full panel. Otherwise, we consider all reviewer comments to be important to improve the manuscript.

Reviewer #1 (Recommendations for the authors):

Comments on GWAS:

– I noticed that sometimes, the data cleaning procedure appears to retain pixels which are not part of the root system (e.g., Ting-1 in GWA-5; Tomegap2 in GWA-1 and GWA-5 or Dra2-1 in GWA3).

We thank reviewer #1 for this comment. It indeed happens that the data cleaning procedure does not remove all background pixels. The algorithm cannot distinguish background pixels from root pixels and it needed the mentioned adjustments to retain as much root signal as possible and remove as much background noise as possible. Since we wanted all images to be denoised in the same way, we had to compromise with some images showing remaining background noise. A list of images with these problems has been curated and is included in Supplementary File 3. Certain traits, such as average angle, have been weighted by vector length so tiny measurements from the noise do not have a strong influence on the trait value. Further, we expect noise to not be systematically associated with certain genotypes, rather it will be randomly scattered across replicates within a genotype and across all the genotypes. This measurement error would therefore result in lower heritability estimates, making the conclusions we draw naturally more conservative.

– The authors use “depth” to refer to the vertical distance covered by the overall root system, and “length” to refer to the sum of the length of all vectors that represent the root system. In the extreme case of a single main root growing straight down depth would be equal to cumulative length and depth should never be larger than cumulative length, yet there are cases where depth is larger than the total length of the root (Kz9, GWA3 and Kor3 GWA2). The length measurement in several accessions suffers from unconnected vectors, while the overall shape of the root system is probably correct (depth, width, hull). See for example accession Ge-0 GWA3, GWA5 or GWA6, where the root appears to be made up of many very short vectors with gaps. This raises the question if the length variable is suitable for quantitative genetic analysis. I also wonder how this issue affects the accuracy of the vector angles with respect to the true root system.

We thank reviewer #1 for addressing this point. They are correct, the total length measurements come with caveats due to the nature of the unconnected vectors and therefore the unaccounted length between detected root segments. This caveat about the length measurement and other potential caveats have been included in the text on page 6/7. For the vector angles, each vector angle is weighted by the length of the segment, so angles computed from very small short segments are weighted differently than longer connected segments. While it is true that the angles for fully connected roots could differ, this is the most objective and robust approach we could devise and, given the high heritability of the trait, we believe the signal is well captured.

– Breeding values were computed for the quantitative genetic analysis. It should be clarified whether “breeding value” refers to the estimated model coefficients (p. 8), or to the fitted values (p. 54). All traits were overall using the same formula. It is not really clear why the authors opted for this approach instead of using different models for the different traits. Some traits, such as width (see Figure S5) do not display a linear relationship with time.

Average angle per day appears to extremely dynamic and I am not sure if the breeding values for this trait are meaningful with regards to the original data, raising the question if the fitted values for this trait are a good input for GWAS. Is it possible that the overall observation that p-values of some SNPs decrease over time simply reflects the linear relationship with time specified in the model, as fitted values increase over time?

What is the advantage of using fitted values at different timepoints over using the model coefficients for each accession as the input?

Since the model assumes a linear relationship, which appears to not capture the dynamics of all traits, would it not be reasonable to simply compute the trait mean per day per accession and use these values as input for GWAS, instead of using fitted values?

We thank the reviewer for the point regarding breeding values. We realize this was not explained in enough detail. The reviewer is correct about the estimated per-genotype trait deviations and the terminology in the text was unclear. For the sake of clarity, the terminology has been changed to only use “fitted values''. Using our temporal model, the slope and intercept coefficients were used to calculate the expected genotype value at any timepoint. The text on page 8 (now page 9) was clarified.

Both using the accessions’ coefficients as well as the predicted value at a given timepoint are fair and interesting approaches to study the differences in root traits among ecotypes, and each provides different conclusions. Simply using the coefficients for GWA will tell us about the genes responsible for the main trends in the data, but these are less intuitive. Conducting GWA on the predicted values at each timepoint gives us an understanding about when a trait is most relevant, which is an important conclusion of our efforts to create a longer time series of trait data than previously available.

Another question from the reviewer regards the different possible functions used to model traits; this is a matter of statistical tractability and trade-offs between different models. While higher order polynomials or sigmoid curves may fit well or even better than a linear relationship of the trends, these models were not tractable in our experimental and analysis design. Let us give some background: To make the most of the data, the approach we took is that of time series analyses of Bayesian generalized linear mixed models where we can include experimental design factors to control for batch effects and at the same time fit genotype factors to estimate genotype-specific trends (breeding values of the coefficients). This is advantageous because we can control for noise while using all the replicates and all the time points to infer genotypic effects (this is a common approach in complex breeding programmes and designs with temporal datasets to extract breeding values, see: http://cran.nexr.com/web/packages/MCMCglmm/vignettes/CourseNotes.pdf). The limitation is that when temporal functions (polynomial, etc.) become more complex, the model requires fitting extra parameters (number of genotypes * new function parameter). This makes it infeasible to fit very complex models, and indeed when we tried using MCMCglmm the models often did not converge to a solution. One may be tempted to fit a polynomial model for all the data points for each genotype separately, which would be feasible. However, this approach would not enable calculating breeding values for accessions, and would lack control for batch effects. Given these trade-offs of model selection, we decided to sacrifice complexity of the time series function in favor of better control for experimental design.

Reviewer #2 (Recommendations for the authors):

1. While the GLO-Bot system is using potting mix, it will still be very different from soil grown roots. With GLO-roots, the roots are grown in a thin sheet of soil (essentially restricted to 2 dimensions) and the roots will become constrained by the width of the system as well. This will most likely reduce differences between shallow rooting accessions and intermediate accessions. It will be important to discuss the constraints of the system and discuss whether any time point will be closer to a 3D soil growth environment. One way how this might be elucidated would be to look at broad sense heritability at different time points. Such an approach might also indicate an optimal timepoint for the GWAS and the correlations with environmental variables.

We thank reviewer #2 for this very important point. We agree with the reviewer that the dimensions of the growth vessel will impact the architecture of the root system, as will soil type. While measurements of root growth in large volumes of natural soil will improve the ecological significance of the trait measurements, we are not aware of data supporting the view that this will necessarily affect the heritability of these traits. Likely, heritability will be affected by the complexity of the environment that roots grow; conditions such as agar gel will likely exhibit tighter relationships between genotype and phenotype because there is less heterogeneity in the environment. This fundamentally differs from more natural soil conditions where heterogeneity in the environment is intrinsic to the growth condition. We believe that the relationship between genotype and phenotype has greater ecological significance when measured using the GLO-Roots system due to the more physiologically relevant ways in which light, water and nutrient availability interact with the root system and processes such as light, transpiration and gas exchange impact the plant as a whole. Certainly, as soon as the roots touch the edges of the rhizotron, traits like root depth and root width will lose significance, however traits like average angle or length per day will still show differences until the rhizotron is filled with roots.

We conducted an analysis of heritability over time for our traits, although we had not included it in the text. One could conduct MCMCglmm for each timepoint. Another alternative is to study the proportion of variance explained (PVE) of each trait by timepoint by the genome in the GWA analyses. GEMMA provides such PVE, a parameter that is also often called Chip heritability or narrow sense heritability, as the kinship matrix, built to control for genotype relatedness, captures additive differences. Similarly as Figure 6A, we can see how the PVE changes over time, see Author response image 1. It is interesting to note the relatively high amount of variance explained for traits such as weighted average angle and weighted average angle per day. These data are interesting leads to continue working on in the future.

Author response image 1

2. Page 8: "For each of these traits, we used a generalized linear mixed model in the R package MCMCglmm to account for replicate and block effects noises and estimate the "breeding value" of each accession for a given trait (Wilson et al. 2010; Mrode 2014); (Hadfield 2010) (Supplemental Figure 5). The denoised trait values for each genotype at 0, 48, 96, 144, 192, 240, 288, and 336 hours after the start of imaging (equivalent to 14, 16, 18, 20, 22, 24, and 28 DAS) were used for further analyses."

It remains a bit unclear what the process was (how were the block effects defined, etc.) and what the rationale was to use "breeding value". Also, most readers will not know what a "breeding value" is. Are the "denoised" trait values breeding values? Why do the authors think that the raw trait values are not helpful for further analysis (after all, many of the BSH seems reasonable)?

We realize that the rationale of the MCMCglmm approach was not described in enough detail. We have now included a much more extended version in the trait processing section of the methods. Further, we clarified the description highlighted in this comment on page 9 to be more consistent with the rest of the text. We have removed the term “breeding value” and instead use the general term “fitted values” to represent all of the fitted genetic values for each genotype: the slope, intercept, and the expected values at each time point for each genotype.

3. Page 9ff: The candidate gene follow up is very incomplete. For instance, the authors show in their Manhattan plots that some associations are linked to other SNPs in that region. It therefore would be good to also include genes in proximity of the linked SNPs in the table. Also, the SNPs identified with a GWAS are not necessarily the causal SNPs (in particular, this is an issue if not the full genome sequence was used for GWAS). The causal SNPs could be linked or partially linked to the GWAS SNPs (e.g. due to allelic heterogeneity). It would therefore be helpful if for the genes discussed in the text, the genomic region could be shown with haplotypes or just using the 1001 genome project SNPs in accessions with alternative GWAS SNPs.

We thank the reviewer for this point. We know the SNPs associated with the trait are not necessarily the causal ones. We have included the list of all genes within a 20 Kb region of the significant SNPs (10 kb upstream and 10kb downstream from the SNP) in Supplementary File 4 also listing the reference and alternative alleles. This has been clarified in the text, as well.

4. Page 5: Root Angle output is unclear (also in the Suppl. Table). This is important as it is a major trait in the manuscript. Is this the average of all roots? Do all roots need to be detected for that? What are the two points that are used to measure the angles? What is a root segment? How many are considered?

Root angles are computed on sections of roots visible in the image (aka “segments”). These segments do not always correspond to a biological structure, but rather the visible readout of the root system structure within the GLO-Bot platform. These root system structures visible on the GLO-Bot images are indeed disconnected in many locations, due to variability in the reporter’s intensity and the presence of physical obstructions such as soil particles. We actually used this disconnection to our advantage to compute more faithful root angle measurements within the image, as we use the “segments” as individual elements (see Author response image 2). We have added a detailed clarification to the text on page 6.

Author response image 2

5. Page 6: How many root systems were compared with SmartRoot? Is each point in Figure 2B an entire root system? Is the comparison similarly accurate for very different RSAs? What is the advantage of GloRoot vs Smartroot (i.e. how much time does it take to trace a GLO-bot RSA manually)?

The growing root systems of 5 plants were traced with SmartRoot. Since each root was imaged over a series of days and occasionally imaged multiple times within one day each plant has between 15 and 22 root images. In total, 95 root images were traced. Each point in the plot is a time point. Manually tracing a single GLO-bot RSA image with SmartRoot can take days, depending on the size of the root system. SmartRoot relies on continuous structures to follow the root. Due to the nature of the GLO-bot RSA images, this becomes a very tedious process of tracing each root system. We added more details to the main text and the figure legend of Figure 2B, stating: “Manual SmartRoot tracing of a single GLO-Roots image can take multiple days, thus making analysis of large-scale time course experiments unfeasible.”

6. The trait data is not deposited with the manuscript or on a public repository. This would be very important to allow the community to build on the work of the authors.

The trait data is available on Zenodo (https://doi.org/10.5281/zenodo.5709009) as stated at the end of the manuscript. A link has also been inserted into the text for clarity.

Reviewer #3 (Recommendations for the authors):

The comments below may be helpful to the authors.

Since there weren’t line numbers, I give page numbers below and tried to ‘quote’ key words that should allow the authors to find what my comments pertain to.

1: ‘above ground’ is not hyphenated here but is elsewhere.

We thank reviewer #3 for this comment and corrected this mistake. It now reads; ”The plant kingdom contains a stunning array of complex morphologies easily observed above-ground, but more challenging to visualize below-ground.”

1: 'climate variables of the accessions respective origins' or something as to not be confused with growth conditions, etc.

We have added in ‘of the accessions’ respective origins’ for clarity.

2: 'between' should be among.

Thank you, we changed “between” to “among”.

2: "Most RSA studies…" – can the most here really be substantiated? If not, "many" may suffice. I am not sure that most are in gel or in Arabidopsis.

We thank the reviewer for this valuable input and agree that “most” should be replaced by “many”.

3: rhizotrons are often called rhizoboxes. If rhizotron is preferred, perhaps "(or rhizobox)" can be added at first mention to orient some readers.

Thank you for this suggestion. We prefer rhizotrons to avoid confusion with the growth boxes that hold the rhizotrons in our setup. We added “rhizoboxes” to the text. It now reads; “In addition to GLO1, the frame contains a growth area for seven black polyethylene growth boxes (Figure 1 B), each of which holds 12 rhizotrons (also known as rhizoboxes) arranged in a two by six grid (Figure 1 C, Figure 1—figure supplement 1B).

5: 'detection and translation of subsequent images' – I believe this describes what is commonly called registration so the word may be useful to add here for clarity.

We agree with the reviewer and adapted the text. It now reads: “This plugin works by user selection of a region present in all images followed by automated detection and translation of subsequent images, a method that is also known as image registration.”

5: 'invert the images' may be expressed as 'invert the colors of the images' to be precise.

For clarity, this has been revised to: “…invert the gray-scale values of the images, thus making the root system itself black and the background white for downstream whole root system analyses.”

5: 'root segment' needs defined.

Root angles are computed on sections of roots visible in the image (aka “segments”). These segments do not always correspond to a biological structure, but rather the visible readout of the root system structure within the GLO-Bot platform. These root system structures visible on the GLO-Bot images are indeed disconnected in many locations, due to variability in the reporter’s intensity. We actually used this disconnection to our advantage to compute the root angle distribution within the image, as we use the “segments” as individual elements (see figure in comment 2.4). This also has been clarified in the text on page 6.

6: 'series of vectors' – seems like RSML and some would ask if RSML is supported.

We choose not to, since in the analysis we lose most of the topological information that is central in the RSML structure.

6: 'half of which' could be confused for accession or replicate (I know it's replicate) but consider saying "where five replicates" instead.

Thank you for this valid concern, we changed the text accordingly. It now reads; ”Each accession was grown in ten replicates, where five replicates were imaged daily from 9 days after sowing (DAS) until 31 DAS and the other five replicates from 21 – 31 DAS.”

6: for the 6 accessions, you can add a bit more quantification by using a calculation of heritability or repeatability. You could have also used the additional replication for power analysis to determine the need for replications numbers, for example boot strapping from 2 – 10 reps to see if replicates improve measures or heritability.

We thank reviewer #3 for these comments. We calculated heritability for the 6 accessions but chose not to include it since the small number of genotypes precludes confident heritability determination. We recognize the point about determining replicate numbers, however this was not something we did at the time. While the 6 accessions provided insight into some of the methods (e.g. luciferin watering speed and frequency, optimal imaging time window, and ground truthing image analysis methods), they did not inform factors such as replicate number in the experiment. Based on reviewer comments, we are choosing to not only emphasize the 6 accessions as a foundation for the diversity panel, but also highlighting them as an alternative example of how the GLO-bot system can be used (i.e. few accessions vs many accessions, high imaging frequency vs every other day imaging, raw data vs. fitted data, etc.).

7: in the last paragraph, what type of randomization used should be specified (maybe in methods, I don't like this format where some methods are in results but not all, but I know it's because results are first). On page 8 it is mentioned that an R package is used for block. I think this is leading me to conclude the results simply have to methods mixed in.

The 93 accessions of the GLO-Roots diversity panel were selected balancing genetic similarity and environmental variation along with reporter intensity, and the number of independent transformation events recovered. Since we imaged the set of 93 accessions in two parts (part 1 and 2), we mostly split the accessions alphabetically and added two Col-0 controls per imaging session. The rhizotrons of each imaging session were then swapped once between boxes to avoid rhizotron building effects. For the second imaging session (part 2), we made sure to include a Col-0 control to the bin that grew in the shelf outside of GLO-Bot (part 2, Bin D) to examine position effects. We now list the location of every accession in Figure 4—figure supplement 2. We added this information to the methods section on page 55.

We removed the comment about the block effects from the text since the final model did not take positional effects into account. It is true that there are more methods mixed into the results than usual. Given the “tools and resources” style of paper, we tried to include key pieces of methodology to help the reader understand the results.

8: So is it right that a subset of days are used? Why?

The six accessions demonstrate how GLO-bot can be used for continuous imaging of a small number of accessions with many replicates. With this in mind, the rhizotrons are imaged every day for a long period of time. Since the diversity panel had so many more accessions the plants were imaged for a subset of the days to account for possible hardware errors as they arose. We added this important information to the text on page 9.

8: PCA over time implies time dependent traits are just sized related.

We agree with this and added in “and plant size”.

9: ‘fitted values’ could be clarified.

We agree with the reviewer and have clarified the language regarding the fitted values. The description of fitted value calculation on page 9 has been revised to explain how the values were calculated and the manuscript has been reviewed to use consistent terminology. Figures have also been updated to indicate whether fitted values or raw values were used.

10: ‘published bioclimatic variables’ could include for clarity ‘…from the origin sites of the respective accessions.’

We appreciate this suggestion and have added the phrase in for clarity.

11: ‘progress in machine-learning’ – I remember I think one of the original rationales for this system was to image on both sides for a more complete root system? Was that in the original paper? Is that done here? With regards to bioluminescence – is it still clear that this method is advantageous to RGB imaging or other forms of imaging? For example, the Smith et al. 2020 paper shows many examples of identifying roots from complex backgrounds. The RootPainter software has been used for a rhizobox study:

https://www.biorxiv.org/content/10.1101/2021.08.24.457510v1

The reviewer is correct that imaging from both sides of the rhizotron is an important aspect of the GLO-Roots/GLO-bot methodology. This imaging itself in the GLO-Bot system is performed in the same manner as the original paper, so imaging from both sides does indeed happen here and continues to be advantageous for extensively capturing the entire root system. Given the size of Arabidopsis roots, the bioluminescence continues to be advantageous. While the field of machine-learning for identifying roots is rapidly progressing, advances seem to be primarily focused on plant systems with much thicker roots. To our knowledge, RootPainter is still not sufficient for observing Arabidopsis roots in RGB images. We have added the comparison to RGB imaging to the discussion on page 13.

Can the possible advantages be discussed, especially if there are specific examples from images that show color imaging alone is not sufficient?

In the original GLO-Roots paper (Rellan-Alvarez et al., 2015), we scanned the rhizotrons in order to measure the moisture level of the soil. Since Arabidopsis has very thin roots, RGB imaging will not reveal the whole root system (see for example Ogura et al., 2019, https://www.cell.com/cell/pdf/S0092-8674(19)30684-1.pdf). Furthermore, manual tracing is required since the contrast between the roots and the soil is not strong enough. This approach would not allow the throughput and time-lapse aspect that we were aiming for. The advantage of luminescence is that it makes fine roots appear thicker, that it shines through thin layers of soil and that it improves the contrast between root and soil signal. We used Arabidopsis because of its big collection of natural variation, genetic resources, plant size and generation time in order to possibly identify genes that are involved in root system architecture growth. We aimed at visualizing development from young root systems to mature root systems to observe if root traits like average angle change over time. For species with bigger roots, however, RootPainter would be a very interesting pipeline to use and would save the step of plant transformation, reduce imaging time and image processing. We added RGB imaging as a possibility for bigger root systems to the discussion (page 13).

11: sample sizes could be sample numbers I think.

We agree with the reviewer and changed the text accordingly.

Figure 1 – scale needed in all panels.

We have added scale bars to all relevant parts of this figure.

48: Not clear why Materials and methods are in supplemental?

We thank the reviewer for this comment. Since this manuscript is not layouted yet, we decided that the figures and supplemental figures and tables should follow the main text. In the layouted version, however, the figures, the figure supplements and supplementary files will be included in the text so that the Materials and methods would follow the discussion. We agree that this is important since the methods section is essential, especially in an article that describes new methods.

51: "Accessions were always planted in the same locations relative to each other so the population was treated identically in each replicate." Strictly speaking – this is a bad choice, somehow you decided to intentionally NOT randomize which goes against statistical training and a strict reviewer may say it compromises the experimental design where all the measurements could be artifacts of position. Is there a good reason to not randomize? It seems odd in what otherwise seems like such careful design.

We thank reviewer #3 for this very valid criticism. Before we started imaging the GLO-Roots diversity panel, we carefully discussed possible strategies. We discussed if all 6 replicates should be imaged within one experiment (imaging run) or if all 93 accessions should be included. Furthermore, we discussed if we should randomize the positions of rhizotrons or if we should keep the same order. We decided for the latter to be able to examine positional effects. Since we included 4 Col-0 controls distributed among the 8 rhizotron boxes, we believed that we could examine positional effects more efficiently. Rhizotron positions for the diversity panel are mapped out in the new Figure 4—figure supplement 2.

51: 'four raw images' – root and shoot from each side is four? Maybe make clear.

Within the imager, there are two cameras on top of one another, taking an image of the top part of the rhizotron and one of the bottom part of the rhizotron. After those images were taken, the stage turns the rhizotron and the cameras take another top and bottom image, yielding 4 images per rhizotron. This was clarified in the text on page 55.

51: what is the 'correct blank images' – is this process entirely automated? Sounds like a user walks the processing through several fairly automated steps but not quite automated from image to data it seems.

This assumption is correct. One blank image was taken for each replicate. The blank image for the respective replicate is manually selected and then the automated script runs and subtracts the inherent background noise. This was clarified in the text on page 55.

https://doi.org/10.7554/eLife.76968.sa2

Article and author information

Author details

  1. Therese LaRue

    1. Department of Biology, Stanford University, Stanford, United States
    2. Department of Plant Biology, Carnegie Institution for Science, Stanford, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Heike Lindner
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6805-5118
  2. Heike Lindner

    1. Department of Plant Biology, Carnegie Institution for Science, Stanford, United States
    2. Institute of Plant Sciences, University of Bern, Bern, Switzerland
    Contribution
    Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Therese LaRue
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5913-9803
  3. Ankit Srinivas

    Department of Plant Biology, Carnegie Institution for Science, Stanford, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  4. Moises Exposito-Alonso

    1. Department of Biology, Stanford University, Stanford, United States
    2. Department of Plant Biology, Carnegie Institution for Science, Stanford, United States
    Contribution
    Software, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5711-0700
  5. Guillaume Lobet

    1. UCLouvain, Faculty of Bioengineering, Louvain-la-Neuve, Belgium
    2. Forschungszentrum Jülich, Agrosphere Institute, Juelich, Germany
    Contribution
    Resources, Software, Funding acquisition, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  6. José R Dinneny

    Department of Biology, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    dinneny@stanford.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3998-724X

Funding

U.S. Department of Energy (DE-SC0008769)

  • José R Dinneny

U.S. Department of Energy (DE-SC0018277)

  • José R Dinneny

National Institutes of Health (T32GM007276)

  • Therese LaRue

Deutsche Forschungsgemeinschaft (LI 2776/1-1)

  • Heike Lindner

National Institutes of Health (1DP5OD029506-01)

  • Moises Exposito-Alonso

U.S. Department of Energy (DE-SC0021286)

  • Moises Exposito-Alonso

Deutsche Forschungsgemeinschaft (EXC-2070 - 390732324)

  • Guillaume Lobet

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Work in the JRD lab was funded by the US Department of Energy’s Office of Biological and Environmental Research (DE-SC0008769 and DE-SC0018277) and the Carnegie Institution for Science Endowment. TSL was funded by the National Science Foundation Graduate Research Fellowship and the National Institutes of Health Predoctoral Training Grant (T32GM007276). HL was funded by a German Research Foundation (DFG) research fellowship (LI 2776/1-1). MEA is funded by the Carnegie Institution for Science, an NIH Early Investigator Award (1DP5OD029506-01), and a DOE BER grant (DE-SC0021286). GL is funded by the German Research Foundation under Germany’s Excellence Strategy, EXC-2070–390732324 (PhenoRob). Part of the computational analyses were done at the Carnegie High-Performance Computing clusters Memex and Calc.

We thank Peter Sand and Zoltan DeWitt from Modular Science for the design and maintenance of GLO-Bot. We thank Wolfgang Busch for advice on the initial selection of Arabidopsis accessions. We thank Chris Ballard from Sutter Instrument for maintenance of GLO1. We thank Theo van de Sande and Ismael Villa for maintenance of the growth room. We thank Michael Raissig and Josep Vilarrasa-Blasi for constructive suggestions on the manuscript.

Senior Editor

  1. Jürgen Kleine-Vehn, University of Freiburg, Germany

Reviewing Editor

  1. Caroline Gutjahr, Technical University of Munich, Germany

Reviewers

  1. Niklas Schandry, Ludwig Maximilians University Muenchen, Germany
  2. Larry York, Oak Ridge National Laboratory, United States

Publication history

  1. Preprint posted: November 13, 2021 (view preprint)
  2. Received: January 11, 2022
  3. Accepted: August 28, 2022
  4. Accepted Manuscript published: September 1, 2022 (version 1)
  5. Version of Record published: September 21, 2022 (version 2)

Copyright

© 2022, LaRue, Lindner et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 712
    Page views
  • 272
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Therese LaRue
  2. Heike Lindner
  3. Ankit Srinivas
  4. Moises Exposito-Alonso
  5. Guillaume Lobet
  6. José R Dinneny
(2022)
Uncovering natural variation in root system architecture and growth dynamics using a robotics-assisted phenomics platform
eLife 11:e76968.
https://doi.org/10.7554/eLife.76968

Further reading

    1. Ecology
    2. Plant Biology
    Laura Stefan, Nadine Engbersen, Christian Schöb
    Research Article

    By capitalising on positive biodiversity-productivity relationships, intercropping provides opportunities to improve agricultural sustainability. Intercropping is generally implemented using commercial seeds that were bred for maximal productivity in monocultures, thereby ignoring the ability of plants to adapt over generations to the surrounding neighbourhood, notably through increased complementarity, i.e. reduced competition or increased facilitation. This is why using monoculture-adapted seeds for intercropping might limit the benefits of crop diversity on yield. However, the adaptation potential of crops and the corresponding changes in complementarity have not been explored in annual crop systems. Here we show that plant-plant interactions among annual crops shifted towards reduced competition and/or increased facilitation when the plants were growing in the same community type as their parents did in the previous two generations. Total yield did not respond to this common coexistence history, but in fertilized conditions, we observed increased overyielding in mixtures with a common coexistence history. Surprisingly, we observed character convergence between species sharing the same coexistence history for two generations, in monocultures but also in mixtures: the six crop species tested converged towards taller phenotypes with lower leaf dry matter content. This study provides the first empirical evidence for the potential of parental diversity affecting plant-plant interactions, species complementarity and therefore potentially ecosystem functioning of the following generations in annual cropping systems. Although further studies are required to assess the context-dependence of these results, our findings may still have important implications for diversified agriculture as they illustrate the potential of targeted cultivars to increase complementarity of species in intercropping, which could be achieved through specific breeding for mixtures.

    1. Cell Biology
    2. Plant Biology
    Alexis Lebecq, Mehdi Doumane ... Marie-Cécile Caillaud
    Research Article

    Membrane lipids, and especially phosphoinositides, are differentially enriched within the eukaryotic endomembrane system. This generates a landmark code by modulating the properties of each membrane. Phosphatidylinositol 4,5-bisphosphate [PI(4,5)P2] specifically accumulates at the plasma membrane in yeast, animal, and plant cells, where it regulates a wide range of cellular processes including endocytic trafficking. However, the functional consequences of mispatterning PI(4,5)P2 in plants are unknown. Here, we functionally characterized the putative phosphoinositide phosphatase SUPPRESSOR OF ACTIN9 (SAC9) in Arabidopsis thaliana (Arabidopsis). We found that SAC9 depletion led to the ectopic localization of PI(4,5)P2 on cortical intracellular compartments, which depends on PI4P and PI(4,5)P2 production at the plasma membrane. SAC9 localizes to a subpopulation of trans-Golgi Network/early endosomes that are enriched in a region close to the cell cortex and that are coated with clathrin. Furthermore, it interacts and colocalizes with Src Homology 3 Domain Protein 2 (SH3P2), a protein involved in endocytic trafficking. In the absence of SAC9, SH3P2 localization is altered and the clathrin-mediated endocytosis rate is reduced. Together, our results highlight the importance of restricting PI(4,5)P2 at the plasma membrane and illustrate that one of the consequences of PI(4,5)P2 misspatterning in plants is to impact the endocytic trafficking.