BRANCHING AND ROOTING OUT WITH A CT SCANNER: THE WHY, THE HOW, AND THE OUTCOMES, PRESENT AND POSSIBLY FUTURE

EDITED BY: Pierre Dutilleul and Jonathan A. Lafond PUBLISHED IN: Frontiers in Plant Science

#### *Frontiers Copyright Statement*

*© Copyright 2007-2016 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

> *All copyright, and all rights therein, are protected by national and international copyright laws.*

*The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-791-0 DOI 10.3389/978-2-88919-791-0

## About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

# What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# **BRANCHING AND ROOTING OUT WITH A CT SCANNER: THE WHY, THE HOW, AND THE OUTCOMES, PRESENT AND POSSIBLY FUTURE**

Topic Editors: **Pierre Dutilleul,** McGill University, Canada **Jonathan A. Lafond,** Université Laval, Canada

Sun rays through the foliage in the forest. Image "Sun Rays" taken from FreeDigitalPhotos.net (http://www.freedigitalphotos.net/images/ sun-rays-photo-p298350).

Until recently, a majority of the applications of X-ray computed tomography (CT) scanning in plant sciences remained descriptive; some included a quantification of the plant materials when the root-soil isolation or branch-leaf separation was satisfactory; and a few involved the modeling of plant biology processes or the assessment of treatment or disease effects on plant biomass and structures during growth. In the last decade, repeated CT scanning of the same plants was reported in an increasing number of studies in which moderate doses of X-rays had been used.

Besides the general objectives of Frontiers in Plant Science research topics, "Branching and Rooting Out with a CT Scanner" was proposed to meet specific objectives: (i) providing a non-technical update on knowledge about the application of CT scanning technology to plants, starting with the type of CT scanning data collected (CT images vs. CT numbers) and their processing in the graphical and numerical approaches; (ii) drawing the limits of the CT scanning approach, which because it is based on material density can distinguish materials with contrasting or moderately overlapping densities (e.g., branches vs. leaves, roots vs. non-organic soils) but not the others (e.g., roots vs. organic soils); (iii) explaining with a sufficient level of detail the main procedures used for graphical, quantitative and statistical analyses of plant CT scanning data, including fractal complexity measures and statistics appropriate for repeated plant CT scanning, in experiments where the research hypotheses are about biological processes such as light interception by canopies, root disease development and plant growth under stress conditions; (iv) comparing plant CT scanning with an alternative technology that applies to plants, such as the phenomics platforms which target leaf canopies; and (v) providing current and potential users of plant CT scanning with up-to-date information and exhaustive documentation, including clear perspectives and well-defined goals for the future, for them to be even more efficient or most efficient from start in their research work.

**Citation:** Dutilleul, P., Lafond, J. A., eds. (2016). Branching and Rooting Out with a CT Scanner: The Why, the How, and the Outcomes, Present and Possibly Future. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-791-0

# Table of Contents


Pierre Dutilleul, Liwen Han, Fernando Valladares and Christian Messier

*20 A Comprehensive Approach to Assess* **Arabidopsis** *Survival Phenotype in Water-Limited Condition Using a Non-invasive High-Throughput Phenomics Platform*

Emilio Vello, Akiko Tomita, Amadou Oury Diallo and Thomas E. Bureau


Craig J. Sturrock, James Woodhall, Matthew Brown, Catherine Walker, Sacha J. Mooney and Rumiana V. Ray


Jonathan A. Lafond, Liwen Han and Pierre Dutilleul

# Editorial: Branching and Rooting Out with a CT Scanner: The Why, the How, and the Outcomes, Present and Possibly Future

Pierre Dutilleul <sup>1</sup> \* and Jonathan A. Lafond<sup>2</sup>

<sup>1</sup> Department of Plant Science, McGill University, Montréal, QC, Canada, <sup>2</sup> Département des sols et de Génie Agroalimentaire, Université Laval, Québec, QC, Canada

Keywords: X-ray computed tomography (CT) scanning, leaf canopies and branching patterns, root systems and other below-ground plant organs, graphical vs. quantitative and statistical analyses, summary and perspectives on plant CT scanning

#### **The Editorial on the Research Topic**

#### **Branching and Rooting Out with a CT Scanner: The Why, the How, and the Outcomes, Present and Possibly Future**

At the time this editorial was written, an extremely interesting loop was closing! Indeed, this text was written last after the invitation received in October 2012 and the proposal submitted in August 2013 for a Frontiers in Plant Science research topic (specialty: Plant Biophysics and Modeling) on "plant CT scanning," the preliminary submission of abstracts in December 2013 and the seven peer-reviewed contributions endorsed for publication and published in March-December 2015. A rapid "guided tour" of the Research Topic, finally titled "Branching and Rooting Out with a CT Scanner" (sub-title: "The Why, the How, and the Outcomes, Present and Possibly Future"), is proposed here. The field of research (excluding synchrotron-based X-ray microtomography of plant tissues) is introduced first, followed by the aims and scope of the Research Topic and a brief description of the major achievements made in each of the articles (2 Methods, 4 Original Research, 1 Perspective). Acknowledgements are made separately.

The use of X-ray computed tomography (CT) scanning with plants is part of the incorporation of new technologies into the scientific research work. Since CT scanners were originally designed for medical applications and later for industrial applications, some adjustment to the reality of the plant world has been necessary, in the air or in the soil or water medium depending on the plant structure or material of interest: the branching pattern in a leaf canopy and the leaves themselves, or a root system. To the best of our knowledge, the first published CT scanning application in which plants were involved is the study of Tollner et al. (1987). The plant part of it consisted in visualizing portions of a root system in the surrounding soil biota, through a few 2-D cross-sectional CT images with no attempt to isolate the entire root system from the soil medium and produce a 3-D image. Plant CT scanning applications to study leaf canopies and relate the complexity of branching patterns to light interception in particular came only in the early 2000s (Dutilleul, 2003), but immediately with 3-D image construction because a canopy is surrounded with air, the challenge coming from the separation of leaves from branches.

Until recently, a majority of the applications of X-ray CT scanning in plant sciences remained descriptive; some included a quantification of plant materials when the root-soil isolation or

Edited and reviewed by: Maciej Andrzej Zwieniecki,

University of California, Davis, USA \*Correspondence: Pierre Dutilleul

#### Specialty section:

pierre.dutilleul@mcgill.ca

This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science

Received: 19 December 2015 Accepted: 11 January 2016 Published: 02 February 2016

#### Citation:

Dutilleul P and Lafond JA (2016) Editorial: Branching and Rooting Out with a CT Scanner: The Why, the How, and the Outcomes, Present and Possibly Future. Front. Plant Sci. 7:41. doi: 10.3389/fpls.2016.00041 branch-leaf separation was satisfactory; and a few involved the modeling of plant biology processes or the assessment of treatment or disease effects on plant biomass and structures during growth. In the last decade, repeated CT scanning of the same plants was reported in an increasing number of studies in which moderate X-ray doses had been used. Besides the general objectives of Frontiers in Plant Science research topics, "Branching and Rooting Out with a CT Scanner" was proposed to meet specific objectives: (i) providing a non-technical update on knowledge about the application of CT scanning technology to plants; (ii) drawing the limits of the CT scanning approach, which is based on material density; (iii) explaining with a sufficient level of detail the main procedures used for graphical, quantitative, and statistical analyses of plant CT scanning data, including fractals and statistics appropriate for repeated plant CT scanning, in experiments where the research hypotheses involve key biological processes; (iv) comparing plant CT scanning with an alternative technology targeting leaf canopies; and (v) providing current and potential users of plant CT scanning with up-todate information and exhaustive documentation, including clear perspectives and well-defined goals for future projects.

Starting with the two Methods articles, Dutilleul et al. and Vello et al. can be used to briefly discuss objectives (ii)–(iv). On the plant CT scanning side (Dutilleul et al.), the crowns of various miniature conifers with contrasting leaf types (needlelike vs. scalelike) were reconstructed in 3D from the CT numbers collected at high resolution, and the crown traits measured from 3-D images (leaf areas and volumes, fractal dimensions of branching patterns) were correlated with a shade tolerance index, the conclusion being that mean values and correlations can differ with leaf type. Regarding the use of phenomics platforms with plants (Vello et al.), alternative non-invasive technologies, exploiting sensors for visible, fluorescent and near-infrared lights, are described to accurately screen survival phenotypes in Arabidopsis thaliana exposed to water-limited conditions. Following two drought protocols and a robust analysis methodology, the authors clearly assessed the plant wilting or dryness status at different time points. In the Original Research articles by Subramanian et al. and Sturrock et al., root system architectures of crop seedlings were investigated, respectively, under salt stress over weeks following seed planting (crop: corn) at high resolution and over days after inoculation of a plant pathogenic fungus (crops: wheat, oil seed rape) at ultra-high resolution, both using repeated-measures ANOVA for statistical analysis because of the repeated plant CT scanning. While Subramanian et al. included fractal dimensions of upper- and

# REFERENCES

lower-root systems of corn seedlings among the measured traits, Sturrock et al. assessed disease severity in relation to pathogen DNA quantification in soil using real-time PCR. Very interestingly, Paya et al. used X-ray CT scanning to uncover root–root interactions and quantify spatial relationships between interacting root systems in 3D. More specifically, they assessed the effects of inter- vs. intra-specific interactions of very young tree seedlings (quaking aspen, black spruce) on their development and utilization of 3-D space by roots, with adjustment of the statistical analysis to the unbalanced design. Very originally, Koebernick et al. combined soil and root water flows by simulations depending on the parameterization and description of the root system (crop: faba bean), and showed how developing root architectures derived from temporally repeated X-ray CT scanning can be implemented in numerical soil-plant models. With their Perspective article. Lafond et al. provide a timely summary of the situation in the X-ray CT scanning of root systems and leaf canopies, data collection and analyses, thus fulfilling objectives (i) and (v) of the Research Topic and helping readers to be efficient in their future work.

# AUTHOR CONTRIBUTIONS

For this Editorial, PD prepared a first draft, which was read and edited by JL before and after a meeting in Montréal, Canada on Wednesday, December 16, 2015. The submitted version is the result of the final editing by PD following JL's comments. Including the Editorial, PD signs/co-signs a total of four contributions, while JL signs/co-signs two. As Topic Editors, PD handled and coordinated the review of three of the other submissions that were endorsed for publication, and JL, one.

# ACKNOWLEDGMENTS

A big Thank You to all 28 contributors, from European and North American institutions, for making this outstanding contribution to Plant Science research possible and to independent reviewers, Review Editors and Associate Editors as well as members of the Frontiers in Plant Science Editorial Office for their valuable help and advice at each step of the process, from Project Manager Carina Paraiso to Assistants (in chronological order) Drs. Gearóid Ó Faoleán, Gina Raihani, and Tobias Preuten, and "Plant Biophysics and Modeling" Specialty Chief Editor, Prof. Maciej Zwieniecki.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Dutilleul and Lafond. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Dutilleul, P. (2003). "Spatio-temporal repeated measures of lengthening curves and expanding shapes and the complexity analysis of plant structures using a CT scanner" in Invited Talk at The 2003 Joint Statistical Meetings, San Francisco, CA, August 3–7.

Tollner, E. W., Verma, B. P., and Cheshire, J. M. (1987). Observing soiltool interactions and soil organisms using X-ray computer tomography. Trans. Am. Soc. Agric. Eng. 30, 1605–1610. doi: 10.13031/2013. 30611

# **Crown traits of coniferous trees and their relation to shade tolerance can differ with leaf type: a biophysical demonstration using computed tomography scanning data**

*Pierre Dutilleul <sup>1</sup> \*, Liwen Han1, Fernando Valladares <sup>2</sup> and Christian Messier 3, 4*

*<sup>1</sup> Environmetrics Laboratory, Department of Plant Science, McGill University, Montréal, QC, Canada, <sup>2</sup> Museo Nacional de Ciencias Naturales, Consejo Superior de Investigaciones Cientificas, Madrid, Spain, <sup>3</sup> Département des sciences biologiques, Centre d'étude de la forêt (CEF), Université du Québec à Montréal, Montréal, QC, Canada, <sup>4</sup> Département des ressources naturelles, Institut des Sciences de la Forêt tempérée (ISFORT), Université du Québec en Outaouais, Ripon, QC, Canada*

#### *Edited by:*

*Sanna Sevanto, Los Alamos National Laboratory, USA*

#### *Reviewed by:*

*Lars Hendrik Wegner, Karlsruhe Institute of Technology, Germany Yusuke Onoda, Kyoto University, Japan*

#### *\*Correspondence:*

*Pierre Dutilleul, Department of Plant Science, McGill University, Macdonald Campus, 21,111 Lakeshore Road, Ste-Anne-de-Bellevue, QC H9X 3V9, Canada pierre.dutilleul@mcgill.ca*

#### *Specialty section:*

*This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science*

> *Received: 25 November 2014 Accepted: 03 March 2015 Published: 24 March 2015*

#### *Citation:*

*Dutilleul P, Han L, Valladares F and Messier C (2015) Crown traits of coniferous trees and their relation to shade tolerance can differ with leaf type: a biophysical demonstration using computed tomography scanning data. Front. Plant Sci. 6:172. doi: 10.3389/fpls.2015.00172* Plant light interception and shade tolerance are intrinsically related in that they involve structural, morphological and physiological adaptations to manage light capture for photosynthetic utilization, in order to sustain survival, development and reproduction. At the scale of small-size trees, crown traits related to structural geometry of branching pattern and space occupancy through phyllotaxis can be accurately evaluated in 3D, using computed tomography (CT) scanning data. We demonstrate this by scrutinizing the crowns of 15 potted miniature conifers of different species or varieties, classified in two groups based on leaf type (10 needlelike, 5 scalelike); we also test whether mean values of crown traits measured from CT scanning data and correlations with a shade tolerance index (STI) differ between groups. Seven crown traits, including fractal dimensions (FD1: smaller scales, FD2: larger scales) and leaf areas, were evaluated for all 15 miniature conifers; an average silhouette-to-total-area ratio was also calculated for each of the 10 needlelike-leaf conifers. Between-group differences in mean values are significant (*P <* 0*.*05) for STI, FD1, FD2, and the average leaf area displayed (*A*¯ D). Between-group differences in sign and strength of correlations are observed. For example, the correlation between STI and FD1 is negative and significant (*P <* 0*.*10) for the needlelike-leaf group, but is positive and significant (*P <* 0*.*05) for the miniature conifers with scalelike leaves, which had lower STI and higher FD1 on average in our study; the positive correlation between STI and *A*¯ <sup>D</sup> is significant (*P <* 0*.*05) for the scalelike-leaf group, and very moderate for the needlelike-leaf one. A contrasting physical attachment of the leaves to branches may explain part of the between-group differences. Our findings open new avenues for the understanding of fundamental plant growth processes; the information gained could be included in a multi-scale approach to tree crown modeling.

**Keywords: plant light interception and shade tolerance, conifer crowns, needlelike vs. scalelike leaves, leaf area and volume, silhouette-to-total-area ratio (STAR), branching pattern complexity, fractal dimensions (FD), computed tomography (CT) scanning**

# **Introduction**

Individual plant light interception has been the subject of many studies in the last decades (see, e.g., Dutilleul et al., 2008; Duursma et al., 2012), because the amount of diffuse radiation intercepted for photosynthetic utilization by individual plants, and vegetation as a whole, plays an important role in growth and development, and hence productivity, as well as in atmospheric CO2 recycling and total carbon uptake (Ackerly and Bazzaz, 1995; Roderick et al., 2001; Sage and Coleman, 2001; Thornley, 2002). From both the agronomic and ecological perspectives, the relationship between light penetration (or its complement, light interception) and leaf area has been modeled with the Beer-Lambert law for light penetration into translucent media (Monsi and Saeki, 1953, 2005); see, e.g., Foroutan-pour et al. (2001), Duursma et al. (2012), and the references therein. A concept, or property of some plant species, related to light interception, is shade tolerance, originally proposed as the "capacity to endure shade" by Shirley (1943). While it is now generally acknowledged that shade tolerance indicates the degree to which a plant can survive and grow in low light conditions (Kobe et al., 1995), the survival, development, and reproduction of a plant species at a particular light level does not mean or imply that this species is at its physiological optimum (Humbert et al., 2007).

Crown traits important for plant light interception efficiency may be the same that influence shade tolerance. Among those traits, some can characterize the geometric structure and complexity of the branching pattern (fractal dimensions), and others, the amount of leaves (volume in 3D, number, areas in 2D). That is why, focusing on trees and coniferous species in particular, one of the objectives of our study was to investigate the existence of links and assess the strength and sign of correlations between various crown traits and a shade tolerance index, taking into account that there is no "classical" attachment of leaf blades to branches via petioles in conifers and leaf type for them can be of two types: needlelike or scalelike. In a detailed discussion of the architecture of terrestrial plants and their modular nature, the structural determinants of light capture and the 2-D and 3-D geometries of foliage arrangement within the crown, among other themes, Valladares and Niinemets (2007) emphasize that conifers have extensively aggregated foliage, citing Oker-Blom and Smolander (1988), Niinemets (1997), and Stenberg et al. (2001) for this.

Working at the whole-tree scale while collecting sufficiently fine data to measure crown traits accurately and thoroughly represents a challenge that we have addressed by applying highresolution X-ray computed tomography (CT) scanning to the above-ground structure of miniature conifers, of less than 30 cm in width and height and growing in pots. Indeed, a CT scanner of medical type, like the one at the CT Scanning Laboratory for agricultural and environmental research on Macdonald Campus of McGill University, can be used for such small-size trees, and make indirect measurement of density on parallelepipedal rectangular parts (called voxels, the extension of pixels in 2D) of a size as small as 0.23×0*.*23×0*.*20 mm3 (width: 12 cm; height: 10 cm). Thus, the criterion of high resolution of Ketcham and Carlson (2001, Table 1) is satisfied, and more than 100 million data points (CT numbers) can be obtained for each individual tree. Recent applications of CT scanning technology in the plant sciences have been investigating the structural complexity of root systems (see, e.g., Gregory et al., 2003; Anderson and Hopmans, 2013), more than that of canopies or crowns (Dutilleul et al., 2005, 2008), and the fundamental question of tree growth using stem sections (see Dutilleul et al., 2014 and references therein).

Our objectives with the study reported here were multiple. Basically, we wanted to upgrade (analytically speaking) and expand (botanically speaking) the work of Dutilleul et al. (2008), who studied the developing crowns (branching patterns and leaf canopies) of four pyramidal (non-miniature) cedars (*Thuja occidentalis* Fastigiata). Leaving developing crown aspects aside but including plant light interception efficiency when possible, we aimed to (i) compare crown traits measured from CT scanning data between miniature conifers with needlelike vs. scalelike leaves, covering as wide a range of shade tolerance as possible; (ii) test for statistical differences between mean values of the crown traits and a shade tolerance index depending on leaf type; (iii) analyze correlations between crown traits and shade tolerance; and (iv) discuss possible biological meanings of (ii) and (iii). So doing, we tested the hypothesis of multi-fractality of branching pattern (Stewart, 1988; Prusinkiewicz and Lindenmayer, 1990), and provided supplementary data for the bottom right part of Figure 7 in Duursma et al. (2012), relating the average silhouetteto-total-area ratio (STAR) to the number of leaves; this number (*N*) is known to be large in conifers. All the abbreviations of variables that we have analyzed are defined, with their unit when they have one, in **Table 1**.

# **Materials and Methods**

### **Plant Material**

The scientific and common species names of the 15 miniature conifers that we studied are listed in **Table 2**, where they are classified according to their leaf type (i.e., 10 needlelike, 5




**TABLE 2 | Scientific and common species names of the 15 miniature conifers, their leaf type and shade tolerance index value, and the estimates and associated standard errors of the two fractal dimensions of their branching pattern.**

\**Depicted in Figure 1.*

scalelike; see **Figure 1** for examples). These conifers, which are horticultural varieties of indigenous species, were grown in the Canadian Province of British Columbia (Pacific Northwest Propagators Inc., Rosedale) and purchased through a Montréal (Québec) plant shop in the month of May, prior to the scanning of their crowns in June-July (see below). Their height and width then varied, respectively, from 7.9 to 22.2 cm and from 10.3 to 27.5 cm, respectively. Trees were about 2 years old at the time of purchase. The rooting medium was organic, fertilized with Polyon NPK 16-3-13 plus minors.

### **Shade Tolerance Index**

As mentioned in the Introduction, plant tolerance to low light conditions, for survival and development, must be distinguished from optimum sun exposure, as recommended in horticulture for example. Accordingly, the shade tolerance index values for the varieties studied in **Table 2** are equal or close to the values reported for indigenous species in forest ecology data sources (e.g., *Sylvics of North America*, Burns and Honkala, 1990; see also Kwantlen Polytechnic University, 2012). Values of our index range from 1.0 (requires Full Sun) to 5.0 (well-adapted to Full Shade), by increments of 0.5; the median value of 3.0 thus corresponds to: needs Half Sun/accommodates Half Shade.

#### **Computed Tomography Scanning**

The 15 crowns of miniature conifers were CT scanned at the CT Scanning Laboratory for agricultural and environmental research on Macdonald Campus of McGill University in Ste-Anne-de-Bellevue (Québec, Canada), which is equipped with a helical high-resolution CT scanner XVision (Toshiba Corporation, Medical Systems Division, Tokyo, Japan). CT scanning sessions were distributed over several days, with two or three trees CT scanned per day. Basic CT scanning configuration parameters were the same for all 15 trees: tube current, 50 mA, and tube

voltage, 120 kV. The X-ray beam width (i.e., a parameter to be distinguished from the thickness of CT images) was the same (1 mm) for all but one (the widest tree, for which 2 mm was used).

Depending on the width of the crown, a different field of view was used: SS, very small (18 cm in diameter); S, small (24); or M, medium (32). A zoom factor was used to improve spatial resolution horizontally; for example, for the crown of *Chamaecyparis pisifera* Golden Pin Cushion (Sawara Falsecypress), which had a width of 10.3 cm, a zoom factor of 1.5 was used with the SS field of view. Among the 11 miniature conifers that were CT scanned with the SS field of view, a zoom factor (ZF) was used for six (two times with *ZF* = 1*.*5 and four times with *ZF* = 1*.*2). Vertically, the thickness of CT images was 0.2, 0.3, or 0.4 mm, depending on the height of the tree; its width was also taken into account in order to have the same resolution as much as possible in all three dimensions (e.g., 0.23×0.23×0.20 mm<sup>3</sup> for Sawara Falsecypress). Between 400 and 600 CT images, each made of 512 × 512 CT numbers, were produced. Of all the CT images produced, a small portion (the top ones) corresponded to pure air and another small portion (the bottom ones) contained information about the base of the tree and surface soil and roots.

Prior to CT scanning, the equipment was calibrated with the appropriate phantoms, so that the CT numbers for air and water corresponded to −1000 and 0 Hounsfield units (HU), respectively. Following CT scanning, the raw data files, which contained between ca. 100 and 150 million CT numbers in 512 × 512 matrices, were transferred to a Windows 7 Dell workstation for graphical and numerical analyses in MATLAB R2014a (The MathWorks Inc.).

#### **Fractal Dimension Estimation**

For reasons already made clear when complexity of the aboveground structure was analyzed from photographs of plants from which leaves had been removed manually (i.e., the thickness of branches introduces a bias; Foroutan-pour et al., 1999a), fractal dimension estimation in our study was performed on skeletal branching patterns, prepared in a customized MATLAB graphical unit interface by tracing branches using the 3-D array of CT scanning data collected for the crown of a miniature conifer. The fact that branching patterns were skeletal means that their thickness was one voxel (i.e., the 3-D extension of one pixel in 2D); in our case, a voxel has two opposite square faces, and is the smallest volumetric unit for which a CT number is produced.

A cube-counting procedure was used to estimate fractal dimensions; for reasons that will appear clearly in the Results, two fractal dimensions (denoted FD1, FD2) were estimated for each tree. In the cube-counting procedure for fractal dimension estimation of an object or a structure in 3-D space, the object or structure of interest (e.g., a skeletal branching pattern) first needs to be included in the smallest cube that can contain it. In our framework, the length of the sides of that cube is given by the larger of two quantities: the number of horizontal sections containing the 3-D image of the skeletal branching pattern and the diameter (in voxels) of the smallest circle containing its vertical projection (along the Z-axis, onto the X-Y plane). For the two examples which will be detailed, that length is a perfect 400 (which can easily be divided by powers of 2, without rounding) and 438 (resulting in 219, 109.5, 54.75, 27.375, 13.6875, 6.84375, 3.421875, 1.7109375 after successive divisions by 2, rounded to 219, 110, 55, 27, 14, 7, 3, 2). These nine decreasing cube sidelengths, denoted "*s*" hereafter, provide as many scales, larger or smaller. The number of cubes with sidelength *s* intersecting a skeletal branching pattern was thus counted for 9 scales; **Figure 2** illustrates the cube-counting procedure for the three cube sidelengths or scales corresponding to the divisions by 2, 4, and 8. In our customized MATLAB program, cube counting for a given scale was repeated 8 more times, by moving the "smallest cube containing the entire structure of interest" one voxel on the left/right, or one voxel in the front/back, or the two together. The minimum count for cube sidelength *s* (denoted "*C*(*s*)" below) was retained for a more accurate estimation, in accordance with similar procedures in 2D (Foroutan-pour et al., 1999b) and 3D (Lontoc-Roy et al., 2006); see also Li et al. (2009). A fractal dimension estimate is then the estimated slope of the straight line fitted by least squares in the biplot of log(*C*(*s*)) against log(1/*s*) over a number of scales, where log(.) is the natural logarithm and *C*(*s*) denotes the number of cubes with sidelength *s* intersecting the skeletal branching pattern:

$$
\log(\text{C(s)}) = k + \text{FD}\,\log(1/s) \tag{1}
$$

As we shall see in Section Bifractality of Conifer Branching Patterns, the *R*2-values in these fittings, depending on the range of scales covered, are very important.

#### **Plant Light Interception Efficiency**

Besides the shade tolerance index (STI, obtained independently) and the two fractal dimensions (FD1, over smaller scales; FD2, over larger scales), five crown traits were evaluated from the raw CT scanning data or the derived 2-D or 3-D images, for all the 15 miniature conifers. These five crown traits are: absolute leaf area (LA, in mm2) in the vertical projection of the crown; relative leaf area (PCT, in %) in the smallest disc including the vertical projection of the crown; total leaf volume (*V*L, in mm3), obtained from all the leaf voxels (i.e., voxels with a CT number smaller than the branch threshold, and greater than −980 HU to discriminate them from air; see **Figure 3** for examples); leaf volume-to-branch

**FIGURE 3 | (A)** From left to right, one of the 428 gray-tone CT images (thickness: 0.4 mm) constructed for the crown of the miniature coniferous tree (*Picea abies* Tompa) in **Figure 1A**; the branch part delineated from the corresponding 512 × 512 matrix of CT numbers, colored in brown; and the leaf part obtained by subtraction (excluding air) and colored in green. **(B)**

Similar information for the miniature coniferous tree (*Microbiota decussata* Gold Spot) depicted in part in **Figure 1B**; the gray-tone CT image displayed is one out of 450 with thickness of 0.3 mm used to analyze the crown in this case. Note: The two skeletal branching patterns are contained in 400 and 438 CT images, respectively.

volume ratio (*V*L/*V*B); and average leaf area displayed (*A*¯ D, in mm2; Pearcy et al., 2011; Duursma et al., 2012), over four classes of azimuth (the X- and Y-axes in the CT scanning framework and their two directions) × 20 solar elevation classes (every 4.5 degrees starting from 0), plus the 90-degree solar elevation.

For each of the 10 miniature conifers with needlelike leaves, we applied a cylindrical model (height, *h*; radius, *r*) to 25 leaves sampled in the crown on computer, using a random stratified design (1 stratum = 20% of the height of the tree, from bottom to top; 5 leaves randomly sampled and measured per stratum) and the 3-D array of CT numbers and the 3-D image of leaf voxels. This application consisted in equaling the measured individual leaf volume from CT scanning data with the volume of a cylinder, π *r*<sup>2</sup> *h*, using the largest distance calculated in MATLAB between two voxels of the leaf for *h*, and solving the equality for *r*, which then allowed the calculation of an individual leaf area, <sup>π</sup>*r*<sup>2</sup> <sup>+</sup> <sup>2</sup><sup>π</sup> *r h*, excluding the bottom area of the cylinder (where the needlelike leaf is attached to the branch). Dividing the total leaf volume by the mean of the 25 individual leaf volumes thus provided an estimate of the total tree leaf number, *N*, which multiplied by the mean of the 25 individual leaf areas provided in turn an estimate of the total tree leaf area, *A*<sup>L</sup> (in mm2). Finally, an average silhouette-to-total area ratio, STAR (mm2 mm−2) was calculated by dividing *A*¯ <sup>D</sup> by *A*<sup>L</sup> (Oker-Blom and Smolander, 1988; Duursma et al., 2012).

No standardization was applied to *N* in Figure 7A of Duursma et al. (2012), for the investigation of a negative relationship with STAR. Nevertheless, it is a natural question to ask whether a standardization, e.g., by expressing size trait data (LA, *V*L, *N*, *A*¯ D, *A*L) per unit stem basal area (BA), would change correlations. For our miniature conifers, we thus measured BA at a height of about 2.5 cm, starting from the first CT image in the 3-D skeleton of the branching pattern, as an equivalent to breast height in non-miniature trees; as specified in the legend of **Figure 4**, the height of our 15 miniature conifers ranges from 7.9 to 22.2 cm. The number of stem voxels found in the CT image at ca. 2.5-cm height, multiplied by the horizontal area of a voxel in mm2, provided the BA measure. The calculation of BA was made in the same way whether the miniature conifer had one or several stems; three of the 15 miniature conifers had several stems.

#### **Statistical Analyses**

Normality of the distribution of sample data was tested per leaf-type group for a given variable. It was accepted at 5% after arcsine-square-root transformation for PCT, after logtransformation for BA and without transformation for the other variables for which the mean values were compared statistically between leaf-type groups. Accordingly, these comparisons were carried with parametric *t*-tests, using effective degrees of freedom when the sample variances could not be pooled. Spearman's

rank-based coefficient was used for correlation analyses because it can capture non-linear relationships, i.e., it is not restricted to linear relationships like Pearson's sample correlation coefficient. A *t*-test for paired observations was performed to compare the mean values of FD1 and FD2 over all the 15 miniature conifers and per leaf-type group. The SAS 9.3 (SAS Institute Inc.) procedures UNIVARIATE (option NORMAL), TTEST and CORR (option SPEARMAN) were used for the normality tests, *t*-tests for comparisons of means and correlation analyses, respectively.

# **Results**

### **Conifer Crown Image Processing**

The 15 skeletal branching patterns, which were traced from CT scanning data first, are displayed in 2D (front view) in **Figure 4**. Their actual 3-D structure can be better viewed with a customized MATLAB graphical unit interface, but differences in structural complexity among some of the crowns can already be anticipated prior to any quantification with fractal analysis (see Subsection Bifractality of Conifer Branching Patterns). After the appropriate number of layers was added to the branch skeletons, using −650 HU as threshold for most trees to delineate (from the CT numbers) branches from leaves attached to them, entire "digital branches" were obtained (see middle images in the detailed examples of **Figure 3**) and thereafter, by subtraction (excluding air) remained leaves (see right images in the detailed examples of **Figure 3**). This way of proceeding at the whole-tree scale provided the complete crown renderings of **Figure 5**, with a semi-transparency option to allow the eye to penetrate the leaf canopies. Again, prior to any quantification through leaf areas and volumes (see Subsections Differences in Means of Shade Tolerance Index and Conifer Crown Traits and Light Interception Efficiency for Needlelike-Leaf Group), differences among some of the crowns can be anticipated in light interception efficiency.

### **Bifractality of Conifer Branching Patterns**

As a preliminary note, it is important to emphasize that in loglog plots such as those of **Figure 6**, where log(*C*(*s*)) is plotted against log(1/*s*), the smaller scales *s* (i.e., smaller cubeside lengths, after divisions by the greater powers of 2 in the cubecounting procedure) are represented by data points on the right, and the larger scales *s* (i.e., larger cubeside lengths, after divisions by a few powers of 2 in the procedure), by data points on the left. This positioning (smaller scales right, larger scales left) follows from the use of the inverse in log(1/*s*) and the fact that logarithmic functions take negative values for positive quantities smaller than 1.0. Accordingly, we number the scales from right to left (i.e., from smaller scales to larger scales) below.

Following Foroutan-pour et al. (1999b), it is not recommended to include the smallest scales (1 and 2 here; see, e.g., the top-right data points in **Figure 6**) and the largest scales (8 and 9; see, e.g., the bottom-left data points in **Figure 6**) in a boxcounting procedure of fractal dimension estimation, because the

**FIGURE 5 | False-colored 3-D renderings of the complete crowns (branches in brown, leaves in green) for the 15 miniature conifers of Figure 4, in same order and size representation. (A)** with needlelike leaves (first two rows) and **(B)** with scalelike leaves (third row).

estimation would be biased if they were included. Remain then the options of using all five middle data points (scales 3-4-5-6-7) and subsets of four and three successive data points (scales 3-4-5- 6, 4-5-6-7 and 3-4-5, 4-5-6, 5-6-7). We tried them all and found that the FD estimates obtained using Equation (1) with five scales and four scales had intermediate values, between those for scales 3-4-5 (smaller scales) and 5-6-7 (larger scales), and were close to the FD estimates for scales 4-5-6. More concretely, for the examples of **Figure 6**, the FD estimates (with the associated *<sup>R</sup>*2-value as measure of goodness-of-fit in parentheses) read as follows (in the same order of subsets of scales as above): 1.2749 (0.9920), 1.1744 (0.9962), 1.3729 (0.9952), 1.0934 (0.9987), 1.2611 (0.9984), 1.4709 (0.9975) in **Figure 6A** and 1.8218 (0.9923), 1.7299 (0.9883), 1.9739 (0.9988), 1.5176 (0.9925), 1.946 (0.9971), 2.0507 (0.9996) in **Figure 6B**.

Because (i) there is a change in direction when following the 5 middle data points and passing by scale 5 in the log-log plots of **Figure 6**, and (ii) scales 3-4-5 and 5-6-7 provided the FD estimates with the highest *R*2-values on average over the 15 miniature conifers, these were chosen for structural complexity analysis and further statistical inference, in which they were respectively denoted FD1 (smaller scales) and FD2 (larger scales); the use of adjusted *R*2-values (adjusted for the number of data points) does not change this reasoning. Using the FD1 and FD2 values listed in **Table 2** (right columns), we found that the difference between sample means of FD2 and FD1 was very similar in the two groups of coniferous trees: 0.6200 (± 0.0491, *n* = 10) for needlelike-leaf and 0.6290 (± 0.0729, *n* = 5) for scalelike-leaf, and the difference from 0.0 was significantly different at 1% in each group. The last result indicates a greater measured structural complexity for conifer branching patterns at larger scales than at smaller scales.

# **Differences in Means of Shade Tolerance Index and Conifer Crown Traits**

The sample means and associated standard errors of the 12 variables that were studied for the two groups of coniferous trees in relation to their leaf type are presented in **Table 3**, together with the results of *t*-tests for the statistical comparisons of means. Differences are significant at 5% for STI (greater mean for the needlelike-leaf group), FD1 and FD2 as measures of structural complexity of the branching pattern (greater mean for the scalelike-leaf group), and *A*¯ <sup>D</sup> prior to standardization by BA, as one of the important crown traits for plant light interception efficiency (greater mean for the needlelike-leaf group). That is, for one third of the variables studied for both groups of conifers, and one or two variables per type of variable. The absence of a significant difference between the two mean values of *A*¯ <sup>D</sup> after standardization is one of a small number of effects of the standardization that we have observed in our study.

### **Relations to Shade Tolerance and Among Crown Traits**

Several of the significant correlations observed were expected, such as the positive ones: (i) between FD1 and FD2 (when structural complexity of branching pattern is higher/lower over smaller scales, it is higher/lower over larger scales), (ii) between *V*<sup>L</sup> and *V*L/*V*<sup>B</sup> (by construction), and (iii) between a size trait and itself, not standardized vs. standardized by BA (assuming or anticipating that the distribution of BA values among the 15 miniature conifers would be almost uniform), and (iv) the negative correlation between STAR and *A*<sup>L</sup> (by construction). Besides those correlations, it is worth commenting on the following relationships found: (i) correlations between STI and

FD1, FD2 are both significant and negative for the needlelike-leaf group, but both significant and positive for the scalelike-leaf group; (ii) correlation between STI and *A*¯ <sup>D</sup> is not significant for the first group, but is positive and significant at 5% prior to standardization by BA for the second group; and (iii) several correlations between branching complexity measures FD1, FD2, and leaf areas (average displayed or in vertical projection of the crown, absolutely or relatively) are significant and positive for the scalelike-leaf group, whereas it is rather with leaf and branch volumetric measures that FD1 and FD2 are correlated significantly and negatively for the needlelike-leaf group. All the significant correlations, at 5 or 10%, are clearly identified in **Tables 4**, **<sup>5</sup>**. The correlations of *<sup>A</sup>*¯ <sup>D</sup> with STI and FD1, which were both positive and significant at 5% prior to standardization of *A*¯ <sup>D</sup> by BA, remained positive (0.5270 and 0.4000, respectively) but lost their statistical significance after standardization of *A*¯ <sup>D</sup> by BA (see **Table 5B**); the sample size of the scalelike-leaf group (*<sup>n</sup>* <sup>=</sup> 5) makes the discussion of this result difficult to pursue, except to say that the BA measures were not uniformly distributed in that group.

# **Light Interception Efficiency for Needlelike-Leaf Group**

The variables total tree leaf number, total tree leaf area and average silhouette-to-total-area ratio were studied for the 10 individuals of the needlelike-leaf group, because individual leaves could be isolated from CT scanning data for them. The corresponding sample means and standard errors are: *N* (not standardized), 3857 ± 724; *N* (standardized), 76 ± 12; *AL* (not standardized), 253,106 ± 38,810; *A*<sup>L</sup> (standardized), 4884 ± 502; and STAR, 20.53 ± 2.42.

Significant correlations of STAR with other variables, at 5 or 10%, are all negative: regardless of standardization, with PCT; before and after standardization by BA, with *V*L, *N* and *A*L. The scattergram of STAR against *N* (not standardized), using the 10 data points obtained in our study, shows indeed a strong negative relationship (**Figure 7**); this could be expected (Duursma et al., 2012, Figure 7A), but had not yet been observed with numbers of leaves in the thousands.

# **Discussion**

### **CT Scanning Technology vs. Other Approaches to 3-D Tree Crown Reconstruction**

The CT scanning of the crown of any of the 15 miniature conifers in our study generated around 125 million CT numbers (Subsection Computed Tomography Scanning). This numerical CT scanning data, which is made of indirect measures of density of all the parts of the crown (stem, branches, and leaves) and the surrounding air, can be mapped in 512×512 CT images with lighter and darker gray tones for pixels with higher and lower densities (see examples in **Figures 3A,B**, left panels). Even more interestingly, 3-D images of complete crowns can be constructed with branches colored in brown and leaves in green (**Figure 5**); skeletons of branching patterns (**Figure 4**) are first extracted from the CT images and then let "grow" in an iterative procedure, only to draw the limits between the end of a branch and the beginning of a leaf or an area with leaves, for which CT numbers are available too. After CT scanning, we did not proceed to destructive sampling, which would have consisted in detaching the leaves from the branches as in Foroutan-pour et al. (1999a) who studied soybean canopies without a CT scanner. In our case, there is an intrinsic difficulty of defoliating conifers with sufficient accuracy, especially those with scalelike leaves.

Other approaches, procedures and techniques have been used for 3-D crown reconstruction, but for deciduous trees (e.g., hybrid poplar, sugar maple, yellow birch), generally potted 2–3 year-old saplings, and a cactus-like euphorbia tree. For example, Delagrange and Rochon (2011) worked with a hybrid poplar clone grown in a nursery for 2 years, to compare the results obtained from high-definition photographs in the Tree Analyser (TA) software using a space carving approach vs. 3-D "point clouds" acquired from terrestrial light detection and ranging (T-LiDAR) scans performed on trees without leaves to reconstruct the lignified structure of the sapling, to which foliage was added using allometric relationships between the number


**TABLE 3 | Group means depending on leaf type and corresponding standard errors for the shade tolerance index and the crown traits that were measured from CT scanning data and CT images for both groups of miniature conifers, together with the result of the statistical comparison of the two means per variable (AH0, Accept the null hypothesis of equality of means; RH0, Reject, at 5%).**

*For the definition of abbreviations, see Table 1. †Not standardized, ‡Standardized by basal area.*

of leaves and maximum leaf length and the length of the current-year shoot. Even though some discrepancy is visible between the crown in the black-and-white picture in Figure 1A of Delagrange and Rochon (2011) and the crown reconstructed from T-LiDAR scans in their Figure 2D, the authors found that T-LiDAR is better than TA in terms of precision and accuracy of the reconstruction. The discrepancy is likely to be the result of performing T-LiDAR scanning of the crown after leaves had been detached from it; this could have implications for the results of fractal analysis (structural complexity) and light interception efficient analysis (space occupancy). Using portable scanning LiDAR data, Hosoi et al. (2013) developed a method to produce a 3-D voxel-based solid model of a tree (voxel size: 0.5×0.5×0.5 cm3), for accurate estimation of the volume of woody material. They applied their model to a Japanese zelkova tree, and found very satisfactory results for the stem and large branches (diameter *>* 1 cm); the error in volume estimates was 0.5%.

Investigating an approach for 3-D data collection on plant architecture that would not be time-consuming and would not require costly experiment, Nock et al. (2013) tested a low-cost, 3- D camera and open-source software for the measurement of stem and branch diameters and lengths. Besides technical and calibration aspects, the authors report that the tested camera is able to accurately capture the diameter of maple branches *>* 6 mm, and a cactus-like euphorbia is well-acquired thanks to the width of its axes. Focusing on tree crown reconstruction from point clouds acquired with terrestrial LiDAR scanning, the study of Delagrange et al. (2014) can be seen as a follow-up to Delagrange and Rochon (2011), with the presentation of an improved skeletal extraction method for use with 104 or 105 points in the cloud for a synthetic tree and real branches of elm and from 3-year-old sugar maple and yellow birch saplings grown under low and high light regimes; some difficulty to detect smaller branches (length *<* 3.5 cm) could potentially affect the results of a fractal analysis of the branching pattern, whereas the authors declare that small branches account for little in terms of total skeleton length; no results for leaves were reported.

#### **Space Occupancy in Crowns of Coniferous Trees**

The skeletal branching patterns constructed from CT scanning data in this study revealed "gaps" in the crowns of coniferous trees that are partially filled to various degrees, depending on species, by needlelike or scalelike leaves (**Figures 4**, **5**). The size of gaps in skeletal branching patterns is related to the level of structural complexity: the more complex the branching patterns (the higher the degree of subdivision of branches), the smaller the gaps, and vice versa; see "Foliage Dispersion" in Valladares and Niinemets (2007, p. 119). We observed this relationship for needlelike-leaf coniferous trees, but not directly on total tree leaf volume. We observed this relationship for that group through a significant and negative correlation between FD1, FD2 and leaf volume-to-branch volume ratio. We did not observe a similar relationship for the scalelike-leaf group, which may be explained by the different type of leaves (i.e., short and pasted on branches vs. longer needles forming a certain angle with the branch) and a different range of FD1, FD2 values (i.e., lower for needlelike-leaf vs. higher for scalelike-leaf; **Tables 2**, **3**). The higher (lower) the value of the fractal dimension parameter, the more (less) complex the structure of the branching pattern.

### **Range of Shade Tolerance Index Values and the Leaf-Type Classification Factor**

Our primary goal in conducting this study was not to cover the whole range (1–5) of possible values for the shade tolerance index in each of two groups of miniature conifers classified according to their leaf type. Instead, it was of biophysical nature—the 15 experimental trees were chosen based on the architecture of their crown, among the species and varieties available at the grower at


**rank-basedcorrelationcoefficients(withprobabilitiesofsignificancebelow,inparentheses)betweentheshadetoleranceindexandtraitsmeasuredfromCT**

Frontiers in Plant Science | www.frontiersin.org March 2015 | Volume 6 | Article 172 |

**TABLE 5 | Spearman's rank-based correlation coefficients (with probabilities of significance below, in parentheses) between the shade tolerance index and crown traits measured from CT scanning data and CT images appropriately processed, by group depending on leaf type (***a***) needlelike (***<sup>n</sup>* **<sup>=</sup> <sup>10</sup>) and (***b***) scalelike (***<sup>n</sup>* **<sup>=</sup> <sup>5</sup>); the correlations statistically significant at 5% are bolded and underlined; those only significant at 10% are simply underlined.**


*For the definition of abbreviations, please see Table 1. †Not standardized, ‡Standardized by basal area.*

the time. The main research objective was to address technological challenges, including the collection and advanced appropriate processing of CT scanning data for tree crown reconstruction at an unprecedented level, while providing insight on the physiological side—through shade tolerance—and thus preparing for future larger studies on the subject after establishment of the technological and analytical protocol and procedures. Accordingly, it is important to keep in mind the ranges of shade tolerance index values in the present study, when commenting the differences observed between groups in mean values of other crown traits and in the correlations between some of the crown traits in the last subsection of the Discussion. The inclusion of miniature varieties of *Pinus* in the needlelike-leaf group would definitely decrease the mean value of the shade tolerance index value of that group; a presently open question that can be answered in a later study concerns the resulting effects on mean values of other crown traits and related correlations for that group.

## **Leaf Type Effects and Observed Differences in Shade Tolerance, Crown Structure and Light Interception**

Based on the relative 2-D and 3-D measures that are the proportion of leaf area in the vertical projection of the crown and the leaf volume-to-branch volume ratio, the two groups of miniature conifers showed very little difference in mean values, while the corresponding absolute measures showed a greater difference which was not close to be significant (**Table 3**). Interestingly, the four variables that show significant differences in mean values between groups also show differences in correlations (**Tables 3**, **4**).

In contrast with the rigidity and proximity of scalelike leaves to branches, needles have a well-identified point of attachment to branches in the form of an alveolus and possess a greater mobility around their point of attachment; this cannot explain *per se* the higher mean value of the index for the needlelike-leaf group (see

the point raised in Subsection Range of Shade Tolerance Index Values and the Leaf-Type Classification Factor), but may provide a leaf canopy more instrumental in capturing the irradiance in low light environments. The size of leaves seems larger for needles, but we have only visual observations from CT images to use for scalelike leaves.

The differences in mean values of FD1, FD2 indicate that: (i) branches tend to be aligned at smaller scales in 3D (i.e., the mean FD1 is close to 1.0, the dimension of a straight line in classical Euclidean geometry) in the needlelike-leaf group, but are laid in a manner mid-way between linear and planar (i.e., the mean FD1 is close to 1.5, between 1.0 and 2.0, the dimension of a plane) in the scalelike-leaf group; and (ii) over larger scales, the spatial distribution of branches is almost planar (i.e., the mean FD2 is close, but not equal to 2.0) for our miniature conifers with needlelike leaves, and slightly more complex than a planar distribution (i.e., the mean FD2 is greater than 2.0) for those with scalelike leaves. These results concern the crown structures, and reflect a greater flexibility and degree of repeated subdivision of the branching patterns of scalelike-leaf trees, likely to adjust for the characteristics of their leaves and succeed in capturing a sufficient amount of light for their survival and development. The difference in mean values of the average leaf area displayed is consistent with the other differences observed in mean values between groups, and could be expected in some way; the mean value of *A*¯ <sup>D</sup> is greater for the needlelike-leaf group, that is, the reverse than for FD1 and FD2.

Combining our results on the differences in mean values of the shade tolerance index and FD1 with those on the change in sign of their correlations depending on leaf type (i.e., negative for needlelike and positive for scalelike), a relationship between the two variables that is quadratic and asymmetric, instead of linear or at the least monotonic, seems possible. In fact, while the shade tolerance index takes its highest value of 4.5 for FD1 *<* 1.1 (two times) in the needlelike-leaf group (minimum FD1: 1.0639, maximum FD1: 1.3465), it takes its highest value of 3.5 for FD1 = 1.5781 after a minimum of 1.5 (two times) for FD1 = 1.2741 and 1.4252 in the scalelike-leaf group (**Table 2**). A larger sample size for the scalelike-leaf group, including a broader range of shade tolerance index values in that group (Subsection Range of Shade Tolerance Index Values and the Leaf-Type Classification Factor), would allow the assessment of such a quadratic relationship more finely.

# **Concluding Remarks**

So far, a high-resolution CT scanner (non-micro, due to the size of the field of view) has been used in a much smaller number of studies on plant canopy architecture, compared to those conducted on plant root systems using micro or non-micro CT scanning. One of our conclusions is that there could or should be more studies like ours (i.e., for crowns of small-size trees), because the resolution achieved with a X-ray CT scanner of medical type such as a Toshiba XVision is very fine and sufficient at the whole-plant scale (i.e., down to 0.23 <sup>×</sup> 0.23 <sup>×</sup> 0.20 mm3 and not exceeding 0.62 <sup>×</sup> <sup>0</sup>*.*<sup>62</sup> <sup>×</sup> <sup>0</sup>*.*4 mm3 in our study). Accordingly, detailed graphical and quantitative information could be gathered for two groups of miniature conifers with different leaf types (i.e., needlelike vs. scalelike), regarding their leaf areas and volumes and the complexity of their branching patterns.

Differences between groups in mean values of crown traits measured from CT scanning data and a shade tolerance index obtained separately were assessed statistically. Significant differences were found for shade tolerance, fractal dimensions and the average leaf area displayed. These differences between mean values had implications for correlations; in particular, shade tolerance was negatively correlated with fractal dimensions for the needlelike-leaf group, and positively correlated with one fractal dimension in the miniature coniferous with scalelike leaves studied. These findings were complemented with the acceptance of the hypothesis of bifractality of the branching pattern over the two groups of miniature conifers and the presentation of new documentation for conifers with needlelike leaves about the strong negative relationship between the average silhouette-tototal-area ratio and the number of leaves, when the latter are in large to very large numbers. In closing, our results here, obtained for crowns of miniature conifers analyzed thoroughly and accurately thanks to CT scanning technology and advanced data processing, could be used for crown modeling of non-miniature indigenous species in situations where the leaf size-to-branch length ratio would justify it, for example for juvenile indigenous trees of a size that fits in the gantry of the CT scanner.

# **Acknowledgments**

The authors are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for the Major Equipment Grant that led to the creation of the CT Scanning Laboratory for agricultural and environmental research on Macdonald Campus of McGill University and to the Canadian federal government and the Québec provincial government for a Canada Foundation for Innovation (CFI) Grant that provided us with the computer equipment. The operating aspects of the reported study were supported by NSERC through a Discovery Grant to PD under Grant No. Nserc Rgpin 138514-10.

# **References**


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Dutilleul, Han, Valladares and Messier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# A Comprehensive Approach to Assess *Arabidopsis* Survival Phenotype in Water-Limited Condition Using a Non-invasive High-Throughput Phenomics Platform

#### *Emilio Vello\*, Akiko Tomita, Amadou Oury Diallo and Thomas E. Bureau\**

*Department of Biology, McGill University, Montreal, QC, Canada*

#### *Edited by:*

*Jonathan A. Lafond, Université Laval, Canada*

#### *Reviewed by:*

*Bettina Berger, University of Adelaide, Australia Christophe Salon, Institut National de la Recherche Agronomique, France*

#### *\*Correspondence:*

*Emilio Vello emilio.vello@mcgill.ca; Thomas E. Bureau thomas.bureau@mcgill.ca*

#### *Specialty section:*

*This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science*

*Received: 12 September 2015 Accepted: 22 November 2015 Published: 15 December 2015*

#### *Citation:*

*Vello E, Tomita A, Diallo AO and Bureau TE (2015) A Comprehensive Approach to Assess Arabidopsis Survival Phenotype in Water-Limited Condition Using a Non-invasive High-Throughput Phenomics Platform. Front. Plant Sci. 6:1101. doi: 10.3389/fpls.2015.01101*

With the rapid rise in global population and the challenges caused by climate changes, the maximization of plant productivity and the development of sustainable agriculture strategies are vital for food security. One of the resources more affected in this new environment will be the limitation of water. In this study, we describe the use of non-invasive technologies exploiting sensors for visible, fluorescent, and near-infrared lights to accurately screen survival phenotypes in *Arabidopsis thaliana* exposed to water-limited conditions. We implemented two drought protocols and a robust analysis methodology that enabled us to clearly assess the wilting or dryness status of the plants at different time points using a phenomics platform. In conclusion, our approach has shown to be very accurate and suitable for experiments where hundred of samples have to be screened making a manual evaluation unthinkable. This approach can be used not only in functional genomics studies but also in agricultural applications.

Keywords: drought, water limited, abiotic stress, phenotype, phenomics, near infrared, high throughput, images

# INTRODUCTION

Climate changes and environmental pollution will have a significant impact in the food production worldwide. The increase of global temperature will cause more frequent drought events (Prasch and Sonnewald, 2013) defined as soil water deficit or low water availability (Harb et al., 2010). Moreover, a persistent expansion of aridity has been observed since the middle of the 20th century and this process will continue according to current projection models (Dai, 2011). In some areas where crop yield reduction is predicted, major improvements in plant breeding programs and agricultural technology have to be developed (Jones and Thornton, 2003).

The use of different phenomics technologies in plants is a key element to improve our knowledge of the genotype–phenotype association of desired agricultural traits (Neilson et al., 2015) such as the response to water deficit. Some of these methods have been taken from medical applications as is in the case of high-resolutions X-ray computed tomography. This has shown to be an excellent tool to analyze the development of root system at high resolution in a non-destructive way (Lontoc-Roy et al., 2005). Other methodologies have been born specifically for plant phenotyping such as the Scanalyzer HTS developed by LemanTec (LemnaTec GmbH, Wuerselen, Germany) to scan small plants using a variety of wavelength cameras. Worldwide, prestigious universities and research institutes have acquired these technologies such as the CT scanning laboratory for agricultural and environmental research and the McGill Plant Phenomics Platform (MP3) at McGill University (Canada), the Australian Plant Phenomics Facility at the University of Adelaide (Australia) and the Arkansas Center for Plant-Powered Production at Arkansas University (USA).

In addition to the availability of the homozygous genome-wide knockout lines (Furbank and Tester, 2011), *Arabidopsis* is part of the mustard family (*Brassicaceae*; Haudry et al., 2013) to which economically relevant crops belong such as the edible canola oil, the cabbage vegetables or the biofuel candidate *Camelina sativa* for which the overexpression of the *Arabidopsis* MYB96 has conferred a drought resistance phenotype (Lee et al., 2014). This coupled with its small size and short life cycle (Lack and Evans, 2001) makes of *Arabidopsis* an excellent model organism to explore with the Scanalyzer HTS (Lemnatec GmbH, Wuerselen, Germany).

Despite the progress achieved in methods to detect genotype– phenotype association and quantify plant phenotypes at high resolution (Klukas et al., 2014), there are still limitations. The development of plant growth protocols, new image analysis algorithms and statistical pipelines is essential to exploit the full potential of these phenomics platforms. In this paper, we present a comprehensive approach to assess drought stress survival experiments using the advantages of our phenomic platform. This method is not limited to a specific protocol and it can be easily implemented into other high-throughput phenotyping facilities.

# MATERIALS AND METHODS

# Plant Lines

Two *Arabidopsis thaliana* drought responsive genes were used in this assay: *gtl1–5* (SALK\_044308C: AT1G33240), and *drs1* (SALK\_149366C: AT1G80710) in addition to ecotype Col-0 as wild type (WT). In previous studies, *gtl1–5* mutants have shown a resistant phenotype (Yoo et al., 2010) and *drs1* mutants a sensitive phenotype (Lee et al., 2010).

# Growth Conditions

Two different protocols named in this work "pot protocol" and "pellet protocol" have been used.

### "Pot Protocol"

Plants were grown in 2 1/2 square inch pots containing 20 g of 2:2:1 mixture of vermiculite, perlite, and peat moss (Cheng et al., 2013). After stratifying in the dark for 3 days at 4◦C, pots were transferred to a growth chamber (Conviron TC30) at 22◦C, 16/8-h light/dark photoperiod, 70% relative humidity and a light intensity of 135 mmol m−<sup>2</sup> s<sup>−</sup>1. A number of 10 pots with WT, 10 pots with mutant *gtl1–5*, 10 pots with mutant *drs1* and 2 pots with soil only to control evaporation were put in each tray with a basket. The trays were sitting in an inverted basket for aeration purpose and black mulch (Fabricville, Canada; felt fabric Z048- BLK) was added on top of every pot. This tray configuration was used for the water-limited group and for the control or wellwatered group. The distribution of the pots in each tray was regularly randomized using a computer algorithm as suggested by Skirycz et al. (2011) to minimize growth chamber effects. Here, this process has been done every 3–4 days and has taken 8 min per tray. The fill capacity of the pots was calculated by weighing the individual components, empty pots, mulch, pots with dry soil and 90 min after watering. Pots were weighed regularly during the experiment. After the appearance of the fourth leaf in over 60% of the plants, water was withheld for the drought or water-limited group until the plants exhibited lethal effects of dehydration or clear symptoms of wilting (Skirycz et al., 2011) after which the plants were rewatered (**Figure 1**). The pots (in the pot and in the mulch) and the trays (at side) had been labeled with barcodes to identify and trace the seedlings at each step of the experiment including the data analysis and the management of the randomization process.

## "Pellet Protocol"

The plants were grown in Jiffy-7 pellets (Skirycz et al., 2011) inside basket-pots with diameter of 5 cm inserted in trays with spots separated by 1 mm. The growing conditions are the same as the "pot protocol" but the growth chamber (Conviron A1000). A number of eight pellets for each line including the WT and three "soil only" were put in each one of the two trays (one for the well-watered group and other for the water-limited group). Mulch has not been used and the fill capacity was calculated similarly as the previous protocol. A blue label for identification purpose was added and the samples were also randomized using a computer algorithm. The beginning and the end of the waterlimited period have been done as described in the "pot protocol" (**Figure 1**).

# Data Acquisition

A customized version of the LemnaTec Scanalyzer High-Throughput Screening (HTS; LemnaTec GmbH, Wuerselen, Germany) installed at the McGill Plant Phenomics Platform (MP31 ) was used to carry out the image acquisition. The unit has a robotic arm that holds sensors or cameras and moves to different positions inside the measurement cabinet. In this study, images were taken with the visible light camera piA2400-17gc (VIS), 2454 × 2056 pixels, the fluorescent light camera scA1400- 17gc (FLUO), 1390 × 1038 pixels and the near infrared camera NIR-300PGE (NIR), 320 × 254 pixels. Every plant was imaged regularly during the experiments. The intensity of the NIR is measured from 0 to 255.

# Extraction of the Digital Plant and its Features

After image acquisition, all the copies were transferred to the MP3 server (a Dell R910 server with 512 GB of RAM and two MD1200 storage devices 72 TB). The FLUO images were

<sup>1</sup>http://mp3*.*biol*.*mcgill*.*ca

converted to HSB color space. Pixels having brightness (B) of the HSB color space higher than 25 + c were retained. The constant c is related to the HTS setting. The selected pixels were tagged as foreground to identify the objects using an adapted version of the "combined contour tracing and region labeling" algorithm proposed by Burger and Burge (2008). The objects having an area greater than 1000 pixels and an Euclidean distance lower than 300 to the theoretical pot center were selected. As a result, each plant was represented by one object or "digital plant" from which the area was calculated (Burger and Burge, 2008; Schneider et al., 2012; Camargo et al., 2014). After the FLUO images were resized to 2900 pixels width and proportional height, the pixels in the VIS images were shifted −290 + a pixels in the x-axis and −122 + b pixels in the y-axis coordinates. Then, the positions of the pixels composing the "digital plant" in the FLUO images were used to select the pixels that constitute the "digital plant" in the VIS images. The constant "a" and "b" are parameters related to the HTS deck configuration. The NIR images were processed as the VIS images using 380 to resize and −30 + a, −12 + b to shift the images. The third quartile or 75th percentile of the pixels of every "NIR digital plant" was then calculated (**Figure 2**).

# Color Classification and Clustering

The hue channel of the HSB color space from VIS was equally divided into 65 categories where each pixel of the "digital plant" was classified (color classification; Berger et al., 2012). Each class or category was defined as an interval of intensities and the union of all categories is equal to the full range of the hue channel (0–255). The median hue intensity of each interval with saturation (S) and brightness (B), equal to 255 or its RGB equivalent, was used to identify the class (**Figures 2E,F**). An Euclidean distance matrix of the "digital plants" was then calculated with the the color classes (percentage of pixels) of the VIS images using as input for a hierarchical cluster analysis. The method "ward" of the R function "hclust" was used for clustering. The resulting cluster tree of samples was then divided into two groups using the R function "cutree" (R Core Team, 2013).

# Statistical Analysis

The Cochran-Mantel-Haeszel ("mantelhaen.test" function in R) test was used to assess the effect of the water-limited condition on the survival rate of the lines and the Pearson's Chi-squared test as a goodness-of-fit test ("chisq.test" function in R) to support the matching between the manual inspection, near infrared and clustering results. The area was assessed using the Student's *t*-test ("ttest" function in LibreOffice, mode = 2 two tailed test and type = 3 heteroscedastic). The index on the third quartile (Q3) NIR intensity was calculated as the Q3 value divided Q3 base multiplied by 100.

# Implementation

The custom image analysis algorithm was developed using java 1.8.0-452 and ImageJ library v1.49 (Rasband, 1997–2012; Schneider et al., 2012). The statistical analysis script was written in R v3.0.2 (R Core Team, 2013) and LibreOffice 4.2.8.2 was used as a complement. PostgreSQL v 9.3.1.was used to construct the database. The analyses were implemented on the MP3 server.

# RESULTS

# Comparison of Protocols

In this study, we have implemented two protocols, namely the "pot protocol" and "pellet protocol," to assess the nondependency of our approach to a particular method. The format of the trays, the containers and the soil used in the protocols are completely different (see Materials and Methods). In the first case, the absolute fill capacity was 72.62 (±0.34) g water/pot compared to 30.15 (±0.74) g water/pellet. This is 58% less water available per plant. Under the same environmental conditions, the "pot protocol" took about 30 days to reach 1–2 g water/pot while the "pellet protocol" took about 15 days from 100% fill capacity to 2–3 g water/pallet where visible signs of wilting had appeared triggering the rewatering process (**Figure 1**). In the "pot

<sup>2</sup>www*.*java*.*com

(A,C,E,G) well-watered and (B,D,F,H) water-limited (Mutant line *drs1*, "pellet protocol"). Row images from VIS camera (A,B). "Digital plants" from FLUO after segmentation process (C,D). Color representation of the plants after color classification (E,F). NIR "digital plants" with the indexes based on the third quartile (Q3) of the intensity histogram (G,H).

protocol," 10 plants failed to germinate against one in the "pellet protocol." These were removed from the analysis, but the pots and pellet were kept in the trays during the whole experiment.

# Precision Phenotyping

The accuracy of the phenotypes depends largely on the effectiveness of the extraction of the plant digital representation. The Scanalyzer HTS allows overlapping of images from its different sensors. The separation between the background and the plant may become problematic using the VIS, especially at the end of the drought treatment when the colors of the leaves look similar to the soil. The use of the near infrared camera (NIR) has similar difficulties. The shifting of the intensity in the near infrared spectrum due to the change in the water content is proportional to the plants and to the soil making it difficult to separate the foreground from the background without losing information. Fluorescent light illuminated on the leaves has a high level of reflectance compared to the soil where it is null. This feature makes the extraction of the "digital plant" very accurate using the fluorescent light camera (FLUO) as shown in **Figure 2**. However, the occasional algal growth on the surface of the containers with particular soil mixtures may produce a high level of noise. One of the possible solutions is to use mulch as in the case with the "pot protocol." This was not necessary in the "pellet protocol" since algal growth was not observed. In these two cases, the separation between background and foreground has been easier and more accurate. The FLUO information has then been used to identify the "digital plant" from the images of the other sensors. Hence, the use of the FLUO as a segmentation mask has shown to be the optimal method compared to using a particular method specific to each camera.

# Projected Leaf Area

The projected leaf area showed a clear pattern between both conditions and in both protocols (**Figure 3**). A *<sup>p</sup>*-value lower than 0.01 was obtained from "days after sowing" (DAS) 27 in the "pellet protocol" and from DAS 28 in the "pot protocol" using a Student's *t*-test for each measurement. At the end of the waterwithholding period, we observed a reduction of the area in the water-limited groups. In the "pellets protocol," this was observed at DAS 27 where the water content was at 11%: 3.66 g water/pellet (**Figure 3B**). The line *drs1* showed this effect 1 day prior. In the "pot protocol," the reduction in the area was observed at DAS 46 where the water content was at 2%: 1.31 g water/pot (**Figure 3A**). However, measurements were not taken at DAS 44 and 45. The line *drs1* had also shown a decrease of the area before the others lines at DAS 43. During the recovery period, the area of the waterlimited group increased again in both protocols. On average, the water-limited group area was 8% of the well-watered group at the beginning of the recovery phase and reached 20% after 9 days in the "pot protocol," and 30% at the beginning and 46% after 2 days of recovery in the "pellet protocol."

# Manual Inspection of Plant Health

In order to evaluate the effectiveness of our method, the samples were classified visually into two categories based on signs of wilting or dryness (Supplementary Tables 1 and 2). In addition, other sources of stress were looked at to verify the correct classification of the data. At the last day of the experiments, "alive" or "dead" status was assigned to each sample as its last classification.

In the control (well-watered) group of both protocols, none of the plants were classified as wilted or dried. In the "pot protocol," at DAS 46 (last day of water-withholding), we have observed two of seven WT, two of eight *gtl1–5,* and 6 of 10 *drs1* plants presenting signs of wilting or dryness. Only one WT plant was recovered at the end of the experiment. In the case of the "pellet protocol," eight of eight WT, six of eight *gtl1–5,* and six of eight *drs1* plants were wilted or dried at DAS 28 (last day of waterwithholding). In this case, the recovery rate was higher, seven WT, six *glt1–5,* and two *drs1* plants. In the "pellet protocol," some of the samples had shown signs of water stress before reaching the end of the drought period, one *drs1* at DAS 25, three *drs1* and one WT at DAS 26 and five *drs1*, two *gtl1–5* and four WT at DAS 27. A Cochran-Mantel-Haeszel test on the waterlimited groups on the last day of the recovery periods yielded a p-value lower than 0.01. **Figure 4** shows survival percentages per line.

# Clustering of Color Classes

Here, the samples were clustered using as input the 65 classes of the color classification. Each class is represented by the percentage of pixels of the hue channel of the HSB color space (see Materials and Methods). In this way, the cluster should not be affected by the size of the plant. The resulting tree was then divided into two groups (Supplementary Tables 3 and 4). The expectation is that the samples are classified according to the colors of the leaves. For example, a healthy plant should show a green pattern against the yellow/brown colors exhibited by a dried plant (**Figures 2E,F**).

both treatment groups (DR: water-limited or drought and WW: well-watered) for (A) "pot protocol" and (B) "pellet protocol".

The efficacy of the clustering was assessed comparing the group assigned by the clustering and the health status from the manual inspection of each plant at every DAS (**Tables 1** and **2**). As expected, before clear signs of dryness, there was no strong agreement in any of the two groups since most of the plants were green without any stress. The minimum percentage of agreement was 40%. At the recovery period, we observed 81, 91, and 91% at DAS 28, 29, and 30 in the "pellet protocol" and 98, 100, and 76% at DAS 46, 53, and 56 in the "pot protocol." Looking closer at the mismatch at DAS 56, we observed that this is produced by 10 plants in the wellwatered group where 80% coincide with the samples to which the inflorescences have been cut. This event has created a visible stress on the plants. The goodness-of-fit test (Pearson's Chisquared test) for the "pot protocol" has a *p*-value lower than 0.05 at DAS 25, greater than 0.05 (non-significant) from DAS 28 to 35, a *p*-value lower than 0.05 from DAS 37 to 43, a p-value lower than 1e-10 at DAS 46 and 53, and a *p*-value lower than 1e-03 at DAS 56 (**Table 1**). In the case of the "pellet protocol,"

the same test yielded a *p*-value greater than 0.05 (non-significant) at DAS 16 and 23, a *p*-value lower than 0.01 at DAS 13, 20, 21, 22, a *p*-value lower than 1e-04 from DAS 25 to DAS 28, and the most significant *p*-value (*<*1e-07) at DAS 29 and 30 (**Table 2**).

# Near Infrared Classification

The near infrared has shown a reflectance increase in the water deficit group not exhibited in the control group (**Figure 5**). Here, we have computed the third quartile (Q3) or 75th percentile from


TABLE 1 | Accuracy of the cluster of color classes and the near infrared classifiers compared to the manual inspection during the "pot protocol" experiment.

*p: Pearson's Chi-squared test p-value.* <sup>a</sup>*An extra stress in the control group have reduced the percentage of agreement.*

TABLE 2 | Accuracy of the cluster of color classes and the near infrared classifiers compared to the manual inspection during the "pellet protocol" experiment.


*p: Pearson's Chi-squared test p-value.*

the pixels of each "NIR digital plant" (Supplementary Tables 5 and 6). Indexes on these values were calculated using DAS 25 and DAS 13 accordingly as references. This allows us to homogenize the start point and to monitor the intensity changes. A condition of the "index base reference" was the homogeneity of their subsequent measurements. In the "pellet protocol," there were no changes at DAS 16, and only 2% at DAS 20. In the "pot protocol," we have observed a difference of only 1% at DAS 28 and 2% at DAS 30.

In the last day of withholding water, an increase of about 20% or more in near infrared intensity was registered in both protocols for the samples exposed to water deficit condition (**Figure 5**). The composition of these two groups shows that 100% of the plants having an index higher than 120 in the water-limited group (**Tables 3** and **4**) have presented signs of wilting or dryness in the manual inspection (DAS 53 and 56 for the "pot protocol" and DAS 29 and 30 for "pellet protocol"). The classification of the samples in stress (wilted or dried: index *>* 120) and nonstress (index *<* = 120) shows 100% agreement with the manual inspection in all DAS for the three lines in the well-watered group of the "pellet protocol" (**Table 4**). In the case of the "pot protocol," 100% agreement was observed in all DAS for the lines *gtl1–5* and *drs1*. The WT in this group showed 78% at DAS 37, 89% at DAS 30, 35, 40, 42, 43, 46, 53, and 56, and 100% at DAS 25, 28, 32, and 34 (**Table 3**). As indicated in **Table 3**, 78 and 89% is produced by only three samples at different measurement points. However, they never exhibited an index higher than 127. In the deficit water period, the "pellet protocol" showed 100% agreement from DAS 13 to 23. From DAS 24 to DAS 28 (last water deficit day), the percentage of agreement varies according to the line and DAS with a maximum of 100% and a minimum of 50% (**Table 4**). The plants showing mismatch were classified by the NIR classification as being stressed. While three plants never showed a visible stress phenotype, the other seven presented visible signs of stress at the last day of the water deficit (DAS 28). The "pot protocol" showed 100% agreement in the deficit period except for WT at DAS 37 and 46 and *glt1– <sup>5</sup>* at DAS 46 (86, 86, and 88%, respectively; **Table 3**). As in the previous case, the two plants not showing visible signs of wilting or dryness have an index higher than 120. We never observed that the NIR classifier pointed to a plant as nonstressed when in reality it was. This may suggest that in some cases the NIR could show a pre-wilting stage of the plant before observed under a visual inspection. In the recovery period of the water-limited group, all lines in all protocols have exhibited 100% agreement between the manual inspection and the NIR classifier even 24 h after rewatering (**Tables 3** and **4**). In this case, a goodness-of-fit test (Pearson's Chi-squared test) has shown a *p*-value lower than 1e-09 for the "pot protocol" and a *p*-value lower than 1e-06 for the "pellet protocol" in all DAS having a *p*-value lower than 1e-10 in the recovery period of both protocols (**Tables 1** and **2**).

# DISCUSSION

The aim of this study was to present a comprehensive approach to assess drought stress survival experiments using the benefits of a phenomics platform. One of the advantages of using the HTS at McGill University is the integration of multiple sensor information. Here, this integration was vital to separate the plants from the background using the FLUO to produce a segmentation mask for the images from the other sensors (**Figure 2**; Berger et al., 2010; Klukas et al., 2014). As such, we have provided

more power to our analysis than the use of only one sensor. Another advantage is the configuration of the sensors is kept as metadata in the database and, as such, makes each measurement comparable and reproducible.

The two protocols used in this study have been proven to be effective since both have been appropriate to assess the survivability of the mutant lines and the outputs have been consistent. However, the "pellet protocol" was more efficient because it took roughly half of the time to reach the same results (**Figure 1**). This is critical when a large number of lines has to be screened. In addition, the mulch added in the "pot protocol" keeps the moisture in the soil thereby reducing the rate of water loss (Junker et al., 2015), but was necessary to avoid noise in the segmentation process produced by algae growing on the soil surface. Nevertheless, there is no evidence that this part of the protocol has physiological effects (Junker et al., 2015) or affects the experimental output in our hands. Nevertheless, the use of materials that do not have a fluorescence signature such as the Jeffry pellets is recommended over the use of mulch on top of the pots.

The differences in the projected leaf area between both treatment-groups were detected as the water-withholding period progressed (**Figure 3**). These observations are in concordance with previous studies where significant differences in total leaf


TABLE 3 | Agreement between the near infrared classifier and the manual inspection by line and treatment during the "pot protocol" experiment.

TABLE 4 | Agreement between the near infrared classifier and the manual inspection by line and treatment during the "pellet protocol" experiment.


area of WT (Col-0) were observed between well-watered and water deficit groups (Bouchabke et al., 2008). These differences were significant before the end of the water-withholding period (16 days for the "pot protocol" and 1 day for the "pellet protocol"). This may be explained by the capacity of the soil mixtures to withhold water or by the rate of the drought process. However, we did not observe in the two experiments similar growth patterns for a particular line in any of the groups except for the line *drs1* where the reduction of the area has started earlier compared to the WT and line *gtl1–5*. This may be an indicator of sensitivity to the water-limited condition. The decrease of the projected leaf area is likely a clear sign of a pre-wilting stage and a sensitive line such as the *drs1* may wilt earlier than other lines.

A color trait clustering method has been shown to be a powerful tool as a classifier (Dana and Ivo, 2008). A neural network classifier was also used to successfully identify heatdamaged and green-frost-damaged soybeans (Shatadal and Tan, 2003). The cluster of color classes allowed us to differentiate between the dead and alive plants with an accuracy greater than 90% at the recovery phase using only the percentage of colors (**Tables 1** and **2**). This means that the size of the rosette has not been taken in consideration which is important to avoid a potential bias in the classification. The cluster was always forced to split the samples into two groups. When there was not a visible effect of drought, the division might have revealed an extra source of stress as was the case in the well-watered group of the "pot protocol" (**Table 1**) or might have showed "pre-existent" differences in the plants. However, when the changes of colors were produced as consequence of water deficiency, the cluster classified both groups accurately.

The near infrared light absorption is increased by the presence of water in the leaves (Fahlgren et al., 2015). Previous studies reported a correlation between the near infrared-based indexes and relative water content with severely drought-stressed plants (Berger et al., 2010) as is the case in many survival studies. We have used this property to classify the plants using index numbers based on the third quartile of the NIR pixel distribution. This quartile is more sensitive to an intensity increase than the other two since it is located in the upper part of the "scale." The accuracy of the NIR classification during the recovery phase for both protocols has been at least 98% (**Tables 1** and **2**). Furthermore, the NIR detected in some cases a pre-wilting stage prior to exhibiting visible signs, it was not affected by the "extra stress" subjected to the inflorescences in the "pot protocol" and it never showed a false negative. Evidently, the performance of the NIR classifier was superior from the beginning in both experiments. In addition, this method is not affected by the size of the rosettes (**Figures 3** and **4**) as is the same case in the cluster of color classes.

The two mutant lines, *dsr1* and *gtl1–5*, were included in our experiments to show the applicability of our method. Our results have shown that the line *drs1* has a significantly lower survival rate compared to WT (**Figure 4**). This is in concordance with the literature as this mutant has been identified as drought sensitive (Lee et al., 2010). *gtl1–5* is a drought resistant line (Yoo et al., 2010). In our case, it showed the same survival rate as WT (**Figure 4**). This may be explained by the high recovery rate of the WT. The rewatering point was based on the observation of lethal effects of dehydration. However, in some cases, the prediction of the recovery of a plant is not evident and should be a good subject for further investigations. The increase of

the number of samples with two experimental repetitions using different rewatering points is a possible solution. If the phenomics information is processed in real time, the number of wilted plants could be easily obtained using the NIR classifier to determine exactly these points.

Both classifiers as proxies of the "health status" of the plants have shown to be independent of the size of the rosettes. Passioura (1991) has claimed that plant growth is affected by soil structure in many ways such as the root growth rate or the ability to uptake water and nutrients. This may explain the differences in the projected leaf area ranking of our lines between these protocols. In a survival assay, where a classification of the samples between dead and alive is sought, the projected leaf area is not always a clear index of plant status. Skirycz et al. (2011) pointed that the survivability is not an indicator of growth performance and most of the survival phenotypes in drought are associated with constitutive activation of watersaving mechanisms. Hence, the use of the "NIR" and "the cluster of color classes" classifiers overcomes this limitation of the projected leaf area.

In conclusion, we have shown that our approach is very accurate and can be used with different soil mixtures and containers. The cluster of color classes and the NIR have been shown to be very good classifiers in survival drought experiments. However, the NIR was excellent and efficient during the entirety experiments including the early stages due to its association with water content. When hundreds of samples are tested and analyzed at several time points, the use of a phenomics platform coupled with a bioinformatics approach becomes

# REFERENCES


strictly necessary and this without taking in consideration the objectivity that a human cannot assure.

# AUTHOR CONTRIBUTIONS

EV, AT, and TB designed the project. EV and TEB wrote the manuscript. EV and AT performed the experiments and the phenotyping analyses. AD selected the mutant lines and assisted with the experimental calibration at beginning of the project. Funding was obtained by TB. All authors discussed the results and provided input on the manuscript.

# ACKNOWLEDGMENTS

We would like to thank Zoe Joly-Lopez, Maia Kaplan, and Yang Shao for reviewing the manuscript. We also thank Mark Romer and Claire Cooney for their help with the growth chambers. This project was funded by grants from Genome Quebec, Genome Canada, Canada Foundation for Innovation (CFI) and Natural Sciences and Engineering Research Council (NSERC) of Canada to TB.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal*.*frontiersin*.*org/article/10*.*3389/fpls*.*2015*.*01101


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Vello, Tomita, Diallo and Bureau. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Computed tomography scanning can monitor the effects of soil medium on root system development: an example of salt stress in corn

*Sowmyalakshmi Subramanian, Liwen Han, Pierre Dutilleul and Donald L. Smith\**

*Department of Plant Science, McGill University, Montréal, QC, Canada*

#### *Edited by:*

*Sergey Shabala, University of Tasmania, Australia*

#### *Reviewed by:*

*Jeffrey M. Warren, Oak Ridge National Laboratory, USA Megan Shelden, The University of Adelaide, Australia*

#### *\*Correspondence:*

*Donald L. Smith, Department of Plant Science, McGill University, Macdonald Campus, 21111 Lakeshore Road, Ste-Anne-de-Bellevue, Montréal, QC, Canada donald.smith@mcgill.ca*

#### *Specialty section:*

*This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science*

*Received: 12 December 2014 Accepted: 01 April 2015 Published: 28 April 2015*

#### *Citation:*

*Subramanian S, Han L, Dutilleul P and Smith DL (2015) Computed tomography scanning can monitor the effects of soil medium on root system development: an example of salt stress in corn. Front. Plant Sci. 6:256. doi: 10.3389/fpls.2015.00256* Seeds and young seedlings often encounter high soluble salt levels in the upmost soil layers, impeding vigorous growth by affecting root establishment. Computed tomography (CT) scanning used at low X-ray doses can help study root development in such conditions non-destructively, because plants are allowed to grow throughout the experiment. Using a high-resolution Toshiba XVision CT scanner, we studied corn (*Zea mays* L.) root growth under optimal and salt-stressed conditions in 3D and on a weekly basis over 3 weeks. Two groups of three corn plants were grown in the controlled environment of a growth chamber, in mid-sized plastic pots filled with sieved and autoclaved sand. Seedlings were subjected to first CT scanning 1 week after seed planting. Our main research objectives concerning root systems were: (i) to quantify structural complexity from fractal dimensions estimated on skeletal 3-D images built from CT scanning data; (ii) to measure growth from volumes and lengths and the derived relative rates and increments, after isolating primary and secondary roots from the soil medium in CT scanning data; and (iii) to assess differences in complexity and growth per week and over Weeks 1–3 for groups of corn plants. Differences between groups were present from Week 1; starting in Week 2 secondary roots were present and could be isolated, which refined the complexity and growth analyses of root systems. Besides expected Week main effects (*P <* 0.01 or 0.05), Week × Group interaction (*P <* 0.05 or 0.10), and Group main effects were observed. Graphical, quantitative, and statistical analyses of CT scanning data were thus completed at an unprecedented level, and provided new and important insights regarding root system development. Repeated CT scanning is the key to a better understanding of the establishment in the soil medium of crop plants such as corn and the assessment of salt stress effects on developing root systems, in complexity, volume, and length.

Keywords: corn (*Zea mays* L.), NaCl salt stress, developing root systems, structural complexity, fractal dimension (FD), root volumes and relative growth rates, root lengths and increments, computed tomography (CT) scanning

# Introduction

Plants encounter various abiotic and biotic stressors during their life cycle. Two of the most prevalent abiotic stressors confronting global agriculture are soil salinity and drought. Salinization is one of the more serious agricultural limitations, especially in the arid and the semi-arid regions of the world. Approximately 20% of the world's irrigated land yields a third of the global food stocks. At the same time, about 30% of this highly productive land is affected by salinity (Yan, 2008; Xu et al., 2011). Land clearing and irrigation are among the major contributors to salinization of agricultural lands. Their general impact has been reported (Munns, 2005; Rengasamy, 2008), and this is aggravated by a number of factors, including climate, the degree of water deficit, the inherent salt content of soils, topography, and the underlying geology and hydrology (Wiebe et al., 2007).

During the initial development of a seed into a plant, the germinating seed puts forth the radical that subsequently differentiates into the root system, the fundamentals of which determine plant growth and productivity. A plant's response to salinity is a complex process which affects the plant's tissue and organ development at various stages of growth. Seed germination under saline conditions cause significant reductions in seed germination percentage, thereby causing uneven stands and reduced yield (Foolad et al., 1999). Sodium chloride (NaCl) is a dominant salt in nature, which at sufficiently high concentrations reduces the ability of plants to take up water (water-deficit effect) and other essential nutrients (ion-excess effect; Munns, 2005; Munns et al., 2006). The uppermost soil layers are generally the site of highest soluble salt accumulation due to evaporation and capillary rise of water, so that seeds and young seedlings are frequently confronted with salinities much too high to allow vigorous growth (Almansouri et al., 2001). However, depending on the soil type and irrigation, or the subsoil salinity independently of the water capillary rise, a saline gradient is usually seen in many saline fields that affect root spread of crops (Wiebe et al., 2007; Rahnama et al., 2011).

Salt stress causes changes in plant growth through (1) osmotically induced water stress, (2) specific ion toxicity due to high concentrations of sodium and chloride, (3) nutrient ion imbalance due to high levels of Na+ and Cl−, which reduce the uptake of K+, NO<sup>−</sup>, and PO4 <sup>3</sup>−, and (4) increased production of reactive oxygen species (ROS), which damage macromolecules inside plant tissue, all of which result in plant growth reduction (Greenway and Munns, 1980). For instance, salt stress enhances the accumulation of NaCl in chloroplasts of higher plants, which affects their growth rate, and is often associated with a decrease in photosynthetic electron transport activities (Kirst, 1990). Additionally, in higher plants, it inhibits photosystem-II activity (Kao et al., 2003; Parida et al., 2003), which indirectly reflects upon the below-ground biomass, the roots. Simulation of a salinity gradient using a PVC tube with paper roll soaked in salt demonstrated sensitivity to salt for roots of wheat cultivars with regard to branching (Rahnama et al., 2011). Screening of genotypes of wild and domestic barley for salinity stress suggests large variations in response to salt in Petri plate assays. Increasing salt concentration (100–150 mM NaCl) decreased both shoot and root growth in various aspects in barley cultivars (Shelden et al., 2013), suggesting saline soils substantially alter plant metabolic processes (Levitt, 1980).

Corn is the third major cereal crop produced across the globe (FAO, 2014), and is grown under a very wide range of climatic conditions. The seedling structure of the family Poaceae is unique among monocotyledonous plants (Coudert et al., 2010). Specific terms such as scutellar node, coleoptilar node, mesocotyl, and coleorhizae, coined by Tillich in 1977 (Tillich, 1977; Hochholdinger et al., 2004), have been since used to describe these structures. The root system architecture of corn is complex, and was described by Abbe and Stein (1954). During the plant's life cycle, an embryonic root system and a post-embryonic root system can be distinguished. The embryonic root system consists of a primary root and a number of seminal roots, all of which dominate during the first 2 weeks of seedling development and establishment. It is followed by the early post-embryonic root system in which the primary and seminal roots develop lateral roots 6–7 days after initialization of the embryonic root system. The post-embryonic root system consists of shoot-borne roots formed at the nodes, called the crown roots, 10–14 days after seed germination. Roots developing on the above-ground nodes are called the brace roots. All the lateral roots that developed from the embryonic and post-embryonic systems influence the architecture of the whole system through the branching patterns, including the secondary and higher-order roots; they govern the anchorage as well as the nutrient and water uptake, and are sensitive to environmental factors (Hochholdinger et al., 2004; Lynch, 2013).

Differences in root architecture allow the crops to explore various soil domains at different intensities, in coordination with other environmental factors (Schwinning and Weiner, 1998; Postma and Lynch, 2012). The study of a root system's architecture is of importance to plant breeders because genetic variation and a suite of quantitative trait loci are involved in its development and functioning (de Dorlodot et al., 2007). The plasticity and dynamics of root system architecture are equally important for the manipulation of crop agronomic traits (Richards, 2008; Zhu et al., 2011), since a proper understanding is required to develop and breed crops for targeted environments (Smet et al., 2012). In root system architecture studies involving field-grown corn, the secondary and higher-order roots that developed after the primary and seminal roots contribute significantly to the total root number and total root length, although the root length of specific orders can vary according to soil types (Wu and Guo, 2014).

Corn is a C4, cross-pollinated, and highly polymorphic crop, with variable salinity tolerance across genotypes. When grown under saline conditions, it can show decreased growth and yield, several of its growth variables being affected at early seedling and growing stages, with the roots being affected the most (Carpici et al., 2009). Nuclear magnetic resonance studies suggest that corn root tips accumulate sodium rapidly (Spickett et al., 1993). Differences in the pattern of root solute potential were observed in corn treated with NaCl as a salt accumulation treatment vs. a salt shock treatment with the administration of 100 mM NaCl as a single dose; the sudden shock caused rapid inhibition of root extension, accompanied by decreased root solute and potential (Rodríguez et al., 1997).

Computed tomography (CT) scanning was originally designed as a medical diagnostic tool in the early 1970s (Kalender, 2000), and is now applied in a variety of non-medical fields, such as archeology, biology, the soil, material and earth sciences, the timber industry, industrial inspections, and aviation security, owing to its non-invasiveness combined with high spatial resolution based on indirect matter density measurement (van Kaick and Delorme, 2005). Very importantly, the low dose of X-rays used in studies involving a medical-type CT scanner, compared to industrial CT scanners which use much higher doses, has been proved to be adequate to study developing plant structures (Dutilleul et al., 2005, 2008; Lontoc-Roy et al., 2005; Han et al., 2008, 2009). A recent study on rice root variables and the associated microbial communities suggests that there were no significant differences between non-CT scanned and CT scanned samples (Zappala et al., 2013).

The use of CT scanning technology with plants was initiated in the late 1990s in an approach alternative to destructive characterization and rhizotron-based observation of root branching patterns. Thus, Heeraman et al. (1997) were able to visualize, non-invasively in 3D, 0.8 cm of bush bean roots (*Phaseolus vulgaris*), but much remained to be done in terms of graphical, quantitative, and statistical analyses of plant CT scanning data. The comparison of destructive vs. CT scanning-based characterization of a root system was among the research questions initially investigated. Actually, a reliable comparison between the two procedures, applied for the same root system, is difficult and even practically impossible, as the CT scanning must precede the destructive characterization; accordingly, the two will never be applied at exactly the same time and in exactly the same experimental conditions (e.g., root moisture), as destructive characterization usually involves root washing. Some studies (chickpea: Perret et al., 2007; cereals: Flavel et al., 2012) suggest an underestimation of total root length with CT scanning, which would mean the incompleteness of root isolation by CT scanning has been larger than the root loss by soil separation and washing in destructive characterization. Nevertheless, Lontoc-Roy et al. (2006) clearly showed that the 3-D image of a corn root system reconstructed from CT scanning data collected in an appropriate soil-moisture condition for the type of soil medium used provides a more reliable basis for fractal analysis and the estimation of a fractal dimension (FD) as measure of structural complexity, than 2-D photographs of the same root system extracted from the soil and analyzed using the fractal analysis module of the WinRhizo software (Regent Instruments Inc., Québec City, QC, Canada). Definitely, CT scanning technology helps capture details of root systems, such as lateral root growth and its orientation, the variables of which cannot be studied using a destructive method (Perret et al., 2007; Han et al., 2008).

Therefore, we have used a high-resolution X-ray CT scanner to study the architecture of developing root systems of corn variety 19K19 under optimal and salt-stressed conditions, at an unprecedented level of graphical, quantitative, and statistical analyses. Our main research objectives concerning root systems were: (i) to quantify structural complexity from FDs estimated on skeletal 3-D images built from CT scanning data; (ii) to measure growth from volumes and lengths and the derived relative rates and increments, after isolating primary and secondary roots from the soil medium in CT scanning data; and (iii) to assess differences in complexity and growth per week and over Weeks 1–3 for groups of corn plants. Hence, our main biological research objective is to further our knowledge and understanding of the below-ground response of corn to salinity stress at the initial stages of plant development.

# Materials and Methods

#### Plant Material, Soil Preparation, and Salt Stress Imposition

Seeds from one of the highest yielding organic corn varieties, 19K19 Blue River, procured from Doug Shirray Seeds and Ag supplies (Tavistock, ON, Canada), were used in our experiments since these are among the most easily available untreated corn seeds. The growth medium was sand, which was sieved to 2 mm to homogenize the rooting media, autoclaved, and kept dry for at least 1 week before potting. In black plastic pots with an 18-cm diameter at the top, a 0.1-mm side-wall thickness, and 17-cm total height (Classic-R 400; Plant Products Co. Ltd., Laval, QC, Canada), filled with sieved and autoclaved sand, one corn seed was placed in the center at a depth of 2.5 cm. Three such pots were prepared for control and three more for the salt stress. Thereafter, the pots were uniformly watered, and the seedlings were allowed to emerge. One day after emergence, the pots were given <sup>1</sup> */*<sup>2</sup> Hoagland solution (Hoagland and Arnon, 1950) for control and <sup>1</sup>*/*<sup>2</sup> Hoagland + 100 mM NaCl for salt stress as a one-step salt shock. According to our results from other experiments with this cultivar of corn, 100 mM NaCl imposition was a salt stress for which this cultivar expressed a response.

The plants were grown in a growth chamber (Conviron Model No. PGR15, Controlled Environments Ltd, Winnipeg, MB, Canada), set at 25 ± 2◦C (day temperature) and 22 ± 2◦C (night temperature), with a photoperiod of 14/10 h day/night cycle, 60–70% relative humidity and photosynthetic irradiance of 350–400 µE m−<sup>2</sup> s<sup>−</sup>1. Subsequently, the pots were given <sup>1</sup> */*2 Hoagland once a week. The watering was scheduled so as to reduce moisture content of the sand rooting medium at the times of CT scanning, thus following Lontoc-Roy et al. (2006); this way of proceeding has the positive effect of increasing the contrast between roots and soil (in the case of sand) in the CT images, and consequently improves the CT scanning data analyses.

### CT Scanning Configuration, Data Collection and Processing, and Skeletal Root Image Construction

The CT Scanning Laboratory for agricultural and environmental research on Macdonald Campus of McGill University (Ste-Annede-Bellevue, QC, Canada) contains a Toshiba XVision highresolution CT scanner (Toshiba Corporation, Medical Systems Division, Tokyo, Japan). Our experimental corn seedlings were CT scanned in this facility, once a week for 3 weeks; more specifically, it is the plant–soil materials in the pots that were CT scanned. The first sessions of CT scanning (May 7–8, 2012) corresponded to 7 or 8 days after the seeds were planted in the pots. It was not possible to CT scan all six plants in 1 day, so the order for CT scanning was randomized the first week and repeated the next 2 weeks; on each day, one plant from one group and two plants from the other group were CT scanned; between CT scanning sessions and until completion of the experiment, the potted corn seedlings were maintained in the growth chamber. Each potted plant was positioned horizontally on the couch of the CT scanner, and entered the gantry 'stem first,' for a top-to-bottom CT scanning of its content. Earlier tests and previous experiments suggest this approach to be better for tracing the roots embedded in the soil during the procedure of CT scanning data analysis (Lontoc-Roy et al., 2005; Han et al., 2008). The helical CT scanning mode was chosen with an image reconstruction interval length of 0.4 mm along the *Z* axis. Other configuration parameter values were based on the experimental conditions, such as the density of the rooting medium used and the size of the object to be scanned. Hence, the field of view was set at SS (18-cm diameter), with no zoom factor (value of 1.0); the tube voltage, at 120 kV; and the tube current, at 150 mA. Every CT scanning session for each root system comprized of 300 cross-sectional CT images covering a depth of 12 cm.

The raw CT scanning data were obtained using the FC70 function, and consisted in CT numbers (CTN), expressed in Hounsfield units (HU). By definition, a CTN is an average relative measure of the density for a pixel in a CT image or for the corresponding volume (voxel) with equal lengths and a width as third dimension (0.35 mm × 0.35 mm × 0.4 mm in this study):

$$\text{CTN (HU)} = \frac{\mu\_{object} - \mu\_{water}}{\mu\_{object} - \mu\_{air}} \text{ times } 1000 \tag{1}$$

where *µobject* = mean value of the linear attenuation coefficient for the voxel; *µwater* = linear attenuation coefficient of pure water; and *µair* = linear attenuation coefficient of a standardized air sample. The CTN scale is linear and is centered at 0 to correspond to water, and the CTN for air is calibrated to −1000 HU. Hence, densities greater or less than water correspond to positive and negative CTNs, respectively.

The raw CTN data thus collected (300 512 × 512 matrices of CTNs) were first displayed on the CT console and further processed on a 4-core Dell i3 computer. Using MATLAB 2013b and 2014a (The MathWorks Inc., Natick, MA, USA), the data for a given corn seedling CT scanned in a given week were converted into an internal 3-D array, and scrutinized with a Graphical Unit Interface application program created to view the cross-sectional CT images and produce digitally on computer the 3-D architecture of the first-order roots and as many secondary roots of the system as possible. The construction of skeletal 3-D images of the corn root systems was performed in this environment.

#### Fractal Dimension Estimation

As already noted by Foroutan-pour et al. (1999a) for fractal analysis from photographs of branching patterns of soybean seedlings, it is very important to skeletonize the image of the structure of interest, whether it is a shoot branching pattern or a root system, in order to perform an unbiased estimation of the FD; by "skeletonization," we mean reducing the thickness of a branch or a root to 1 pixel in the 2-D image, or 1 voxel in 3D. Based on the methodological results of Foroutan-pour et al. (1999b), it is also very important not to use the full range of box sizes in the boxcounting procedure of FD estimation in 2D, or the full range of cube sidelengths in the cube-counting procedure in 3D; instead, it is recommended to use, in 3D, a subset of cube sidelengths that does not contain the smallest and largest lengths. This follows the approach of Han et al. (2008), who analyzed 3-D skeletal images of potato root systems; we have followed and applied their FD estimation procedure.

Among the nine cube sidelengths available in our study, after trying all possible subsets including 3, 4, or 5 middle length values, we retained the FD estimates obtained after discarding the three smallest and three largest cube sidelengths, because this was found to provide the highest *R*2-values in the log–log plot

$$
\log[N(s)] = k + \text{FD}\log(1/s) \tag{2}
$$

where log(*.* ) is the natural logarithm, *N*(*s*) denotes the number of cubes with sidelength *s* intersecting the skeletal root system, and the straight line is fitted by ordinary least squares.

#### Root Volumes and Relative Growth Rates

In preparation of our experiment, we had made preliminary tests using root systems of corn seedlings other than the experimental ones but of same variety, grown in same rooting medium with same Hoagland nutrition solution as the control group. One of our goals in doing this was to photograph corn root systems at same developmental stages (1, 2, and 3 weeks after seed planting) immerged in a container filled with water (**Figure 1**), for later comparison with 3-D images constructed from CT scanning data. Another goal was to make manual measurement of thickness on those corn roots once washed, so the procedure described below was applied with confidence to construct non-skeletal 3-D images of the experimental corn root systems from CT scanning data. Our procedure also takes into account the difference in size between primary and secondary roots and the spatial resolution of the CT scanning data collected (voxel dimensions: 0.35 mm × 0.35 mm × 0.4 mm).

The skeletal primary roots were 'grown' by up to four layers (one layer a time), in 3D, if the CT numbers of the neighboring voxels did not exceed 850 HU; for secondary roots, the threshold used for CT numbers was 950 HU and the maximum number of layers added to the skeleton was two. From the primary and secondary roots thus 'grown,' three types of root volumes were calculated for each of the six experimental corn seedlings in each of the 3 weeks: for lower roots alone, for upper roots alone, and for lower and upper roots combined. The stem and seed were not included in the evaluation of root volumes. Relative growth rates between successive Weeks *t* and *t*+1 were calculated accordingly from the estimated root volumes (denoted Root vol\_*t* and Root vol\_*t*+1 below):

Relative growth rate\_Week *t* + 1*,*

Week *t* = *(*Root vol\_*t* + 1 − Root vol\_*t)/*Root vol\_*t* (3)

#### Root Lengths and Length Increments

photograph is shown on a black background. The three corn seedlings

MATLAB, in combination with the image analysis toolset ImageJ (National Institutes of Health, Bethesda, MD, USA), were used for root length measures, in the following sequence of steps and operations. First, the 3-D array for the skeleton of root system of a given corn seedling CT scanned in a given week was loaded in MATLAB. Using the MATLAB function *imwrite*, each slice, out of 300 available per corn seedling per week, was then exported as an 8-bit gray image with a given format (i.e., BMP), into a designated folder. An image stack was then built from all the images after these were imported into ImageJ. The corresponding 3-D image was skeletonized with the ImageJ skeletonization procedure, and a customized 3-D analysis plugin (https://github.com/fiji/AnalyzeSkeleton/) was used to perform the root length measurements. Finally, the output was summarized to obtain total lengths of lower and upper roots separately. Say Root length\_*t* and Root length\_*t*+1 are measures of total root length for a given corn seedling in Weeks *t* and *t*+1; the corresponding increment was then calculated as

Root length increment\_Week *t* + 1*,*

Week *t* = Root length\_*t* + 1 − Root length\_*t* (4)

#### Statistical Analyses

Sample means per group (Control, optimal vs. Salt, salt stress) and the corresponding SEs were computed and plotted on a weekly basis. Furthermore, analysis of variances (ANOVAs) for temporal repeated measures (ANOVARs; Crowder and Hand, 1990; Dutilleul, 2011) were performed on the data tables of FD estimates, root volumes and lengths (1 weekly observation = 1 temporal measure), as well as on those of fractal dimension ratios (FDR) between the initial FD estimate and the next two, relative growth rates (Eq. 3), and root length increments (Eq. 4). Univariate ANOVARs using the sample variance–covariance matrix in the modified *F*-ratio tests were carried out in SAS 9.3 PROC GLM, option REPEATED (SAS Institute Inc., Cary, NC, USA), because the small sample sizes did not allow a modeling of the variance–covariance structure in a mixed-model approach. For the same reason, we considered three significance levels, namely 0.01, 0.05, and 0.10, in our hypothesis testing. Differences among treatments were only considered meaningful when occurring at one of these significance levels, and when these

are discussed, the *P*-value is bolded. The between-group homogeneity of variance was tested, and rejected only once (i.e., for volume of lower roots in Week 1) in 15 tests, with no consequence for our results.

those of Weeks 1–3 (from A–C) for the control corn seedlings in our experiment.

#### Results

Differences in structural complexity and space occupancy of the developing root systems of corn seedlings, grown in sand under optimal condition vs. salt stress, were investigated in 3D. Structural complexity was measured through FDs of skeletal 3-D images, and space occupancy, through root volumes. Betweengroup differences were assessed absolutely and relatively in and over the 3 weeks of the experiment. Results are presented below and in tables and figures.

#### Three-Dimensional Image Analyses

From the skeletal 3-D images of root systems constructed from the CT scanning data collected on a weekly basis for individual plants in the control and salt-stressed groups (**Figures 2A,B**), it is clear that germinated seeds of the control group in Week 1 (left panels of **Figure 2A**) show the development of embryonic roots with two subsets of roots called "upper roots" and "lower roots" here, whereas the onset of upper roots is delayed in the salt-stressed plants and upper roots in this group are present and visible only from Week 2 on (middle and right panels of **Figure 2B**). Also, the branching patterns of root systems in Weeks 2 and 3 show less prominent lateral branching in saltstressed plants, compared to control plants. **Figures 3A,B** show the root volumes that it has been possible to isolate from the plant–soil CT scanning data, by expansion from the initial skeletons of root systems. Graphically, all of this strongly suggests that during the initial stages of plant growth, the salt abiotic stressor (100 mM NaCl) has negative effects on the development of the root system, and weakens the below-ground establishment of the corresponding seedling relative to one grown under optimal conditions.

#### Fractal Analyses and ANOVARs on FD Estimates

Fractal analysis was restricted to lower roots, in the absence of upper roots for all three salt-stressed plants in Week 1 and for

group. In each panel, the horizontal plane in the 3-D rendering represents the X-Y plane in CT scanning terminology; this plane is perpendicular to the couch of the CT scanner. The small filled spheres locate the connection points of the primary upper and lower roots to the below-ground part of the stem.

FIGURE 3 | Non-skeletal 3-D images of the corn root systems depicted in skeletal form in Figures 2A,B for (A) optimal and (B) salt-stress treatments, respectively; the below-ground part of the stem was false-colored in green, the upper roots in light brown, and the lower roots in dark brown, with a slight glossy effect. See text for the procedures of expansion from the skeletal images and attachment of corn root volumes to the skeletons, using the CT scanning data collected and applying appropriate thresholds to them.

one of them in Week 2. Combining upper roots, when present, with lower roots in the same fractal analysis would not be recommended anyway, because the two types of roots form distinct systems by themselves (**Figures 2** and **3**). Overall, the observed FD values are low to moderately high, that is, above 1.0 (which classically corresponds to a straight line) but well below 2.0 (a plane), and actually below 1.5; this could be expected for young plant material such as root systems of corn seedlings only a few weeks old. On a weekly basis (mean ± SE), FD estimates ranged from 1.084 ± 0.026 (Week 1) to 1.284 ± 0.037 (Week 3) in the control group vs. from 1.132 ± 0.006 (Week 1) to 1.279 ± 0.030 (Week 3) in the salt-stressed group, with 1.231 ± 0.044 vs. 1.188 ± 0.045 in Week 2 for the two groups, respectively. The crossing of curves between the two groups from Week 1–2 in **Figure 4A**, which is the result of a higher FD mean value in Week 1 for the salt-stressed group and a higher FD mean value in Week 2 for the control group, led us to analyze the ratios of FD estimates in Weeks 2 and 3 relative to the corresponding FD estimate in Week 1 (**Figure 4B**), and eventually, perform two ANOVARs instead of one.

The ANOVARs indicate non-significant (*P* ≥ 0.10) betweengroup differences in the FD mean values for the 3 weeks taken separately as well as in the two FDR mean values, considering one ratio at a time. However, on average over the 3 weeks (see Between-subject effects in **Table 1**), the between-group difference in the FD overall mean value is significant (*P <* 0.10). Additionally, the Week main effects (see Within-subject effects in **Table 1**) are highly significant (*<sup>P</sup> <sup>&</sup>lt;* 0.01) for the FDs, and are only significant (*P <* 0.10) for the FDR. To be complete, the difference in the FD mean value is significant (*P <* 0.05) between Weeks 1 and 2, and highly significant (*P <* 0.01) between Weeks 1 and 3.

#### ANOVARs on Root Volumes and Relative Growth Rates

As shown in **Figure 4C**, the weekly means of lower-root volume for the control group are systematically above those for the salt-stressed group; furthermore, this is also the case for the weekly means of upper-root volume, due to the delayed onset of upper roots for salt-stressed plants (see "Three-Dimensional Image Analyses" and **Figures 2** and **3**).

FIGURE 4 | Means and SE of corn root system traits measured for optimal (Control) and salt-stress (Salt) treatment groups: (A) fractal dimensions (FD) in Weeks 1–3 (only for lower-root systems; see Table 1); (B) fractal dimension ratios (FDR; i.e., against Week 1); (C) root volumes in Weeks 1–3 (for lower- and upper-root systems); (D) relative growth rates (i.e., relative to the previous week and calculated using volumes of lower and upper roots combined); (E) root lengths in Weeks 1–3 (for lower- and upper-root systems); and (F) length increments (i.e., between successive weeks and calculated using root lengths of lower- and upper-root systems combined).

TABLE 1 | Repeated measures analysis of variance (ANOVAR) results for fractal dimensions (FD) and fractal dimension ratios (FDR) for corn lower-root systems.


†*The Error df is 4 (df: number of degrees of freedom).*

Relative growth rates of lower-root volume between Weeks 1 and 2 are almost equal in the two experimental groups; for all the other relative growth rates that we considered, the mean value for the control group is greater than that for the salt-stressed group, with a difference of 10% or more (**Figure 4D**).

Aside from the obvious difference in means between the two groups for the upper-root volume in Week 1, the ANOVAR performed on root volumes (**Table 2**) provides: (i) significant (*P <* 0.10) between-group differences for the three types of root volumes considered (lower roots alone, upper roots alone, lower and upper roots combined) in Week 3; (ii) significant (*P <* 0.10) Group main effects (Between-subject effects) for



†*Except for N/A, the Error df is 4 (df: number of degrees of freedom).* ‡*Corn plants of the salt-stressed group had no upper root in Week 1. It follows: (i) the non-relevance of performing a statistical test to compare groups for their upper-root volumes in Week 1; (ii) within-subject effects with 1, 1, and 4 df instead of 2, 2, and 8 df, due to the inclusion of 2 weeks (Weeks 2 and 3) instead of the 3 weeks in the ANOVAR for upper-root volume.*

§ *0.3776 if a t-test with an effective number of df is used, to adjust for the sole case of significant (P < 0.05) between-group heterogeneity of the variance in our study.*

the upper-root volume averaged over Weeks 2 and 3; (iii) highly significant (*P <* 0.01) Week main effects (Within-subject effects) for the three types of root volumes; and (iv) a significant (*P <* 0.10) or highly significant (*P <* 0.01) Week-by-Group interaction, indicating an increasing difference in weekly means between groups from the first or second week to the last. As for relative growth rates (**Table 3**), the ANOVAR only found highly significant (*P <* 0.01) Time main effects (Withinsubject effects) for the upper-root volume and the lower and upper-root volumes combined (i.e., no significant Week-by-Group interaction); thus, relatively, space occupancy by the root systems increased at the same pace in the two groups of corn seedlings. To be complete, the ANOVA found a significant (*P <* 0.10) difference between groups in their mean relative growth rates for upper-root volume from Week 2 to Week 3.

#### ANOVARs on Root Lengths and Length Increments

Graphically, the weekly means of lower- and upper-root lengths and the corresponding increments varied over time and differed between the control and salt-stressed groups, in a way similar to the weekly means of lower- and upper-root volumes and the relative growth rates (see **Figures 4E,F** vs. **Figures 4C,D**). Quantitatively and statistically, a smaller number of significant effects were found for lengths than for volumes, with two significant (*P <* 0.10) effects per type of root when analyzed separately, and combined; nevertheless, five of the six Withinsubject effects maintained their statistical significance. On length increments, significant differences between groups are found for lower roots (Weeks 1–2, *P <* 0.10) and upper roots (Weeks 2–3, *P <* 0.01).

TABLE 3 | Repeated measures analysis of variance results for relative growth rates derived from volumes for corn lower- and upper-root systems and the two combined.


†*Except for N/A, the Error df is 4 (df: number of degrees of freedom).* ‡*Corn plants of the salt-stressed group had no upper root in Week 1, so relative growth rates based on upper-root volume could be calculated only from Weeks 2–3 for them, limiting the comparison between groups to a classical type of ANOVA (i.e., without repeated measures).*


TABLE 4 | Repeated measures analysis of variance results for length measures for corn lower- and upper-root systems and the two combined.

†*Except for N/A, the Error df is 4 (df: number of degrees of freedom).* ‡*Corn plants of the salt-stressed group had no upper root in Week 1. It follows: (i) the non-relevance of performing a statistical test to compare groups for their upper-root lengths in Week 1; (ii) within-subject effects with 1, 1, and 4 df instead of 2, 2, and 8 df, due to the inclusion of 2 weeks (Weeks 2 and 3) instead of the 3 weeks in the ANOVAR for upper-root length.*

TABLE 5 | Repeated measures analysis of variance results for length increments for corn lower- and upper-root systems and the two combined.


†*Except for N/A, the Error df is 4 (df: number of degrees of freedom).* ‡*Corn plants of the salt-stressed group had no upper root in Week 1, so length increments based on upper-root length could be calculated only from Weeks 2–3 for them, limiting the comparison between groups to a classical type of ANOVA (i.e., without repeated measures).*

# Discussion

### Experimental Approach

Traditionally, studies of root system architecture under laboratory conditions mostly use platforms such as WinRhizo that can help generate data for plants grown in solid medium. For example, a study conducted on a number of *Arabidopsis* accessions, using EZ-Rhizo to screen for root system architecture related to a quantitative trait locus, suggests natural variations across the accessions (Armengaud et al., 2009). Although such an approach is non-invasive, plants grown in Petri plates for such analyses have a critical time frame of study, since the enclosed environment can eventually lead to stressful conditions for plant growth

Frontiers in Plant Science | www.frontiersin.org April 2015 | Volume 6 | Article 256 |

within the Petri plates. Other techniques to screen root system architecture, such as the use of PVC pipes (Shelden et al., 2013) and gel chambers (Bengough et al., 2004) for fast screening of seedlings, or of a Rhizobox or Rhizotron facility (with PVC boxes of various sizes) to study growth of visible roots along the transparent sides of a box, are also useful for short-term studies (Neumann et al., 2009). However, these studies are restricted to 2-D scanning and cannot predict the 3-D growth of roots.

Our study is one of a few in which root systems and their surrounding soil medium have been repeatedly CT scanned (e.g., Lontoc-Roy et al., 2005; Han et al., 2008, 2009). Repeated CT scanning of plant structures in successive stages has multiple advantages, including those of following the development of the same structure over time and capturing the changes. It also has constraints and limits, starting with the use of low X-ray doses because a root system is living plant material. As pointed out by Dutilleul et al. (2005), the radiation output of a CT scanner increases strongly with tube voltage, but it is the product of tube current (mA) by scan time (s) by number of scans that actually represents the radiation level delivered during exposure time. Helical scanning reduces exposure time substantially, while allowing the reconstruction of several images from CT scanning data acquired in one rotation around the sample installed on the couch. In root system studies involving CT scanning, since the plant structure of interest is generally surrounded by a soil medium (i.e., an exception is provided by hydroponic culture) it is equally important that the X-ray dose be sufficiently high to allow radiation to penetrate the soil column throughout. In fact, CT numbers are function of X-ray attenuation coefficients (Eq. 1), and too-low doses would result in imprecise CT numbers for the part of the CT scanned volume near the center of the column in the case of a soil as dense as sand. All in all, some balance must be found, such as the use of 120 kV as tube voltage and 150 mA as tube current for plastic pots with an 18-cm top diameter in our study. Working with smaller containers may be an option to consider depending on the plant species, but when chosen, this option can be at the expense of the last stage up to which the development of root systems can be followed appropriately, since root tips might reach the edges of small containers more rapidly which would alter root growth.

The soil type × moisture content combination is another very important factor to take into account in the repeated CT scanning for underground plant structures. Lontoc-Roy et al. (2006) have shown that sieved homogeneous sand allows a better isolation of corn roots from the plant–soil CT scanning data when it is relatively dry, than when it is water-saturated at the time of CT scanning; for loamy sand (by mass: 78.0% sand, 14.4% silt, and 7.6% clay), it is the contrary. Furthermore, it is recommended that soil moisture content be as much as possible the same at all the times of CT scanning for all the plants, to keep as much as possible the same plant–soil contrasts in the CT scanning data and images and obtain comparable results; this was achieved by watering the corn seedlings after CT scanning, leaving the sand dry before CT scanning, in our experiment. Otherwise, some data transformation can accommodate the situation (Han et al., 2008, 2009). The constraints above are important, but are readily eclipsed by the tremendous insight

gained by the follow-up made possible week after week, thanks to the repeated CT scanning, on developing root systems and their non-destructive characterization, with no such extra variability as the one that would be introduced if different seedlings were used for analysis at different times. Of course, as discussed in Subsection "Analytical Aspects," the repeated-measures nature of FDs and root volumes, for example, must be incorporated in the statistical tests, for them to be valid, in a study like ours.

#### Analytical Aspects

Three types of analytical aspects are discussed here: (i) graphical, in relation to spatial resolution, and (ii) statistical, concerning (ii.1) the recommended assessment procedure for effects involving time, and (ii.2) the question of power of the tests concerning treatment effects (e.g., optimal vs. salt stress). Our choice of CT scanning configuration parameter values (field of view: SS or 18 cm, zoom factor: 1.0; see Section "Materials and Methods") provided a spatial resolution sufficiently fine to isolate from plant–soil CT scanning data and reconstruct all primary roots and portions of secondary roots for the corn seedlings which were 1, 2, and 3 weeks old in our experiment; the voxel dimensions were 0.35 mm × 0.35 mm × 0.4 mm. This would be coarse for smaller root systems, like wheat and rice (Flavel et al., 2012; Zappala et al., 2013); in these two root CT scanning studies, the voxel dimensions were, in fact, smaller. Thus, everything is relative; for a high-resolution CT scanner such as the Toshiba XVision, the scale of observation is in dm and the scale of resolution is 0.1–1.0 mm, while for a ultra-high-resolution CT scanner, they are in cm and 0.01– 0.1 mm, respectively (see Table 1 in Ketcham and Carlson, 2001).

The repeated-measures nature of the characteristics (e.g., FDs, root volumes) derived from 3-D images of reconstructed plant structures in a repeated CT scanning study like ours has implications. The presence of temporal heterogeneity of variance and temporal autocorrelation, two common properties of temporal repeated measures, contradicts the ANOVA assumptions of homogeneity of variance and independence, and the classical ANOVA *F* tests may be invalid. A multiplicative factor (called "Box's epsilon"; Crowder and Hand, 1990; Dutilleul, 2011) is then applied to the numbers of degrees of freedom of the *F* distribution of the ANOVA test statistic, providing the ANOVAR *F* tests. This generally results in an increase of the *P*-values, thus decreasing the statistical significance of the effects related to the temporal repeated measures (e.g., Week and Week × Group). The modification was slight with three temporal repeated measures in our case, because autocorrelation and heteroscedasticity were then weak, and it is never required with only two repeated measures (i.e., our FDR and relative growth rates).

Now that the experimental approach is well established (see Subsection "Experimental Approach"), the CT scanning of larger numbers of plants per day can be envisaged. Besides a better representation of the treatment effect in time, larger sample sizes (i.e., numbers of individuals in a group) will enhance the power of the statistical tests, meaning that existing effects will be detected more often. From the results reported in **Tables 1–3** and calculations of degrees of freedom made aside, it can be seen that with a few more units in each group, most of the *P*-values between 0.10 and 0.20 in **Tables 1–3** could become smaller than 0.10 or even 0.05, all other things being the same. This would mean a larger number of significant (*P <* 0.10 or *P <* 0.05) effects, especially those related to the treatment (e.g., salt stress), with greater insight into differences in structural complexity and space occupancy of the developing root systems. That said, the number of significant (*P <* 0.01, 0.05, or 0.10) treatment or time-related effects found with three plants per group is remarkable and very encouraging for larger experimentation.

#### Biological Significance

Corn salinity tolerance studies have most often focused on the physiological, biochemical, phytohormone, transcriptional, and proteomic responses and in comparison with model plants *Arabidopsis*, rice, and tomato. Salinity responses vary between and within plant organs, growth stages, and are genotype-specific. Normally, roots can differentiate between wet surroundings and air pockets in its environment, and the collective response is referred to as hydropatterning. The immediate environment of the root dictates root hair patterning, positioning, and development of aerenchyma and distribution of anthocyanins and auxin, and is independent of abscisic acid signaling. In *Arabidopsis*, genes such as TRYPTOPHAN AMINOTRANSFERASE OF *ARABIDOPSIS* 1 and PIN-FORMED3 are necessary to control auxin to induce lateral root formation under high water availability. This also dictates the position where the lateral root founder cells need to be positioned and activated (Baoa et al., 2014). However, when a 100-mM NaCl stress is imposed in corn, a 10-fold increase in abscisic acid (ABA) accumulation in roots, as compared to a onefold increase in shoots, has been observed (Jia et al., 2002). Salt-resistant hybrid SR03 was found to have increased indole-butyric acid (IBA) and ABA levels in the shoots, while the roots maintained increased indoleacetic acid (IAA) levels upon 100-mM NaCl stress (Zörb et al., 2013). Quantitative differences in responses to salt stress were observed in salt-sensitive corn cultivar Trihybrid 321 and salttolerant cultivar Giza 2. Salt stress decreased fresh weight, dry weight, and relative growth rates of both shoots and roots. An increased accumulation of Na+ and Cl− and a decrease in K<sup>+</sup> and Ca2<sup>+</sup> were observed in both shoots and roots in cultivar Giza 2 (Mansour et al., 2005). A study on the root growth direction of *Arabidopsis* to salt stress suggests that the roots encountering salt stress have reduced gravity response, and this seems to be modulated by the quantity of salt present in its environment (Sun et al., 2008). The phenomenon of decreased root growth is evident in our study, graphically (**Figures 2** and **3**) as well as quantitatively (**Figure 4**), and statistically (**Tables 1–5**). Salinity stress has an early effect during corn seedling establishment that is more pronounced in the volume aspects than in the structural complexity of the root system.

Given the above understanding of salinity stress on roots, and shoots, of corn, the CT scanning results that we presented add a new dimension to the understanding of root architecture patterns in growing corn seedlings, unstressed, and in salinity stress condition. This root repeated CT scanning experiment was to quantify and illustrate *in situ* the effects of salinity on germinating and growing corn seedlings at optimal growth temperature, as salinity stress is a very simple abiotic stress to simulate under laboratory conditions.

As possible future perspectives or studies, we can see and propose what follows on the basis of the results obtained in our study. Furthering this study to root growth patterns under cold stress will be very insightful because cold stress delays the onset of roots and the growth pattern is slower as compared to optimal conditions, thereby allowing CT scanning for two or three additional weeks using the experimental approach discussed in Subsection "Experimental Approach." Slowing down the root growth can add to the refinement of the study of root branching patterns, volumes and lengths, and derived growth rates and increments. Such a study could be extended to other agriculturally important

## References


crops and their commercial genotypes, especially crops cultivated in temperate regions, with larger sample sizes.

#### Acknowledgments

The authors are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for the Major Equipment Grant that led to the creation of the CT Scanning Laboratory for agricultural and environmental research on Macdonald Campus of McGill University and to the Canadian federal government and the Québec provincial government for a Canada Foundation for Innovation (CFI) Grant that provided us with the computer equipment. The operating aspects of the reported study were supported through the Regroupements stratégiques Grant from the Fonds de recherche du Québec – Nature et technologies (FQRNT) to Centre SÈVE, of which PD and DS are members.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Subramanian, Han, Dutilleul and Smith. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Effects of damping-off caused by *Rhizoctonia solani* anastomosis group 2-1 on roots of wheat and oil seed rape quantified using X-ray Computed Tomography and real-time PCR

*Craig J. Sturrock1, James Woodhall2, Matthew Brown1, Catherine Walker1, Sacha J. Mooney1 and Rumiana V. Ray1\**

*<sup>1</sup> School of Biosciences, University of Nottingham, Sutton Bonington Campus, Loughborough, UK, <sup>2</sup> The Food and Environment Research Agency, Sand Hutton, UK*

*Rhizoctonia solani* is a plant pathogenic fungus that causes significant establishment and yield losses to several important food crops globally. This is the first application of high resolution X-ray micro Computed Tomography (X-ray μCT) and real-time PCR to study host–pathogen interactions *in situ* and elucidate the mechanism of *Rhizoctonia* damping-off disease over a 6-day period caused by *R. solani*, anastomosis group (AG) 2-1 in wheat (*Triticum aestivum* cv. Gallant) and oil seed rape (OSR, *Brassica napus* cv. Marinka). Temporal, non-destructive analysis of root system architectures was performed using RooTrak and validated by the destructive method of root washing. Disease was assessed visually and related to pathogen DNA quantification in soil using real-time PCR. *R. solani* AG2-1 at similar initial DNA concentrations in soil was capable of causing significant damage to the developing root systems of both wheat and OSR. Disease caused reductions in primary root number, root volume, root surface area, and convex hull which were affected less in the monocotyledonous host. Wheat was more tolerant to the pathogen, exhibited fewer symptoms and developed more complex root systems. In contrast, *R. solani* caused earlier damage and maceration of the taproot of the dicot, OSR. Disease severity was related to pathogen DNA accumulation in soil only for OSR, however, reductions in root traits were significantly associated with both disease and pathogen DNA. The method offers the first steps in advancing current understanding of soil-borne pathogen behavior *in situ* at the pore scale, which may lead to the development of mitigation measures to combat disease influence in the field.

Keywords: *Rhizoctonia solani*, X-ray Computed Tomography, qPCR, wheat, oil seed rape, fungi, soil

*Edited by: Pierre Dutilleul, McGill University, Canada*

#### *Reviewed by:*

*Pierre Dutilleul, McGill University, Canada Timothy Paulitz, U.S. Department of Agriculture – Agricultural Research Service, USA Norman Gerhard Uhlmann, Fraunhofer Institute for Integrated Circuits IIS, Germany*

#### *\*Correspondence:*

*Rumiana V. Ray, School of Biosciences, University of Nottingham, Sutton Bonington Campus, Loughborough, Leicestershire, LE12 5RD, UK rumiana.ray@nottingham.ac.uk*

#### *Specialty section:*

*This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science*

*Received: 31 October 2014 Accepted: 10 June 2015 Published: 24 June 2015*

#### *Citation:*

*Sturrock CJ, Woodhall J, Brown M, Walker C, Mooney SJ and Ray RV (2015) Effects of damping-off caused by Rhizoctonia solani anastomosis group 2-1 on roots of wheat and oil seed rape quantified using X-ray Computed Tomography and real-time PCR. Front. Plant Sci. 6:461. doi: 10.3389/fpls.2015.00461*

# Introduction

*Rhizoctonia solani* Kühn (teleomorph = *Thanatephorus cucumeris* Donk) is a ubiquitous soil-borne plant pathogenic fungus which causes significant yield losses in many agriculturally important crops (Verma, 1996; Paulitz et al., 2006). Individual isolates of *R. solani* are classified into anastomosis groups (AGs) based on their hyphal incompatibility and their host specificity (Anderson, 1982). For example, AG2-1 and AG4 are associated with stem and root rot diseases in dicotyledonous crop species belonging to *Brassicaceae* (Gugel et al., 1987; Sneh et al., 1991; Tewoldemedhin et al., 2006) whilst isolates of AG8 cause 'bare patch' or root rot on monocotyledonous crops from *Poaceae* (Paulitz et al., 2002).

The predominant population of *R. solani* causing severe seedling diseases associated with establishment losses of up to 80–100% and final yield loss of up to 30% of oil seed rape (OSR, *Brassica napus*) worldwide belongs to AG2-1 (Tahvonen et al., 1984; Kataria and Verma, 1992; Khangura et al., 1999). Highly virulent isolates of AG2-1 cause pre- and post-emergence damping-off, stem, and root rot with characteristic water soaked lesions on the root and hypocotyl, stunting of plant growth, root necrosis and cortex tissue maceration, and subsequent death in OSR (Yang et al., 1992). Recent soil surveys, carried out in USA (Schroeder et al., 2011) and UK (Brown et al., 2014) on fields growing winter wheat (*Triticum aestivum*) have revealed the most common pathogen present in soils of increased rotational frequency with OSR is *R. solani* AG2-1, shown in *>*69% of fields (*n* = 90) in England.

Whilst the pathogenicity and aggressiveness of AG2-1 to OSR have been previously studied (Yitbarek et al., 1987; Kranz, 1988), less is known of the impact of this group of pathogens on wheat roots. AG2-1 isolates have been shown to be pathogenic to cereals to varying degrees. Tewoldemedhin et al., (2006) reported AG2- 1 isolates were weakly pathogenic to barley and wheat roots. In contrast, Roberts and Sivasithamparam (1986) reported AG2-1 isolates from wheat roots in 'bare patch' in Western Australia were highly pathogenic to wheat causing an 80% disease index which was similar to disease caused by AG8 isolates. Thus, at present, the ability of AG2-1 to cause significant damage to the root system of seedlings of monocotyledonous crops such as wheat remains unclear.

The etiology of soil-borne diseases caused by pathogens such as *R. solani* on plant seeds and roots below ground has until recently been difficult to study. Traditionally, assessment of disease incidence and severity has involved the use of visual observations of symptoms of infection on affected plant organs following the physical extraction of plants from the ground (Kranz, 1988). However, the inherently destructive nature of visual disease inspection means that it is not possible to monitor temporal disease development and effects on root traits and system architecture. Furthermore, destructive sampling in the field often results in an incomplete root system extraction and loss of the most severely infected or severed primary/secondary roots.

Non-destructive methods for imaging plant roots *in situ* in soil, such as X-ray micro Computed Tomography (X-ray μCT), have become an important tool for quantifying plant root system architecture development in three dimensions (see review by Mooney et al., 2012). However, to date the application of X-ray μCT to investigate the impact of root rot pathogens has been relatively limited to Han et al. (2008) who studied the effects of common potato scab caused by *Streptomyces scabies* on tubers in soil. This was the first use of medical X-ray CT in a phytopathological study to successfully segment root structures from CT images and demonstrated diseased plants had significantly less complex root systems, in addition to delayed root growth and branching. A subsequent study by the same researchers using CT showed the effects of common potato scab on the density of seed and peripheral organs of potato plants in soil over a 10 week period (Han et al., 2009). Interestingly, an early application of a medical CT system to soil science by Grose et al. (1996) measured moisture content in bulk soil and in the soil around roots to predict suitable growth conditions for both *R. solani* and *Gaeumannomyces graminis*. Although at relatively coarse resolutions (200 μm) compared to the resolution achievable on modern systems for similar sized pots (6 cm diameter), the study successfully quantified heterogeneous moisture gradient in the vicinity of the plant roots and demonstrated the potential of the technique for investigation of environmental factors on the soil–plant–microbe system. Recent advances in the sensitivity of X-ray detectors within industrial μCT systems have facilitated much faster acquisition times (minutes rather than hours) facilitating easier repeated scanning of the same sample to visualize the temporal dynamics of plant root systems in undisturbed soil a (Tracy et al., 2012; Zappala et al., 2013b)

Microbiological methods for detection and quantification of target AGs of *R. solani* in soil are highly labor intensive and time consuming, involving the use of soil baiting methods that are often inefficient in detecting and isolating *R. solani*, and microscopy (Sneh et al., 1991). Furthermore, low population densities of *R. solani* in the soil and the lack of selective isolation media for the species make quantification difficult and unreliable. In the last decade, several conventional or real-time quantitative polymerase chain reaction (qPCR) assays have become an established tool for rapidly quantifying fungal pathogens including targeted AGs of *R. solani* at low detection limits in both soil and infected plant tissues (Filion et al., 2003a,b; Sayler and Yang, 2007; Okubara et al., 2008; Budge et al., 2009; Woodhall et al., 2013). We propose that the combination of these two powerful techniques, qPCR and X ray μCT, can allow improved new insight into the temporal host–pathogen interactions and provide quantitative data on the impact of soilborne pathogens on root architectural systems of crop plants grown in soil. The main aim of this study was to elucidate the mechanism of disease caused by AG2-1 of *R. solani* on root traits and system architecture of two different crops, the monocot, wheat, and the dicot, OSR.

# Materials and Methods

#### Soil, Plant, and Inoculum Preparation

The experiment was designed as a factorial block with two main factors, host and inoculation with two levels. The host crops were wheat, (*Triticum aestivum* cv. Gallant) or OSR (*Brassica napus* cv. Marinka) which were either non-inoculated or inoculated with *R. solani* AG2-1 (Isolate 159/8, Goll et al., 2014). The isolate was previously determined to be weakly pathogenic to wheat and pathogenic to OSR. There were nine replicates of the treatment combinations resulting in a total of 36 columns.

Soil columns (30 mm diameter × 70 mm length) were uniformly packed to a bulk density of 1.1 Mg m−<sup>3</sup> with a Newport series loamy sand soil (sand 72.6%, silt 13.2%, and clay 14.2%; pH 6.35; organic matter 2.93%) collected from the University of Nottingham farm at Bunny, Nottinghamshire, UK (52.52◦N, 1.07◦W). Prior to packing, the soil was air-dried, sieved to *<*2 mm and sterilized by γ-irradiation at 27 kGy (Isotron, Daventry, UK). The pathogen treated soils were inoculated with five, 5-mm diameter plugs of actively growing *R. solani* mycelium equally distributed in the vertical direction of the soil during packing of the columns. Seeds of cv. Gallant and cv. Marinka were pre-germinated for 48 h on moist filter paper in petri dishes before being planted at 10 and 5 mm below the soil surface, respectively. The columns were then saturated, drained for 2 days (to a notional field capacity which represents the moisture content of the soil after free drainage had ceased) and placed in a growth room under conditions of 14◦C day/night with an 8 h photoperiod and a photosynthetic photon flux density (PPFD) at plant level of 1000 μmol m−<sup>2</sup> s−1. A transparent plastic unheated seed propagator was used to maintain high relative humidity levels and avoid surface drying of the soil during seedling establishment in the growth room. Three replicates for each treatment combination were randomly selected and destructively harvested via root washing and scored for disease 2, 4, and 6 days following inoculation (dfi). Root disease severity was assessed at each destructive sampling point on soil-free plants on scales from 0 to 5; 0 = no lesions, clean roots; 1 = small lesion on tap root; 2 = necrosis of upto 30%; 3 = necrosis covering 31–60% of the tap root; 4 = necrosis covering 61–99% of the tap root; 5 = completely severed tap root (Khangura et al., 1999). In addition, the three replicates selected for harvest at 6 dfi were also scanned using X-ray μCT at 2, 4, and 6 dfi to permit non-destructive quantification of root system development. Root architecture of the washed roots was assessed using WinRHIZO<sup>R</sup> 2002c scanning equipment and software on each harvest day. The images collected were used to compare with the X-ray μCT images. Soil from the columns was further used for DNA extraction and pathogen quantification.

### X-Ray Micro Computed Tomography (**µ**CT)

The replicate subset allocated for destructive sampling at 6 dfi (12 columns), were scanned at 2, 4, and 6 days using a Phoenix Nanotom<sup>R</sup> (GE Measurement & Control Solutions, Wunstorf, Germany) X-ray μCT scanner. The scanner consists of a 180 kV nanofocus X-ray tube fitted with a tungsten transmission target and a 5-megapixel (2304 × 2304 pixels, 50 × 50 μm pixel size) flat panel detector (Hamamatsu Photonics KK, Shizuoka, Japan). A maximum X-ray energy of 110 kV, 140 μA current and a 0.15 mm thick copper filter was used to scan each sample which consisted of 1300 projection images acquired over a 360◦rotation. Each projection image was the average of three images acquired with a detector exposure time of 500 ms in 'Fast CT mode.' The resulting isotropic voxel edge length was 19 μm (i.e., spatial resolution) and total scan time was 35 min. The total X-ray dose for each sample was calculated as 25.2 Gy over the three scans, which is below the 33 Gy threshold reported by Johnson (1936) which no detrimental effects of postgermination plant growth following exposure to X-ray radiation were observed (Zappala et al., 2013a). Reconstruction of the projection images was performed using the software datos| rec (GE Measurement & Control Solutions, Wunstorf, Germany) to produce 3-D volumetric data sets with dimension 30 × 30 mm (diameter × depth).

#### Image Processing and Analysis

Plant root systems were non-destructively segmented using the *Region Growing* selection tool in VG StudioMAX<sup>R</sup> 2.2 software as described by Tracy et al. (2012). To summarize, the *Region Growing* tool, allows the user to select connected structures within the data that have the same distribution of X-ray attenuation based on their gray values. The user assigns all root material to a region of interest which is then extracted as a separate binary image stack for measurement of root system architecture in RooTrak software. RooTrak software (Mairhofer et al., 2012) permits quantification of descriptive traits on root system architecture, such as total volume, surface area, maximum length and width, convex hull (relates to the space filling in 3D of an object), and centroid Z (relates to the center of mass of a 3D object). Due to small scales differences in seed depth in the reconstructed volumetric data, the measurement field of view was standardized to 30 mm × 25.80 mm (diameter × depth). Therefore, the maximum possible value for root length measurements is limited to 25.80 mm.

Soil porosity (total and incremental with depth) was quantified in FIJI image analysis software (Schindelin et al., 2012) using a modified method of Tracy et al. (2012). To summarize, a resized 16 bit image stack of dimensions 17.1 mm × 17.1 mm × 19 mm (900 × 900 pixels × 1000 images) was first prepared to exclude the area outside of the soil column (i.e., the container and the surrounding air space). Images were binarised to define the air filled pore space with a value of 0 and the 'solid' soil with a value of 1 using the isodata threshold algorithm which performed the best in an evaluation study. Soil porosity for each slice image was calculated based on the percentage of air to the total volume of the resized stack.

### Real Time Quantitative PCR for AG2-1 of *R. solani*

DNA was extracted from soil as described in Woodhall et al. (2012), except sample size was reduced to 45 g and then added to a 250 ml Nalgene bottle with 3 ml antifoam B with six 25.4 mm stainless steel ball bearings and 90 ml grinding buffer (120 mM sodium phosphate buffer pH 8, 2% cetrimonium bromide, 1.5 M sodium chloride). Real-time PCR was undertaken using a 7500 real-time PCR system. Environmental Master Mix 2.0 (Life Technologies, USA) was used for all real-time PCR and consisted of half the total reaction volume of 25 μl, whilst 5 μl consisted of the DNA sample. Primers (MWG Biotech, Germany) and hydrolysis probe specific for AG2-1 (Budge et al., 2009) were used and added to a final concentration in the reaction of 300 and 100 nM, respectively, with the remaining volume made up with molecular grade water. Cycling conditions consisted of 50◦C for 2 min, 95◦C for 10 min, and 40 cycles of 95◦C for 15 s, and 60◦C for 1 min. Each sample was tested in duplicate and an average Ctvalue was determined. Target DNA in soil samples was quantified by including six DNA standards on each PCR run. The standards consisted of a DNA sample of known concentration taken from culture of AG2-1 (Isolate 2023, Food and Environment Research Agency, UK) which was used to produce a dilution series of five 10-fold dilutions. The amount of DNA was then determined by linear regression.

#### Statistical Analysis

Root growth and architecture traits were analyzed using analysis of variance (ANOVA) for repeated measures and corrected for degrees of freedom for all time related effects with Greenhouse– Geisser Epsilon factor. Architecture traits were root volume, surface area, convex hull volume, maximum width, and length. Pathogen DNA data were analyzed by ANOVA containing sampling time, crop, and inoculation as interacting factors in the treatment structure. Regression analysis was used to investigate the relationships between root traits, disease score, and pathogen DNA, using a simple linear model for each crop separately. All analyses were performed in Genstat 15, version 15.1.0.8035.

#### Results

#### Disease Development and Pathogen DNA Accumulation in Soil

No symptoms of root disease were observed in the noninoculated treatments (control) for either crop species (**Figure 1**). OSR plants developed visible lesions on roots as soon as 2 dfi. The symptoms rapidly progressed from moderate (necrosis covering 31–60% of the root, disease score 3) to severe (completely severed taproot, disease score 5) by 4 dfi resulting in complete maceration of root tissue by day 6 (**Figure 1**). Wheat plants exhibited significantly lower disease severity compared to OSR plants (*P* = 0.011), with symptoms classified as slight (small lesions on the primary roots, disease score 1) which were first detected at 6 dfi (**Figure 1**).

DNA of *R. solani*was not detected in the soil of non-inoculated plants at 2 dfi, but was quantifiable at 4 and 6 dfi at low concentrations (0.008 and 0.019 ng g−1) in two soil columns. In contrast, DNA in inoculated soils of both crops at 2 dfi was above 100 ng g−<sup>1</sup> (**Figure 2**). The trend of DNA accumulation over the duration of the sampling period was similar for the two crops showing an increase in pathogen DNA by day 4 followed by a plateau by 6 dfi (**Figure 2**). The mean pathogen DNA in the OSR treatment at 4 dfi was approximately 45% higher than in the wheat treatment (*P* = 0.063) although no differences were observed between crops for 2 or 6 dfi.

lesion on tap root; 2 **=** necrosis of up to 30%; 3 **=** necrosis covering 30–60% of the tap root; 4 **=** necrosis covering 61–99% of the tap root; 5 **=** completely severed tap root) assessed 2 days following inoculation (dfi), 4 and 6 dfi on wheat and oil seed rape

No disease symptoms were shown in the control treatment for both crops. Bar shows standard error of difference (SED) for the interaction between sample time (T) at 2, 4, or 6 dfi and crop (C) species (wheat or OSR).

## Impact of *R. solani* AG2-1 on Root System Architecture of Wheat and OSR

Visual assessment of X-ray μCT 3D images and WinRHIZO<sup>R</sup> images suggested major differences in root system architecture under the experimental factors, inoculation and crop (**Figure 3**; Supplementary Videos 1 and 2). Control OSR plants had a characteristic single tap root that developed lateral roots by 6 dfi. Typically, wheat plants developed between 3 and 5 primary roots with no lateral roots by the end of the experiment. Initial root

FIGURE 4 | Root system volume and surface area over time (T) for crop (C), (A,B) and inoculation (I) with *R. solani* AG2-1 (C,D). Interactions for surface area and volume were detected using repeated measures ANOVA with degrees of freedom (df) corrected by Greenhouse–Geisser epsilon factor. Bars show SED for (1) comparing means for treatment combinations; (2) comparing means with the same level of C; (3) for comparing means for the same level of I and species (wheat or OSR); (4) comparing means with the same level of I.

growth of OSR plants was inhibited in soils inoculated with AG2- 1 of *R. solani* and resulted in complete maceration of root tissue by 6 dfi. Disease effects were less obvious on wheat roots from inoculated soils with *R. solani* (**Figure 3**).

There were significant temporal differences for root volume and surface area measured using X-ray μCT between crops (**Figures 4A,B**; *<sup>P</sup> <sup>&</sup>lt;* 0.001) and between inoculated and noninoculated plants (**Figures 4C,D**; *<sup>P</sup> <sup>&</sup>lt;* 0.001). The absence of interactions between crop and inoculation suggested root volume and surface area were affected mainly by intrinsic differences in root system characteristics of individual crop species and the presence of the pathogen in the soils. Inoculation significantly reduced root volume and surface area in both crops, however, the effects were greater in OSR, where these traits were affected immediately following inoculation and there were relatively small changes over time in trait parameters (**Figure 4**).

Root system traits for which significant temporal interactions between crop and inoculation were detected are shown in **Table 1**. The root system of wheat increased in length and width in time, despite inoculation, to a maximum of 25.8 and 29.3 mm, respectively (**Table 1**). A similar trend was observed for the control OSR plants with the root system length and width reaching 25.8 and 13.5 mm, respectively, by the end of the experiment. However, for the OSR plants inoculated with *R. solani*, root growth was inhibited from day 2, slight increases in length and width were observed by day 4 but ultimately at 6 dfi roots of inoculated plants were 96% shorter and 78% thinner than the controls (0.97 and 2.90 mm, respectively).

Both inoculation treatments in wheat displayed a significant increase in centroid Z (an indication of root structure with depth) after 4 days incubation with a mean value of 16.04 and 16.47 mm for the control and inoculated plants, which then reduced to 14.86 and 14.5 mm, respectively, after 6 dfi. Control OSR plants displayed a sustained increase in centroid Z from 1.07 mm at 2 dfi to 18.52 mm at 6 dfi. Centroid Z remained consistently low throughout the experiment for the *R. solani* treated OSR plants (1 mm). (**Table 1**; time <sup>×</sup> crop <sup>×</sup> inoculation; *<sup>P</sup>* <sup>=</sup> 0.010.)

Convex hull (an indication of the volume of soil explored) increased in all treatments except in OSR inoculated plants, where it remained the same after 4 dfi and for wheat was significantly higher compared to OSR (*P* = 0.001). Inoculation with *R. solani* resulted in smaller rates of increase in convex hull in both plants (**Table 1**). The control wheat treatment showed a significantly higher convex hull which was almost twice the volume compared to the *R. solani* inoculated treatment with values of 4123 and 2038 mm3, respectively, after 6 dfi. The control OSR had a lower convex hull compared to wheat with a mean of 413 mm3. *R. solani* treated OSR exhibited the lowest convex hull with a mean of 7 mm<sup>3</sup> remaining the same at 4 and 6 dfi (**Table 1**; time <sup>×</sup> crop <sup>×</sup> inoculation; *P* = 0.048).

Inoculation with *R. solani* AG2-1 had a major effect on primary root number of both crops and resulted in significant reductions throughout the experiment demonstrated by the absence of significant interactions between experimental factors and time (**Figure 5A**). The number of primary roots was significantly higher in wheat compared to OSR plants which produced just one taproot (**Figure 5B**). Production of primary root numbers in wheat ceased at 4 dfi with no further significant increases being detected (**Figure 5B**). In OSR plants primary root number decreased at each sample time associated with effects of inhibition by the pathogen on root development and digestion of root tissue in time (**Figure 5B**).

Comparison of the WinRHIZO<sup>R</sup> and RooTrak measurements supported all observations and displayed strong significant relationships for comparable root system traits such as volume (*P <* 0.001, *R*<sup>2</sup> = 0.97) and surface area (*P <* 0.001, *R*<sup>2</sup> = 0.97). The relationship for root length measured by the two methods was also significant (*P* = 0.024) but weaker than previously mentioned traits accounting for only 39% of the variance.

#### Relationship of Pathogen DNA and Root System Traits

Linear regression analysis with groups for individual crops was carried out to test the fitted data for the measured traits, pathogen DNA and visual disease symptoms for position and parallelism (**Table 2**). There was a significant relationship between disease score and pathogen DNA accounting for 82% of the variance, however, the data fitted separate lines for each crop, with different slope and intercept indicating a positive relationship between pathogen DNA in soil and disease expression on plant roots for OSR only. Data fitted separate lines for each crop for root length measured by <sup>μ</sup>CT on both disease (*<sup>P</sup> <sup>&</sup>lt;* 0.001, *<sup>R</sup>*<sup>2</sup> <sup>=</sup> 0.96)



and pathogen DNA (*<sup>P</sup> <sup>&</sup>lt;* 0.001, *<sup>R</sup>*<sup>2</sup> <sup>=</sup> 0.77) with the same directionality showing negative relationships (**Table 2**). Similarly regressions (*P <* 0.001) of surface area and root length, measured by WinRHIZO<sup>R</sup> , on disease score accounted for more than 96% of the variance. Fitted separate lines with the same directionality for wheat and OSR suggested that the magnitude of effects on developing traits of the different root systems of individual crops were related to the expression of disease symptoms. All other measured traits by different systems fitted parallel lines for disease expression indicated that the final effects were similar but dependant on intrinsic differences between crops (**Table 2**).

#### Analysis of Soil Porosity

Total mean soil porosity, limited to an extent by the spatial resolution of the scans, was consistent for all soil columns across all treatments (Mean, 15.4%, SEM 1.5). However, measurement of the porosity with depth within a column showed regions of variable porosity indicative of layering created during soil packing which varied between 8 and 50% (**Figure 6C**). Furthermore, there was evidence of higher porosity at the interface of the emerging seedling and the surrounding soil in some of the samples, where the highest porosity values of 50% were recorded. This was particularly evident in one of the OSR replicates treated with *R. solani* AG2-1 showing hypocotyl tissue maceration and decay in the area of high soil porosity (**Figures 6D,E**; Supplementary Video 3). However, there was only

of OSR by *R. solani* AG2-1. (A) 3D X-ray CT image of soil and root (yellow). (B) Image showing only root tissue (white solid arrow indicates maceration of tissue. (C) 2D cross-section (zx plane) image showing high porosity around OSR root (scale bar = 2 mm). (D) Magnified view of image shown in (C), showing necrosis of root cortex (scale bar = 1 mm). (E) 2D cross-section (xy plane) image showing preservation of the stele (solid arrow) but complete necrosis of cortex tissue (scale bar = 0.5 mm).

weak regression between DNA concentration and soil porosity (*R*<sup>2</sup> <sup>=</sup> 0.21).

## Discussion

This work provides the first example of X-ray μCT used for the non-destructive detection of below ground symptoms and impact of *R. solani* on the developing root systems of monocotyledonous and dicotyledonous plants. *R. solani* AG2-1 causes significant pre- and post-emergence damping-off characterized by the inhibition of seed germination, root elongation, and ultimately the digestion of the root and hypocotyl of *Brassica* species (Kataria and Verma, 1992). We found moderate symptoms in OSR as early as 2 dfi and severe disease developed by 4 dfi. In contrast, only mild symptoms developed in wheat plants by 6 dfi for similar initial inoculum in the soil quantified using qPCR as pathogen DNA at 2 dfi. The difference in disease development and severity on the two crops is in agreement with previous reports on the virulence and aggressiveness of AG2-1 to OSR demonstrating that isolates belonging to this group are highly pathogenic to *Brassica* species (Gugel et al., 1987; Verma, 1996).


TABLE 2 | Linear regression models for disease score (y) on pathogen DNA (x) and WinRHIZO<sup>R</sup> , and X-ray CT based measurements of root system architecture traits (y) on disease score (x), and pathogen DNA (x) for each crop.

*\*Denotes values not applicable.*

The delay in symptom development on wheat suggests that AG2- 1 is unable to cause significant symptoms on wheat confirmed by others in their investigations of pathogenicity of *R. solani* AG2-1 to cereals (Khangura et al., 1999; Oros et al., 2013). The effect of the primary host crop, OSR, on *R. solani* development was evident in the more rapid increase of pathogen DNA, reaching maximum of 300 ng g−<sup>1</sup> in soil by 4 dfi in contrast to a twofold less DNA in soils from wheat grown plants (data not shown). This fast DNA accumulation in the soil from OSR, compared to wheat, is most likely related to the differences in the rate of infection and digestion of the emerging radicle and hypocotyl of the primary host species, manifested by the numerous lesions (visualized in this study) inhibited growth and ultimately the complete seedling necrosis by 6 dfi. The plateau of soil pathogen DNA at 6 adfi may be due to an exhaustion of available nutrients from the host plants and return of the pathogen to saprophytic phase of survival. The temporal dynamics of the pathogen during the development of wheat or OSR in field rotations are currently unknown. However, Brown et al. (2014) found no significant differences in pathogen DNA of *R. solani* AG2-1 accumulation in English field soils of wheat following wheat or OSR, suggesting that short wheat/OSR rotations are unlikely to be effective in reducing inoculum concentrations for either crop.

Visualization of the 3-D root system of the two crops grown in soil showed how the contrasting root systems of the monocot and dicot species reacted to the pathogen infection. Differences in the impact of the pathogen appeared to be related to the intrinsic complexity of the architectural root systems of the two crops and their ability to compensate on specific traits. Using time series μCT data importantly revealed that although the infection in the monocot, wheat, appeared asymptomatic, it contrasted the severe symptom expression in the dicot, OSR. *R. solani* AG2-1 was capable of causing significant damage on important developing root architectural traits of both crops including primary root number, root volume and root surface area that were affected less in the monocotyledonous host. Furthermore, the ability of both hosts to explore soil via their developing root system, indicated by the convex hull, was reduced. However, traits such as root length and centroid Z were not affected in the monocot. Both inoculated and control wheat plants developed 3– 4 primary roots that were thicker and longer by 4 dfi compared to OSR plants. In contrast, OSR plants were mostly dependent on the development of strong taproot and subsequent lateral roots for the acquisition of resources, thus early damage to the developing taproot by *R. solani* diminished significantly the ability of the plant to establish or recover from the disease. Wheat was able to compensate by producing more than one primary root (seminal roots) and it is likely that uninfected or less severely infected roots by the pathogen were able to escape the disease and thus compensate for resource use. *R. solani* AG2-1 is most aggressive to young seedlings and host resistance to infection increases with age (Verma, 1996). Therefore faster developing OSR cultivars are more likely to escape the disease and traits related to early germination and establishment, such as seed size will be important for breeding new varieties that are more likely to tolerate *R. solani* infection (Hwang et al., 2014).

Disease score and pathogen DNA were both strongly related to changes in the measured root traits. However, the transiency of these effects in particular in the maturing wheat plant is unknown. The relationship between disease and pathogen DNA was different for the two crops and disease was only predicted successfully for OSR. This has implications in terms of assessment and prediction of disease in the field in relation to individual crop species as clear symptoms were not exhibited in wheat and not related to DNA concentrations. Furthermore, both crops suffered from *R. solani* at the seedling stage thus it is important to elucidate if the disease caused by AG2- 1 is associated with significant yield loss of wheat in the field. Understanding the relationships between initial inoculum concentrations and final yield loss for the two crops can assist in the development of new strategies for prediction of risk and yield loss based on qPCR of soil prior to planting.

From the measured root traits, only root length showed poor correlation between the two imaging approaches which can be attributed to the way the trait was measured. RooTrak root length measurements were limited to a maximum soil depth of 25.80 mm compared to the entire 30 mm column length due to the field of view possible in a single μCT scan. However, as RooTrak can quantify novel root traits such as convex hull, there is potential to measure crop species specific descriptors to define root structure, e.g., differences between the single tap root of OSR versus primary and seminal root system of wheat. A crucial advantage of the μCT imaging is that not only can the developing root systems be quantified non-destructively and temporally but as we have shown changes in the soil microstructure can also be considered. Although, our initial soil conditions were designed as in most repacked column studies to be uniform, verification of the microstructure by imaging showed localized variations in porosity when measured at high resolutions especially at the root surface. This zone, i.e., the rhizosphere, is a crucial interface, where knowledge about the structural arrangement in particular is lacking. Variations in structure as we have revealed here will influence soil moisture availability considering the relationship between matric suction and pore size. Soil bulk density and moisture content are known to significantly influence hyphal growth and disease severity caused by *R. solani* (Glenn et al., 1987; Gill et al., 2001, 2004) but the impact at the pore scale is less well understood. Furthermore, it is generally accepted that the key limiting factor in hyphal proliferation is the availability of air filled pores within the soil (Glenn et al., 1987; Otten et al., 1999; Harris et al., 2003). We found OSR seedlings displaying the highest porosity around the seedling also had the lowest disease severity and longest root and shoot growth (**Figure 6**). This finding is in agreement with Gill et al. (2000) who found that although saprotrophic growth was higher in more porous soils, the disease severity was lower highlighting the potential of X-ray μCT in the study of the physical effects of soil structure on soil-borne pathogenic fungal diseases. This has potential implications for soil management practices, such as conventional and zero tillage as these may

### References

Anderson, N. A. (1982). The genetics and pathology of *Rhizoctonia solani*. *Annu. Rev. Phytopathol.* 20, 329–347. doi: 10.1146/annurev.py.20.090182. 001553

have very different soil structures (Mangalassery et al., 2014). For example, plowing could potentially reduce soil-borne disease severity, by increasing the porous structure of soil, physical disruption to fungal hyphal networks and increasing background microbial activity. Indeed, the most effective cultural control method for soil-borne *Rhizoctonia* root patch in wheat is via tillage practice of soil disturbance by cultivation which destroys established fungal hyphal networks and can increase microbial activity (Paulitz et al., 2002). The effect of tillage on soil-borne pathogens in OSR has received less attention, however, it is likely that reduced or zero tillage maximizes disease and inoculum potential by allowing infected crop residues to remain on the soil surface and preserving hyphal networks in close proximity to the host (Kharbanda and Tewari, 1996). Although soil structure can routinely be imaged at high resolutions (i.e., *<*100 μm), it is still not possible to visualize fungi *per se* using X-ray μCT due to their very low X-ray attenuation (Gleason et al., 2012). However, indirect modeling approaches have been useful to aid understanding of the behavior and functioning of fungi in both real (Pajor et al., 2010; Falconer et al., 2011) and artificial soil microstructures (Otten et al., 2012). These combined approaches may be of value in the future to facilitate further understanding of plant pathogenic fungi in the soil environment.

This study has successfully quantified the impact of *R. solani* on crop root system traits and development through the combined use of X-ray μCT and qPCR. X-ray μCT offers more promise than destructive methods as the development of disease symptoms on the root can be monitored non-destructively in soil. We have shown that disease symptoms developed rapidly in OSR within 2 dfi, whereas wheat displayed a higher tolerance with only mild symptoms present after 6 dfi. Differences in the impact of the pathogen on the two hosts were related to complexity and developmental rates of the different root architectural types of the monocot, wheat, and the dicot, OSR.

### Acknowledgments

The authors would like to acknowledge support from the University of Nottingham and Syngenta for the provision of AG2-1 isolate. JW received support from the Defra Horizon Scanning and Technology Implementation Fund. M. Brown receives a Ph.D. studentship from Syngenta. The authors thank Jim Craigon for his helpful advice on data analysis.

# Supplementary Material

The Supplementary Material for this article can be found online at: http://journal*.*frontiersin*.*org/article/10*.*3389/fpls*.*2015*.*00461


in soil. *Plant Pathol.* 58, 1071–1080. doi: 10.1111/j.1365-3059.2009. 02139.x


in tomato (*Solanum lycopersicum*) by X-ray Micro-Computed Tomography. *Ann. Bot.* 110, 511–519. doi: 10.1093/aob/mcs031


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Sturrock, Woodhall, Brown, Walker, Mooney and Ray. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# X-ray computed tomography uncovers root–root interactions: quantifying spatial relationships between interacting root systems in three dimensions

Alexander M. Paya<sup>1</sup> , Jesse L. Silverberg<sup>2</sup> , Jennifer Padgett <sup>3</sup> and Taryn L. Bauerle<sup>1</sup> \*

*<sup>1</sup> School of Integrative Plant Science, Cornell University, Ithaca, NY, USA, <sup>2</sup> Department of Physics, Cornell University, Ithaca, NY, USA, <sup>3</sup> School of Electrical and Computer Engineering, Cornell University, Ithaca, NY, USA*

#### Edited by:

*Pierre Dutilleul, McGill University, Canada*

#### Reviewed by:

*Lars Hendrik Wegner, Karlsruhe Institute of Technology, Germany Pierre Dutilleul, McGill University, Canada Keith Jones, Brookhaven National laboratory, USA*

#### \*Correspondence:

*Taryn L. Bauerle, School of Integrative Plant Science, Cornell University, 134A Plant Science Bldg., Ithaca 14853, NY, USA bauerle@cornell.edu*

#### Specialty section:

*This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science*

> Received: *24 June 2014* Accepted: *06 April 2015* Published: *29 April 2015*

#### Citation:

*Paya AM, Silverberg JL, Padgett J and Bauerle TL (2015) X-ray computed tomography uncovers root–root interactions: quantifying spatial relationships between interacting root systems in three dimensions. Front. Plant Sci. 6:274. doi: 10.3389/fpls.2015.00274*

Research in the field of plant biology has recently demonstrated that inter- and intra-specific interactions belowground can dramatically alter root growth. Our aim was to answer questions related to the effect of inter- vs. intra-specific interactions on the growth and utilization of undisturbed space by fine roots within three dimensions (3D) using micro X-ray computed tomography. To achieve this, *Populus tremuloides* (quaking aspen) and *Picea mariana* (black spruce) seedlings were planted into containers as either solitary individuals, or inter-/intra-specific pairs, allowed to grow for 2 months, and 3D metrics developed in order to quantify their use of belowground space. In both aspen and spruce, inter-specific root interactions produced a shift in the vertical distribution of the root system volume, and deepened the average position of root tips when compared to intra-specifically growing seedlings. Inter-specific interactions also increased the minimum distance between root tips belonging to the same root system. There was no effect of belowground interactions on the radial distribution of roots, or the directionality of lateral root growth for either species. In conclusion, we found that significant differences were observed more often when comparing controls (solitary individuals) and paired seedlings (inter- or intra-specific), than when comparing inter- and intra-specifically growing seedlings. This would indicate that competition between neighboring seedlings was more responsible for shifting fine root growth in both species than was neighbor identity. However, significant inter- vs. intra-specific differences were observed, which further emphasizes the importance of biological interactions in competition studies.

Keywords: micro-computed tomography, interspecific interactions, belowground competition, Picea mariana, Populus tremuloides, 3D root system images

# Introduction

Plants sharing a finite amount of space will inevitably interact with each other either above or belowground in the pursuit of essential resources. The outcomes of these interactions can range from positive (facilitation) to negative (competition), and are therefore highly relevant for the development of agricultural and ecological management practices (Grime, 1979; Tilman, 1987). Traditional parameters that quantify the effect of belowground interactions on root growth dynamics include diameter class, spatial/temporal deployment, growth rate, and fine root abundance (Casper and Jackson, 1997; Eissenstat and Yanai, 1997; Eissenstat et al., 2000; Hodge, 2004, 2009; Kembel et al., 2008). While parameters such as these differ across species, accurate in situ observations are inherently limited by the opaque and heterogeneous nature of soil matrices, and generally require a destructive harvest of roots (Joslin and Henderson, 1984; Steingrobe et al., 2000), or visualization along a two dimensional (2D) surface (Gross et al., 1993; Majdi, 1996; Eissenstat et al., 2000).

However, recent advances in three dimensional (3D) imaging technology such as ground penetrating radar, laser imaging, nuclear magnetic resonance imaging (MRI), neutron radiography (NT), and X-ray computed tomography (CT) have made the observation of undisturbed root systems possible (Macfall et al., 1991; Butnor et al., 2001; Gregory et al., 2003; Kaestner et al., 2006; Perret et al., 2007; Tracy et al., 2010; Moradi et al., 2011; Mairhofer et al., 2012). Innovations in software such as Rootviz, Root track (Tracy et al., 2010; Mairhofer et al., 2012), RootReader3D (Clark et al., 2011), and Avizo (Saoirse et al., 2010), and specific filtering algorithms (Perret et al., 2007) have helped improve 3D image resolution and stream-line the quantification of anatomical parameters such as lateral root length, lateral root number, root-system surface area, and volume of undisturbed root systems. However, accurately isolating roots from root-soil data is complicated by the continuum of water within the root itself, at the root-soil interface, and between soil particles (Lontoc-Roy et al., 2006). As methods for isolating roots improve, steady technological advancements will, and have already increased the scope of viable research questions and objectives. For example, studies have already begun to explore the 3D spatial distribution of fine and coarse roots in forests (Pierret et al., 1999; Butnor et al., 2001), mechanical buckling in plant roots (Silverberg et al., 2012), and water uptake at the root-soil interface (Moradi et al., 2011).

As 3D imaging technologies become more widely available, questions about the occupation and exploration of space by interacting root systems can be better addressed, offering new insights to this important yet problematic component of root growth dynamics. For example, belowground interactions can result in whole root system segregation (reviewed in Schenk et al., 1999), stunted root elongation (Mahall and Callaway, 1996; Falik et al., 2005; Bhatt et al., 2011), and/or over-yielding in response to spatially proximal self-roots (belonging to same plant) and non-self roots (belonging to neighbor) (Gersani et al., 2001; Maina et al., 2002; Falik et al., 2003). An understanding of the mechanisms regulating the growth of roots driven by belowground interactions is still developing, however growing evidence suggests that traditional parameters including root biomass, root surface area, and diameter are insufficient in integrating spatially complex responses.

To our knowledge, the following experiment is the first attempt at observing and quantifying the effect of belowground interactions between two neighboring root systems in 3D. Our research employed micro-CT to capture the spatial distributions of both interacting (inter- vs. intra-specific) and control (solitary) root systems belonging to 2-month-old tree seedlings. 3D models of root system architecture were developed from annotated CT image slices, and traditional belowground parameters such as root length, surface area, volume, and number of root tip were either measured or counted. Moreover, we also developed a series of belowground parameters that take advantage of skeletonized 3D root systems and binary root system data, and evaluate distances between root tips and characterized the distribution of individual root volumes. Broadly, the goal of this work was to investigate the effects of inter- and intra-specific interactions on belowground parameter values, and compare these parameter values with those obtained from solitary individuals.

# Materials and Methods

# Plant Growth

Acrylic tubes (3.5 mm wall thickness, 64 mm inner diameter, 305 mm length) were covered with fine mesh (0.5 mm) along the base, capped, and secured to provide free drainage. Each tube was filled incrementally with polystyrene beads (1–3 mm), gently tamped throughout the filling process in order to reduce pore size and achieve greater bulk density, and then wrapped in aluminum foil to prevent light penetration. Polystyrene was used in place of peat, sandy loam, sand, or vermiculite based on trail experiments which demonstrated very high water retention in soil or soil like mediums. Additionally, contrast agents such as iodine containing compounds, barium sulfate (BaSO4), gold chloride (Au2Cl6), and cow's milk were used, but sufficient contrast was not achieved.

Black spruce (Picea mariana) and quaking aspen (Populus tremuloides) were selected as "interacting" plant species based on differences in phylogeny, morphology, and the fact that they co-occur across northern latitudes of North America (DeByle and Winokur, 1985). Seeds from each species were germinated for 5–7 days between two sheets of damp cellulose, and then transplanted into pre-wet hydroponic growth plugs (Rapid Rooster Grow plug™, General hydroponics, Sebastopol, CA) following radicle emergence (1–3 mm). Each acrylic tube received a single plug containing one individual of either spruce or aspen (control), or two plugs, each containing one individual, to simulate inter- or intra-specific interactions. A total of 25 tubes were prepared. There was five of each of the following tubes: solitary aspen, intra-specific aspen, interspecific aspen/spruce, intra-specific spruce, and solitary spruce. Containers were randomly arranged on a hydroponic flood table modified to re-circulate nutrient solution for top-down irrigation. To prevent competition for light, a sheet of acetate was placed between interacting seedlings.

Plants were grown under greenhouse conditions (17 C◦ night, 20◦C/day; KPL greenhouses, Cornell University, Ithaca NY) with supplemental lighting (12 h/days) for 60 days (April–June, 2012). Irrigation was fed from a 530 liter (140 gallon standard) reservoir to individual micro-spray emitters focused at the base of each plant (90<sup>0</sup> , 0.5 gph, Hydro Flow™, Redmond, WA). Irrigation provided each tree with a nutrient replete hydroponic solution balanced at 150 ppm N (Peter salts, 21-5-20, 382.8 g/530 L; Epsom salts, MgSO<sup>4</sup> 7H2O, 130.64 g/530 L). Hydroponic Paya et al. X-ray vision

nutrient solution was maintained at a pH of 5.5–5.8, and electro conductivity of 1.6–2.0 throughout the experiment, and automatically controlled by a programmable timer on a rotating schedule. Half way through the growth phase of the experiment, a pump malfunction resulted in the loss of multiple individuals, which produced an uneven number of replicates for each species.

Irrigation was terminated after 2 months of growth. Plants were allowed to transpire residual water remaining in each container for 2 days prior to imaging in order to reduce imaging artifacts. Plants were then transported to Cornell's imaging facility for CT scanning.

#### Micro-CT Scans

Whole seedling's root systems were imaged at Cornell University's Micro-CT facility. Due to a pump failure during the growth stage of the experiment, only 13 out of 25 tubes were imaged: 3 × solitary aspen, 1 × intra-specific aspen, 3 × inter-specific aspen/spruce, 3 × intra-specific spruce, and 3 × solitary spruce. Each scan was performed using a GE CT120 micro-CT scanner (GE Healthcare, London, ON, Canada). Initially, 10 bright-field images were acquired with no objects in the scanner, providing a correction for detector non-uniformity. Calibration and correction for signal non-uniformity was determined from measurement within a SB3 (GE Healthcare) water/bone phantom, scanned with the samples. Resulting image datasets were calibrated to the conventional scale of Hounsfiel radiodensity units (HU), defined so that water and air have HU values of 0 and -1000, respectively. Each scan digitally acquired 720 projections at 0.5◦ intervals over 360◦ using 80 keV, 32 ma, 32 ms exposure time and 100µm x-y-z resolution. The obtained projections were used to reconstruct a CT dataset using a convolution back-projection algorithm implemented in 3D, giving a 70 × 70 × 50 mm<sup>3</sup> volume of image data with 100µm isotropic voxels.

Using a sequential stacking function (MicroView, GE Healthcare), three sequential image stacks (70 × 70 × 50 mm<sup>3</sup> ) were taken from each tube and recompiled into a single 70×70 × 150 mm<sup>3</sup> data set. Using this function, we successfully increased the visible volume (i.e., visible rooting structure) three-fold: from 5 to 15 cm depth (70×70 × 150 mm<sup>3</sup> scan required 1 h of imaging time).

### Destructive Harvesting

Following X-ray scanning, plants were destructively harvested. Leaves/needles and petioles were removed from the main stem and scanned using a photo scanner (Epson Expression 10000XL, 2400 dpi, Epson America Inc., Long Beach, CA). Directly following the removal of aboveground tissues, acrylic containers were inverted and tamped to release the polystyrene medium along with roots, which were gently rinsed under a 0.5 mm sieve. Polystyrene beads still attached to roots were removed using forceps. Individual roots were separated manually to prevent overlapping segments, placed on a photo scanner, and scanned. After scanning, above and belowground tissues were placed in separate paper bags, dried at 55 C for 3 days, and then weighed. Scanned images were analyzed for leaf surface area, root surface area, and total root length using WinRhizo (Winrhizo 2011, Regent Instruments, Canada). The number of root tips were counted manually using the image analysis toolset, ImageJ (National Institute of Health, Bethesda MD). In ImageJ, each 2D root system scan was imported as a TIFF, and using the paintbrush tool, root tips were individually and sequentially numbered to ensure that no root tips were overlooked or counted twice.

#### Image Reconstruction

Projections were exported from the GE CT120 micro-CT scanner as VFF format (Sun TAAC Graphic File) and converted to DICOM format using MicroView's DICOM transfer tool. Image stacks were then imported one at a time into ImageJ using the import, image sequence function. Once the main taproot was found for each seedling, roots originating from a single individual seedling were given an arbitrary color code, and the entire cross sectional area was traced by hand for each root through each CT image slice (70 × 70 × 0.1 mm) (**Figure 1**). This ensured that root diameter, surface area, and volume could be measured in the 3D model/reconstruction, and that root systems belonging to different individuals could be easily differentiated. Once a root came in contact with the container wall, tracing ceased because, (1) the roots were indiscernible from the container, and (2) the roots behaved atypically and tracked the container wall. Color-coded image stacks were then exported as an RGB TIFF stack, and opened in MATLAB <sup>R</sup> 2012b for three dimensional reconstruction and spatial quantification (The MathWorks Inc., Natick MA). In order to 3D render each root system, color codes were identified and isolated, which allowed for the subtraction of non-colored voxels; annotated circles representing root cross sections within each x-y plane were then stacked across the zplane. This process effectively rendered each root system in 3D with little or no constraints on actual root system dimensions. In MATLAB each root cross section was also converted to a 3D binary matrix in order to measure spatially explicit parameters. Entries of these matrices were either 0 or 1, depending on whether that voxel was occupied by the root.

## Morphological and Anatomical Analyses

Root surface area was determined from the 3D data sets by sequentially analyzing each x-y cross-section with MATLAB's bwtraceboundary function. This identified the coordinates of the root perimeter from which we calculated the circumference of all roots passing through the plane. The circumference was multiplied by the cross-sectional thickness (100µm) to estimate root surface area per image slice. This was performed for all crosssectional images and the results summed to calculate root system surface area. Root system volume was calculated by summing the total number of occupied voxels and multiplying by the volume per voxel, 10−<sup>3</sup> mm<sup>3</sup> /voxel.

Root tips were located by scanning through each crosssection and identifying terminal voxels. This generated a list of coordinates that were centered on root tips (**Figure 1**). With this information we were able to determine the depth-wise distribution along the z-axis (complete 3D information). To quantify the radial distribution of root tips (mm<sup>2</sup> ), root tip coordinate data was projected into the x-y plane. An ellipse

whose circumference and orientation represents the occupied x-y area of the root system was then projected over the root tip coordinates. Multiplying π by both the major and minor axes of the ellipse provided the radial distribution of each root system. The ratio of the ellipse's major to minor axis is then a metric that defines how radially symmetric the root distribution is. In particular, if the ratio is 1, then the distribution is circularly symmetric. Values higher than 1 indicate the amount of asymmetric root growth in the plane that passes vertically through the ellipse's major axis.

#### Statistics

In order to validate our 3D rendering protocol, we used simple linear regression to compare destructively harvested (2D) and 3D reconstructed root system data. Surface area and the number of root tips were chosen for regression because both were measurable in 2D and 3D, and could therefore be used to validate our manual tracing procedure of fine root cross sections through CT image slices. The intercept (a) in the simple linear regression model (y = a + bx + error) was constrained to be equal to 0, and thus was not estimated (see **Figure 3**). Differences in mean ranks between solitary, intra- and interspecifically growing plants were analyzed using the Kruskal-Wallis test on the following parameters: biomass, 2D/3D root length, 2D/3D surface area, specific root area (SRA), specific root length (SRL), manual and 3D root tip count, 3D root volume, root–root distance, major/minor axes, radial distribution, and rooting depth weighted by volume (P = 0.05, H<sup>o</sup> = mean ranks are equal). Where the null was rejected, post-hoc analyses were performed using the Wilcoxon each-pair test (non-parametric multiple comparisons, P = 0.01667). It is important to note that the distribution of minimum Euclidean distances between root tips did not distribute normally, nor did the data behave normally post-transformations (e.g., log x, e<sup>x</sup> , or x−<sup>1</sup> ). Instead, we used a two-parameter model for the probability distribution function that fit both species data where the ratio of the sum of squares of regression (SSreg) to the total sum of squares (SStot) was between 0.75 and 0.99. This model was a good predictor of the fraction of root tips f(x) dx separated by a minimum distance between x and x + dx for control, inter-, and intra-specifically growing seedlings, where the pdf is given by:

$$f(\mathbf{x}) = c\_1 \mathbf{x} e^{-\mathbf{x}/c\_2}.\tag{1}$$

Non-linear regression models were fitted using Sigma Plot 11 (Systat Software, San Jose, CA). All other statistics were done using JMP 10.0 (SAS Institute Inc., Cary, NC).

### Results

#### Destructive Analyses

Analyses of each harvested seedling showed that belowground interactions had no significant effect on the aboveground growth of aspen or spruce (P = 0.1487, **Figures 2A,B**). Belowground interactions had no significant effect on either species' root biomass (P = 0.0606, **Figure 2C**). Belowground interactions had a measurable effect on aspen's root surface area, but not spruce (P = 0.0439, **Figure 2D**). On average, inter-specific aspen root systems were reduced to 31% of the control samples surface

Two additional belowground parameters, SRL (cm mg−<sup>1</sup> ) and SRA (cm<sup>2</sup> mg−<sup>1</sup> ), were also calculated for both species. SRL and SRA depict the cost of root construction, and can be highly informative in establishing whether a treatment had an effect on root morphology. Control, inter-, intra-specifically growing spruce seedlings differed in terms of SRA (P = 0.0389), but similar differences were not observed in aspen (P = 0.0783, **Table 1**). Control, inter-, and intra-specific interactions also had a significant effect on the SRL of spruce seedlings (P = 0.0134), but not aspen (P = 0.2082).

#### Validation of 3D Reconstruction

In order to validate our 3D reconstruction protocol, comparisons were made between destructive (i.e., 2D) and 3D root parameters. Surface area and the number of root tips were two parameters measured in both 2D and 3D, and were therefore chosen to validate the 3D reconstruction by way of simple linear regression. We found that 62% of the total number of root tips, and 76% of the total surface area were successfully captured during 3D image reconstruction (**Figure 3**). Examples of 3D root systems form each species combination are presented in **Figure 4**, as well as **Figure S1**.

#### 3D Utilization of Space

The occupation of 3D space by individual root tips, and whole root systems was quantified via a set of five metrics: radial distribution of root tips, directionality (major/minor axes of radial distribution ellipse), minimum root–root tip distance, root system volume as a function of depth, and the vertical position of root tips. The radial distribution of root tips measured the radial expanse of a root system (mm<sup>2</sup> ), i.e., the area of an ellipse that encompassed the x/y distribution of all root tips projected along the z-axis. The radial distribution of aspen roots ranged from 248 to 501, 500 to 522, and 212 to 325 mm<sup>2</sup> in the control, intra-specific, and inter-specific seedlings, respectively (**Table 2**). For spruce, the control, intra-specific, and interspecific seedlings ranged from 16 to 43, 9 to 112, and 36 to 314 mm<sup>2</sup> , respectively (**Table 2**). We found no significant effect of belowground interactions on the radial distribution of roots.

The directional growth of roots was measured by dividing the major (transverse) by the minor (conjugate) axes of an ellipse that encompassed the radial distribution of each root system. Using this approach, we could determine whether a root system was concentric (major/minor = 1), or skewed/directional (e.g., major/minor >> 1). As was the case with the radial distribution of root tips, belowground interactions had no significant effect on the ratio of major/minor axes for either species. However, notable species trends were observed. The major/minor axes of spruce root systems ranged from 1 to 10, with a mean of 4 across solitary and paired seedlings, indicating that spruce roots systems were relatively planar (**Table 2**). The major/minor axes of aspen root systems ranged from 1.1 to 1.6, with a mean of 1.2 across control, inter- and intra-specific seedlings, indicating that aspen root systems were evenly distributed in all directions relative to each root system's center of mass.

The third metric, minimum root–root tip distance, measured the minimum distance between a root tip (x1, y1, z1) and the

area. Lastly, belowground interactions had a significant effect on fine root lengths in aspen (P = 0.0439), but not spruce (P = 0.2120).

post-destructive harvest. Each point (triangle) represents data from a single individual displayed on a log scale (y-axis); \* denotes significant differences among control, intra-, and inter-specific seedlings (Kruskal-Wallis test,

*P* ≤ 0.05).


TABLE 1 | Belowground parameters of destructively harvested aspen (Populus tremuloides) and spruce (Picea mariana) seedlings grown under three experimental conditions: control, intra-specific, and inter-specific.

*Values listed are the median (maximum, minimum).*

nearest neighboring root tip (x2, y2, z2) for every terminal point in a root system. Belowground interactions had a significant effect on the minimum distance between root tips in both aspen (P = 0.0114) and spruce (P = 0.0002). Post-hoc analyses indicated that aspen roots grown intra-specifically (6.8 ± 0.28 mm) had wider distances between root tips compared to the controls (5.9 ± 0.13 mm) (P = 0.0025, **Figure 5B**). As for spruce, post-hoc analyses indicated that the minimum distances between root tips in controls (6.0 ± 0.78 mm) were significantly less than interspecific seedlings (11 ± 1.0 mm, P < 0.0001, **Figure 5D**). The minimum distances between root tips in intra-specific seedlings (7.8 ± 0.92 mm) were also significantly less than inter-specific seedlings (P = 0.0007, **Figure 5D**).

We also quantified differences between control, inter-, and intra-specifically growing plants in terms of their vertical placement of root system volume (**Figures 6A,B**). For both species, control, inter-, and intra-specific interactions were significant predictors of the depth of roots weighted by volume (Kruskal-Wallis test, P < 0.0001). In aspen, post-hoc analyses indicated that the largest difference was observed between control (deep) and intra-specific (shallow) root systems (P < 0.0001). The second largest difference in rooting depth was observed between inter- and intra-specifically growing seedlings (P < 0.0001). The smallest difference was observed between interspecific and control root systems (P < 0.0001). In spruce, posthoc analyses indicated that the largest difference was observed between inter- and intra-specifically growing root systems (P < 0.0001). Inter-specific and control root systems also occupied significantly different rooting depths (P < 0.0001), while no significant difference was observed between intra-specific and control root systems (P = 0.0864). The mean rooting depth for control, intra-specific, and inter-specifically growing root systems for aspen was 56.2 ± 0.6, 23.6 ± 0.6, and 32.7 ± 0.7 mm, respectively, whereas spruce was 41.0 ± 2.4, 34.4 ± 1.2, and 28.6 ± 2.31 mm, respectively.

The vertical distribution of aspen root tips differed significantly between control, intra-, and inter-specifically growing seedlings (P < 0.0001). Solitary aspen tended to distribute their root tips evenly across vertical space, and occupied a mean depth of 58.6 mm ± 1.43 mm. Intra-specifically growing aspen root tips were predominately located between 300 and 1200 mm, and occupied a mean depth of 72.5 ± 3.01 mm. Inter-specifically growing aspen root tips were concentrated between 400 and 1000 mm, and occupied a mean depth of 85.6 ± 2.41 mm. Post-hoc analyses indicated that inter-specific root tips were located significantly deeper than intra-specific root tips (P < 0.0001), as well as control root tips (P < 0.0001). Intra-specific root tips were also located significantly deeper than control root tips (P = 0.0017). The mean depths of spruce root tips did not differ between control, intra- and inter-specific

FIGURE 4 | 3D reconstruction of solitary and paired root systems (axes in 0.1 mm increments). (A) Control spruce (*Picea mariana*) (*n* = 3 containers, 3 individuals). (B) Intra-specific spruce (*n* = 3 containers, 6 individuals). (C) Inter-specific aspen (*Populus tremuloides*) and spruce root systems (*n* = 3 containers, 3 individuals each species). Aspen is highlighted in orange, spruce is highlighted in blue. (D) Intra-specific aspen (*n* = 1 container, 2 individuals). (E) Control aspen (*n* = 3 containers, 3 individuals). Note the differences in x and y-axes. For a complete list of 3D reconstructions, refer to Figure S1.

TABLE 2 | Measurements made in 3D of aspen (Populus tremuloides) and spruce (Picea mariana) seedlings grown under three experimental conditions: control, intra-specific, and inter-specific.


*Values listed are the median (maximum, minimum).*

seedlings, which were located at 45.2 ± 6.56, 50.4 ± 2.86, and 58.2 ± 6.10 mm, respectively.

### Discussion

#### The Effect of Belowground Interactions on Root Growth

In this experiment, each seedling was grown under full nutrient conditions without competition for light to minimize any variation that could be attributed to aboveground resource competition between solitary and paired individuals. The results of this study suggest that belowground interactions between neighboring root systems had measurable effects on root system architecture and 3D distribution of fine roots in the two species. For example in aspen, inter-specific interactions reduced belowground surface area, root length, and SRL when compared to solitary aspen seedlings, which suggests that inter-specific interactions negatively affect aspen root system growth (**Table 1**). In contrast, inter-specifically growing spruce seedlings, when compared to solitary spruce seedlings, increased root surface area, root length, SRA, and SRL (**Table 1**). This trend was not observed when comparing solitary and intra-specifically growing spruce seedlings, which suggests that spruce may over-yield under inter-specific growing conditions (Brassard et al., 2011).

We found that belowground competition shifted the 3D distribution of aspen roots. One example was root–root tip distances; the minimum distance between neighboring root tips under inter-specific conditions (7.5 mm) was greater than controls (5.9 mm, **Figure 5**). While a 1.6 mm difference in spacing may seem small, phosphorus concentrations increase exponentially over a distance of 1 mm from a root's surface (Hendriks et al., 1981). Therefore, relatively small adjustments in the spacing of root tips could result in distinct levels of competition for nutrients in soil (Hodge, 2009). The observed variation in root–root spacing may also result from a lower average number of root tips per inter-specific root system (102 tips) compared to intra-specific (173 tips) or controls (254 tips). Given a constrained volume, a reduction in the number of root tips could result in greater distances between them.

We also observed that the vertical placement of aspen's roots differed between the control, intra-, and inter-specific

seedlings (**Figures 6A,B**). This observation suggests that aspen can respond within 2 months of germination, and with high plasticity, to the presence of a neighboring root system by shifting the vertical placement of its roots. These changes in fine root vertical distribution may result in distinctly different rooting depths between aspen and inter-specific neighbors, thus reducing belowground competition (Schenk et al., 1999; Schenk, 2006). Alternatively, shifts in aspen's vertical root distribution may result in root system overlap with neighboring plants establishing direct competition. Regardless of whether aspen's strategy is to outcompete, avoid competition, or some combination of both in response to a neighboring plant, we found that the presence of both inter- and intra-specific neighbors is sufficient to alter the vertical growth of aspen roots. Care must be taken in the interpretation of these results, however, as specific changes in the occupancy of space by inter- and intra-specifically growing root

(Wilcoxon multiple comparisons, *P* ≤ 0.01667). Note the differences in axes.

FIGURE 6 | Heat map representing root system volume as a function of depth for aspen (A, Populus tremuloides) and spruce seedlings (B, Picea mariana). Control, intra-, and inter-specific interactions were significant predictors of rooting depth weighted by volume for both aspen and spruce (Kruskal-Wallis test, *P* < 0.0001). *Post-hoc* analyses indicated that in aspen, control, inter-, and intra-specifically growing seedlings were all significantly different in terms of their rooting depth weighted by volume (*P* < 0.0001). In spruce, *post-hoc* analyses indicated that inter-specific vs. control plants, as well as control vs. inter-specific plants were significantly different in terms of rooting depth weighted by volume (*P* < 0.0001). Units are in mm3. Each striated column represents the full root volume of a single seedling. Red arrows indicate the mean depth of root tips, while striped black lines indicate the mean depth of root system volume.

systems may result from their being planted side-by-side and not in the center of each container.

Belowground interactions also shifted the vertical positioning of aspen root tips, and to a lesser extent in spruce. In both solitary and paired aspen seedlings, the whole root system's mean rooting depth was shallower than the mean depth of root tips alone. Specifically, inter-specific aspen exhibited the greatest difference between mean root system depth and mean root tip depth (**Figure 6**). Under inter-specific conditions, the mean depth of root tips (85.6 mm) was markedly deeper than the mean depth of the whole root system (32.7 mm). This observable difference between the mean depth of an entire root system, and the mean depth of it's root tips is noteworthy, mainly because not all roots within a root system are functionally equivalent (Pregitzer et al., 1997). Especially in woody plants, a large proportion of a root system is in the form of long-lived, woody roots that anchor the plant and support essential transport functions, as opposed to the most distal part of the root system, the root tips, which are highly metabolically active and demonstrate the highest rates of nutrient and water uptake among all root classes (Pregitzer et al., 1997; Volder et al., 2005). By reducing overlap among "transport" roots and root tips, a root system can occupy an exclusive volume of soil space while simultaneously foraging for resources, and minimizing competition with itself. Therefore, when quantifying root growth dynamics in 3D volumes, either in response to itself or non-self interactions belowground, special attention should be paid to the dynamic growth and placement of root tips independently of whole root systems.

The root system architecture of spruce was dominated by a main taproot with relatively few lateral roots (**Figure 4**), which resulted in very few significant differences between treatments. However, there were some notable trends worth discussing. SRA and SRL tended to be higher under interspecific conditions compared to intra-specific conditions, as well as solitary individuals, indicating a lower cost of construction for spruce roots under inter-specific conditions. Also, inter-specific spruce roots tended to grow deeper, place root tips deeper, and root tips were spaced farther apart when compared to solitary plants—a response that was similar to aspen.

#### Modeling Root–Root Interactions

The use of mathematical models to describe biological phenomena is inherently complicated by the nature of organismal responses to heterogeneously distributed biotic and abiotic cues (Hodge, 2009), though in the context of root systems, both mechanical (Moulia, 2013), and fractal analysis (Tatsumi et al., 1989; Fitter and Stickland, 1992; Ozier-Lafontaine et al., 1999) have been applied with some success. Belowground, root system responses can manifest as a proliferation of roots into a nutrient rich patch (Robinson et al., 1999), altered root morphology (Bolte and Villanueva, 2006), or shifts in the direction of growth (Falik et al., 2005). Accurately modeling this type of non-random growth response is possible (e.g., Godin, 2000), but requires data that is highly resolved, both spatially and temporally. In this study, we demonstrate high spatial resolution for a single point in time, which limits our ability to quantify dynamic growth processes. Nevertheless, 3D structural information, such as that captured by the micro-CT techniques used here, provides insights otherwise inaccessible with 2D destructive imaging. We support that these findings can be used to verify or refute predictions of derived equations that incorporate the interactions between plants.

When modeling root–root distances, we discovered that a phenomenological exponential "growth and decay" model fit well for both species (SSreg/SStot = 0.75–0.99). This model was chosen from a number of mathematical models that were developed to accurately describe the distribution of root–root distances. Alternative models that were generated but not included in this study tended to fit the data somewhat better, but were species specific. This 2-parameter model was adopted because it accurately described the distribution of data for both species across all belowground conditions, whereas other similarly simple models could not. When developing and applying such models, it is important to keep in mind not only the model's fit, but the scope of the experimental question, the complexity of the model, its predictive value, and whether it can be applied broadly across species.

Also, based on equal amounts of nutrient, water, and light supplied to each container housing either one or two individual seedlings, and the equidistant orientation of paired individuals (intra- and inter-), we assumed that belowground interactions (intra- vs. inter-) would have the same effect on each paired seedling. There is no way for us to know with certainty that individuals were experiencing a treatment effect, or simply resource competition, which would result in similar belowground outcomes. However, as previously mentioned, nutrient levels were maintained at an EC of 1.6–2.0 throughout the experiment, reducing the likelihood of resource competition throughout the experiment. Future experiments should aim to parse out the different effects of nutrient competition and non-resource interactions.

#### Comparing Methods and Constraints

Previous attempts to image undisturbed root systems have been met with mixed success. The current benchmark for successfully rendering roots in a 3D data set is set at roughly 90%. For example, Gregory et al. (2003) captured 90% of 7-day-old wheat roots that did not exceed 10 cm in total length. Kaestner et al. (2006) reported that they could successfully capture 90% of Alnus incana (alder) roots larger than 0.18 mm. However, alder roots had to first be removed from their growing medium and packed into quartz sand prior to imaging. In another experiment, Perret et al. (2007) captured around 87% of the total root segments, and 78% of the total root length in 21-day-old chickpea.

In our experiment using P. mariana (spruce) and P. tremuloides (aspen), we successfully rendered between 62 and 76% of the actual root system architecture (**Figure 3**). We believe that roughly 30% of the root systems were lost in the manual root tracing/annotation phase of the methodology because of the criteria we followed for each annotation. Specifically, roots that contacted the container wall were excluded on the basis that these roots behave uncharacteristically, i.e., container tracking and circling. We predicted that based on aspen's larger root system size, a larger proportion of aspen roots would have been lost in the reconstruction phase when compared to spruce. However, data from both species were included in our general linear model (**Figure 3**). We found that both species realized a similar percentage loss of roots, which suggests that any bias introduced as a result of the annotation criteria was similar for both species.

Limitations of our methodology include (1) the use of synthetic growth medium, (2) manual identification of roots in the annotation process, and (3) the relatively small instrument aperture. Regretfully, we could not heed the call of Gregory and Hinsinger (1999), who argued that future advancements in research involving micro-CT and plant roots must focus on using natural soils in place of sand or hydroponics. Distinguishing between water within roots, and water in the medium is an oftenreported limitation—one we experienced early on during method development. We attempted to circumvent this issue by growing plants in hydrophobic, synthetic "sand." This way, the amount of water remaining in the container during imaging would be minimized, facilitating root identification. While this worked well for us, we cannot conclude that either species' growth was unaffected by this growth method. Though residual water was minimized, there were still trace amounts that disrupted several automated root-tracking algorithms that were developed during this experiment. Thus, the data required manually tracing each root through +1400 CT image slices, which required 6–8 h per dataset. Future plant research employing micro-CT should strongly consider developing a robust root-tracking approach that is insensitive to artifacts imposed by residual water in the growth medium.

Lastly, the instrument's aperture for accepting samples greatly limited our container size. Future work that employs micro-CT for phenotyping or quantifying belowground phenomena in an undisturbed space must consider the physical size constraints, and perhaps modify their experimental design to ensure that roots remain unimpeded by the boundaries of the container. The containers used for this experiment were sufficiently long (∼300 mm), but insufficiently wide (max ∼70 mm). Had the container width not been limiting, it is likely that a fewer roots would have been lost during 3D reconstruction.

# Conclusion

In this experiment, we could not conclude with any certainty that intra- and inter-specifically growing seedlings differed in terms of root system architecture and use of 3D space. We showed that, when compared to solitary individuals, interspecific interactions could have the effect of reducing root production, shifting the depth of root tips, increasing spacing between root tips, and altering the distribution of root system volume over vertical space. Because predictable shifts in rooting depths, lateral root placement, and/or root abundance based on neighbor identity may have far reaching implications in terms of ecosystem function (Hooper et al., 2005), species coexistence (Grime, 1997; Stoll and Prati, 2001; Bruno et al., 2003; Kembel et al., 2008), and plant evolution (Myers et al., 2000), interactions at the community level down to the individual and tissue level must be better understood. The future of this technique is in quantifying both very fine and coarse scale morphological and architectural shifts in root system growth. We demonstrate the ability to quantify 3D parameters and track multiple 3D root systems within a shared volume, which is an important advancement in the field of plant imaging. By coupling CT imaging with algorithms tailored to specific experimental conditions, a wide range of relevant architectural, morphological, and 3D parameters can be analyzed, and the effects of belowground interactions better understood. It is our aim that the marriage of CT with novel algorithms will continue to pave the way toward understanding how plants sense, react, and respond belowground to neighboring plants, and shed light on this highly plastic, ecologically significant, and dynamic process that remains almost entirely unnoticed.

# Author Contributions

This work was made possible through intensive, cross-discipline collaboration. Individual contributions are as follows: AP for the experimental design, plant production, data analyses, and lead authorship. JS for 3D metrics, image construction, and major intellectual contributions. JP for intellectual contributions, data processing, and developing novel root-tracking algorithms. TB for intellectual contributions, including experimental design and critical revisions.

# Acknowledgments

We acknowledge the invaluable assistance provided by Cornell's nano- and micro-CT facility run by Mark Riccio. We thank Professor Neil Matson for providing us with reservoirs and tools needed to grow our plants, professor Anthony Reeves for his assistance in developing preliminary root tracking algorithms, and Professor Joseph Fetcho for providing key software. We would also like to acknowledge Annika Kreye and Lindsey Saum for their help in designing and building hydroponic infrastructures, as well as Kim Goodwin and Kendra Hutchins for their help in the greenhouse. This work was partially supported by the Mario Einaudi Center for International Studies, and JS was funded through an NSF GRFP fellowship.

# Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015. 00274/abstract

Figure S1 | 3D reconstruction of all solitary and paired root systems (axes in 0.1 mm increments). (A) Control spruce (*Picea mariana*) (*n* = 3 containers, 3 individuals). (B) Intra-specific spruce (*n* = 3 containers, 6 individuals). (C) Inter-specific aspen (*Populus tremuloides*) and spruce root systems (*n* = 3 containers, 3 individuals each species). Aspen is highlighted in orange, spruce is highlighted in blue. (D) Intra-specific aspen (*n* = 1 container, 2 individuals). (E) Control aspen (*n* = 3 containers, 3 individuals). Note the differences in x- and y-axes.

# References


Medicago truncatula plant roots. Proc. Nat. Acad. Sci. 109, 16794–16799. doi: 10.1073/pnas.1209287109


using X-ray computed tomography. J. Exp. Bot. 61, 311–313. doi: 10.1093/jxb/ erp386

Volder, A., Smart, D. R., Bloom, A. J., and Eissenstat, D. M. (2005). Rapid decline in nitrate uptake and respiration with age in fine lateral roots of grape: implications for root efficiency and competitive effectiveness. New Phytol. 165, 493–502. doi: 10.1111/j.1469-8137.2004.01222.x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Paya, Silverberg, Padgett and Bauerle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Unraveling the hydrodynamics of split root water uptake experiments using CT scanned root architectures and three dimensional flow simulations

Nicolai Koebernick 1 † , Katrin Huber <sup>2</sup> \* † , Elien Kerkhofs <sup>3</sup> , Jan Vanderborght 2, 3 , Mathieu Javaux 2, 4, Harry Vereecken<sup>2</sup> and Doris Vetterlein<sup>1</sup>

<sup>1</sup> Department of Soil Physics, Helmholtz Centre for Environmental Research (UFZ), Halle, Germany, <sup>2</sup> Agrosphere (IBG-3), Forschungszentrum Jülich GmbH, Jülich, Germany, <sup>3</sup> Department of Earth and Environmental Sciences, KU Leuven, Leuven, Belgium, <sup>4</sup> Earth and Life Institute/Environmental Sciences, Université Catholique de Louvain, Louvain-la-Neuve, Belgium

Split root experiments have the potential to disentangle water transport in roots and soil, enabling the investigation of the water uptake pattern of a root system. Interpretation of the experimental data assumes that water flow between the split soil compartments does not occur. Another approach to investigate root water uptake is by numerical simulations combining soil and root water flow depending on the parameterization and description of the root system. Our aim is to demonstrate the synergisms that emerge from combining split root experiments with simulations. We show how growing root architectures derived from temporally repeated X-ray CT scanning can be implemented in numerical soil-plant models. Faba beans were grown with and without split layers and exposed to a single drought period during which plant and soil water status were measured. Root architectures were reconstructed from CT scans and used in the model R-SWMS (root-soil water movement and solute transport) to simulate water potentials in soil and roots in 3D as well as water uptake by growing roots in different depths. CT scans revealed that root development was considerably lower with split layers compared to without. This coincided with a reduction of transpiration, stomatal conductance and shoot growth. Simulated predawn water potentials were lower in the presence of split layers. Simulations showed that this was related to an increased resistance to vertical water flow in the soil by the split layers. Comparison between measured and simulated soil water potentials proved that the split layers were not perfectly isolating and that redistribution of water from the lower, wetter compartments to the drier upper compartments took place, thus water losses were not equal to the root water uptake from those compartments. Still, the layers increased the resistance to vertical flow which resulted in lower simulated collar water potentials that led to reduced stomatal conductance and growth.

Keywords: split-root, R-SWMS, root water uptake, plant root growth, Vicia faba

#### Edited by:

Pierre Dutilleul, McGill University, Canada

#### Reviewed by:

Lars Hendrik Wegner, Karlsruhe Institute of Technology, Germany Stanislaus Josef Schymanski, Swiss Federal Institute of Technology Zurich, Switzerland Pierre Dutilleul, McGill University, Canada Jonathan A. Lafond, Université Laval, Canada

#### \*Correspondence:

Katrin Huber, Agrosphere Institute (IBG-3), Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Straße, D - 52425 Juelich, Germany k.huber@fz-juelich.de

† These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science

> Received: 30 October 2014 Accepted: 09 May 2015 Published: 29 May 2015

#### Citation:

Koebernick N, Huber K, Kerkhofs E, Vanderborght J, Javaux M, Vereecken H and Vetterlein D (2015) Unraveling the hydrodynamics of split root water uptake experiments using CT scanned root architectures and three dimensional flow simulations. Front. Plant Sci. 6:370. doi: 10.3389/fpls.2015.00370

# Introduction

Water scarcity is an important abiotic limitation to plant growth and agricultural productivity. Under water limited conditions, changes in root system architecture (RSA) play a major role to reach locations where water is still present, which is often the subsoil. There is no simple relationship between the amount of roots present in certain locations and the actual root water uptake (RWU) from these sites (Pohlmeier et al., 2008). RWU is repeatedly described as a sink moving down the profile with time, only weakly related to root length density in a certain depth (Hainsworth and Aylmore, 1986; Pierret et al., 2003; Garrigues et al., 2006). In many of these studies change in soil water content in a certain depth is assumed to be synonymous with root water uptake. The illustrative Martini glass analogy first used by Zwieniecki et al. (2002) demonstrates that this assumption is too simple. When drinking a sip of Martini with a straw, the Martini is taken up from the bottom of the glass, but a change in "Martini content" is only observed in the upper layer of the glass due to the very high hydraulic conductivity within the glass. Roots and soil matrix are much more complex than the Martini-glass system; however, in soil-plant system the soil hydraulic conductivity and resulting soil hydraulic redistribution also obstruct the view on the site of root water uptake and its temporal dynamics. This has been known for a long time and a number of strategies have been developed to overcome this problem.

An experimental strategy to prevent soil hydraulic redistribution is to divide the root zone into different compartments, which prevent water flow between compartments to permit controlled heterogeneous distribution of soil moisture (Drew, 1975; Herkelrath et al., 1977). In case of horizontal splits, the split layers should additionally be penetrable by roots, which can be, for example, achieved by applying wax or paraffin. When roots take up water in a given compartment the change in total water content can be directly related to root water uptake from this compartment. This assumption can, however, only be drawn if the split layers are completely hydraulically isolated. In the case of water redistribution through the layers, the leakage rate has to be known. Another problem to determine RWU from a soil compartment arises due to the non-linearity of the soil water retention curve. Water content or soil water potential is usually measured at discrete points in the soil. When roots take up water from the soil, strong gradients in soil water potential can develop around the roots. Thus, an extrapolation between point measurements to the complete soil compartment becomes erroneous. A second experimental strategy is to directly observe water flux in soil as it has been successfully demonstrated by Zarebanadkouki et al. (2012). They imaged water flow into roots using neutron imaging of deuterated water. However, this method is hitherto either constrained to quasi two-dimensions (rhizotrons) or very small root systems and to short time scales.

An alternative approach is to quantify the amount of water being translocated by root or soil hydraulic redistribution. Mechanistic root water uptake models that describe water flow in soil, into, and within roots allow quantifying and locating root water uptake and redistribution of water within the soil and root system. The use of mechanistic models, like R-SWMS (root-soil water movement and solute transport, Javaux et al., 2008), has two prerequisites: (i) that the dominant processes are known and (ii) that the required input parameters are available. To fulfill the latter, dynamic information about RSA as well as hydraulic properties of individual root segments have to be available.

RSA has been obtained in the past using root growth models, i.e., RSA is artificially created based on a set of crop specific parameters and rules (e.g., branching rules, growth rates, etc.) derived from experiments (Clausnitzer and Hopmans, 1994; Lynch et al., 1997; Pagès et al., 2004; Leitner et al., 2010). Mostly, one or several typical realizations of RSA obtained from such models for a plant of a certain age have been used to calculate different scenarios, like root water uptake from saline soils (Schröder et al., 2013), performance of varying root architectural traits under different soil moisture regimes (Leitner et al., 2014), or the impact of stomatal regulation type on root water uptake (Huber et al., 2014).

Root growth models have been used as an alternative to 3Ddata of root systems as these were not available in the past. However, such data are now becoming increasingly accessible with non-invasive methods reaching a level of resolution which is sufficient to visualize most or all of the root system. The most advanced techniques for imaging soil-grown roots include X-ray computed tomography (Mooney et al., 2012), neutron radiography (Oswald et al., 2008), magnetic resonance imaging (Pohlmeier et al., 2008), or transparent soils (Downie et al., 2012). These techniques are of particular interest because they allow for repeated measurements. When ionizing radiation is used, it is however important to choose appropriate scan parameters to minimize potential damage to living tissues (Dutilleul et al., 2005; Zappala et al., 2013). Previous studies clearly demonstrated the potential of X-ray CT to analyze the temporal dynamics of growing roots (Jenneson et al., 1999; Gregory et al., 2003; Lontoc-Roy et al., 2005). While these early studies were limited to young seedlings, more recent work shows that the same is possible for considerably older root systems (Han et al., 2008; Tracy et al., 2012; Koebernick et al., 2014). First modeling approaches based on the use of RSA from non-invasive imaging are available (Stingaciu et al., 2013). The second challenge remains, i.e., the scarcity of data on root hydraulic properties. Measured data are primarily from hydroponically grown very young root systems. Certain assumptions have to be made to separate radial and axial conductivity during the measurements. Nevertheless, there is a wealth of information on how conductivity changes during root development and these have been used to scale the conductivity of individual root segments (Doussan et al., 1998, 2006). As roots age the resistance in the axial pathway typically decreases due to the maturation of xylem vessels, while in the radial pathway resistance increases with the development of apoplastic barriers (Frensch and Steudle, 1989; Bramley et al., 2009).

In order to avoid confounding root water uptake and hydraulic redistribution by the interpretation of local changes in soil water content we have chosen two of the above strategies: (i) an experimental approach of introducing barriers to avoid soil hydraulic redistribution; (ii) a modeling approach which takes soil and root hydraulic redistribution into account.

The objective of the current study is to compare experimental (introducing barriers to avoid soil hydraulic redistribution) and modeling approaches (calculation of soil and root water flow) with respect to their capacity to localize root water uptake in the presence of strong gradients in soil water potential. Local changes in soil water content will be compared to measured and modeled root water uptake.

For the experimental approach we combined a classical set up using wax barriers (Drew, 1975) with quantitative measurement of RSA over time via X-ray CT. This setup allowed the observation of the relation between RSA and water uptake and how it is affected by soil drying. The addition of paraffin layers allowed for the development of strong spatial heterogeneities in soil water potential, as is generally the case under field conditions.

For the modeling approach we used the mechanistic 3D model R-SWMS (Javaux et al., 2008), which enables a detailed description of soil and root water flow. While R-SWMS so far has only been applied for static (non-growing) root systems, mostly created by root architectural models, we now extended the existing model by an additional root development module, which uses the measured CT-data of RSA over time. Doussan's concept of changing axial and radial conductivity with age (Doussan et al., 2006) was included by using his root hydraulic parameterization by assigning these parameters to root age classes derived from the time lapse 3D RSA CT-Data.

Apart from modeling the actual experimental setup, root distributions obtained from split experiments were also used in simulations without splits and vice versa. This approach allowed us to (i) reinterpret measurement results, (ii) show the influence of split layers on plant water potentials that could be linked to differences in plant/root growth and eventually on root water uptake and (iii) show where soil water is taken up during root growth.

# Materials and Methods

#### Experiments

Two subsequent experiments under the same environmental conditions (growth chamber, 23◦C day/18◦C night, 65% relative humidity, photoperiod of 14 h, photon-flux density of 350µmol m−<sup>2</sup> s −1 ) were conducted with Vicia faba L. cv. Fuego.

The first experiment (3 replications), which will be referred to as "NoSplit" in the following, was conducted with homogeneously filled soil columns of 21.5 cm height with unrestricted soil water flow. The second (4 replications), referred to as "Split" was similar to the first one, but paraffin layers at 5, 10, and 15 cm height were established to interrupt soil water redistribution. This method was adopted from Drew (1975), who showed that root growth was unaffected by such layers. Both experiments were conducted consecutively, which explains the differences in the two setups.

#### Experimental Setup **"NoSplit" (without paraffin layers)**

The porous substrate was prepared by mixing quartz particles of different size classes, consisting of 85% sand, 10% silt, and 5% clay (Vetterlein et al., 2007). Additionally 50 g kg−<sup>1</sup> of gravel (2– 3 mm Ø) and 20 g kg−<sup>1</sup> of plastic beads (polypropylene, 2–3 mm Ø) were added to the substrate as internal reference for digital image analysis.

PVC cylinders (inner Ø = 12.5 cm, h = 21.5 cm) were filled up with the substrate by passing it through two sieves of 4 mm mesh size separated by a distance of 10 cm. This procedure was chosen to avoid particle size separation during filling. Resulting bulk density of the substrate was 1.52 ± 0.01 g cm−<sup>3</sup> . The cylinders had porous plates at the lower end (**Figure 1A**), which were connected with plastic tubing to a water source. The soil was gently watered with a nutrient solution (modified from Römheld and Marschner, 1990) by capillary rise from the bottom of the sample (soil water potential ψ = 0 hPa at z = −21.5 cm). Average volumetric soil water content (θ) at the start of the experiment was 31.1 ± 1%. Vicia faba seeds were surface sterilized in 10% H2O<sup>2</sup> solution for 10 min, thoroughly rinsed in deionised water and subsequently imbibed for 1 h in a saturated CaSO4solution. Seeds were placed on wet blotting paper and placed in a dark cabinet at room temperature for 2 days. For each cylinder, one pre-germinated seed was carefully placed in a prepared cavity in the soil at a depth of 1 cm. The soil surface was covered by a 2 cm layer of fine quartz gravel. Until shoot emergence columns were covered with aluminum foil to further minimize evaporation. With the removal of aluminum foil the drying period was initiated (Day 6).

### **"Split" (with paraffin layers)**

The substrate was the same as in the "NoSplit" experiment, however, without the addition of plastic beads as these caused problems in the segmentation procedure (see below). Soil bulk density was slightly higher (1 = 0.12 g/cm<sup>3</sup> ).

For the split layers, molten paraffin was casted and flattened to a thickness of approximately 0.5 mm and cut into a circular shape. At -5, -10, and -15 cm depth a layer of paraffin was placed on top of the soil and sealed to the cylinder walls using molten paraffin (**Figure 1B**). For initial irrigation, we placed rhizon soil moisture samplers (Eikelkamp, Giesbeek, NL) in each soil compartment. Those were connected over night to bottles filled with 150 ml nutrient solution each. Volumetric water content at the start of the experiment was 23.8 ± 0.5% in each compartment. Seed preparation was the same as in the "NoSplit" experiment. To avoid the formation of cracks in the soil due to the placement of large Vicia faba seeds, these were planted in a separate seed compartment: a cylinder (Ø = 6 cm, h = 3 cm) filled with the soil mixture and 20 ml of water. When the roots emerged through the paraffin layer at the bottom of the seed compartment, the small cylinder was placed on the topsoil (Day 0). The remaining bare topsoil was covered with gravel to reduce evaporation. The split samples were initially also covered with aluminum foil, which was removed on Day 4 to start the drying period.

#### Transpiration and Soil Matric Potential

The PVC cylinders were placed on weighing cells (KERN 572, Kern and Sohn GmbH, Balingen, Germany), and grown for 30- 36 days with no additional watering. Weight data were recorded

are given in cm.

every 10 min throughout the experimental period. Four microtensiometers (Vetterlein et al., 1993) were inserted horizontally through sealed boreholes ("NoSplit": -1.5, -6.5, -11.5, and - 16.5 cm soil depth; "Split": -2.5, -7.5, -12.5, -17.5 cm, **Figure 1**) to monitor the soil matric potential (ψm), during drying.

The daily transpiration rate was calculated from weight differences between two subsequent days. Evaporation was assumed to be negligible due to the layer of coarse gravel on the surface and as surface was never rewetted during the experiment. Relative humidity was constant day and night hence dew formation could also be excluded. Only on the seed compartment used in "Split" experiment, there was no gravel layer and hence water applied initially (20 ml) was assumed to be lost by evaporation uniformly within the first 7 days.

Leaf area development was estimated by daily measuring the length and width of the lamina of each leaflet and using the linear model of Peksen (2007):

$$LA = 0.919 + 0.682 \, L \, \* \, W \tag{1}$$

where LA [cm<sup>2</sup> ] is the one-sided leaf area, L [cm] is the length of the lamina, and W [cm] is the width of the lamina. After harvest, we used a flatbed scanner to measure leaf area. The results agreed well with the estimation using Peksen's model. Stomatal conductance was measured at the end of each day using a steadystate porometer (SC-1 Leaf Porometer, Decagon Devices, Inc., Pullman, WA, USA). Two measurements per plant were taken on the abaxial side of the youngest unfolded leaf pair and the mean value of the two measurements was stored.

#### CT Scanning and Image Analysis

All samples from the "NoSplit" and the "Split" experiment were scanned every second day during the night phase with an industrial X-ray micro-CT scanner (X-Tek HMX 225) with a fine focus X-ray tube. The scanning parameters are summarized in **Table 1**. Potential X-ray dose was estimated using the free

TABLE 1 | Table 1 X-ray settings used in the different experimental setups.


online tool Rad Pro Dose Calculator (McGinnis 2002-2009). In the "Split" experiment, which had a higher exposure, cumulative dose at the end of the experiment was 4.8 Gy. This is well below the maximum dose (approximately 30 Gy) suggested for plant CT studies by Zappala et al. (2013). Due to the height of the cylinders separate scans of the upper and the lower part of the sample had to be performed. In the NoSplit setup the mechanism for attaching the porous plate to the soil cylinder at the bottom required an additional plastic ring for sealing reasons which caused photon starvation at the lower end (7 cm), so that not the entire root system could be imaged.

Although the samples were positioned carefully, images scanned at different times were not perfectly aligned. A manual, feature-based method was used to register the images (see Koebernick et al., 2014). The scans from the upper and lower halves of the samples were combined into a single image. The raw images were filtered with a total variation filter (Rudin et al., 1992) to remove small scale noise while preserving sharp edges. We additionally used a pseudomedian filter (Pratt, 1991) to enhance the contrast between roots and soil and to remove beam hardening artifacts. Roots were segmented from the background using a region growing algorithm, similar to the approach of Kaestner et al. (2006). The algorithm used two thresholds to determine, whether a voxel belongs to the root system. The thresholds were chosen manually based on the histogram and visual inspection of the segmentation results. The images were processed with the freely available software QtQuantim (www.quantim.ufz.de). A more detailed description of the technical procedure can be found in Koebernick et al. (2014). In the NoSplit experiment, two samples (NoSplit 1 and NoSplit 3) could not be successfully segmented due to technical difficulties. Due to improved scanning conditions for the Split setup all architectures could be segmented. The segmented images of the root systems are shown in **Figure 2A**. These images contained a number of misclassified voxels (e.g., wall material, paraffin layers, cracks, tensiometers) and roots were disconnected at some points.

For the subsequent simulations, a connected root structure was required. Thus, the binary images had to be manually reconstructed using a three-dimensional virtual reality system, which was initially developed to reconstruct MRI data but can be used for any binarized images (for a detailed description of this method see Stingaciu et al., 2013). Due to the laborintensive manual reconstruction only two replications of the "Split" (Split 1 and Split 3) experiment were reconstructed. We chose Split 1 and Split 3 because these cover the contrasting root architectures in the "Split" experiment. Misclassified regions in the binarized CT images could be excluded by this manual procedure.

For the determination of root age of each segment at each time step, the reconstructed and stored root system of the precedent scan was opened simultaneously with the image of the subsequent scan. Using the overlay of both scans newly grown roots could be identified and added to the existing root structure. The temporal resolution of the growing root architecture was limited by the time interval between two CT scans (2 days). To obtain smoother root growth, the origination time t<sup>s</sup> of a segment s that grew between times t<sup>i</sup> and ti+<sup>1</sup> when a CT scan was made, was calculated using Equation 2:

$$\mathbf{t}\_s = \mathbf{t}\_i + \frac{l\_s}{\Delta l\_s} \left(\mathbf{t}\_{i+1} - \mathbf{t}\_i\right) \tag{2}$$

where 1l<sup>s</sup> [L] is the length of all segments that grew between time t<sup>i</sup> [T] and ti+<sup>1</sup> and that are connected to the same connection point of the root system at time t<sup>i</sup> as the root segment s, and ls is the length of all segments that are closer to the connection point than segment s and therefore should have emerged before segment s. The average length of a manually reconstructed root segment was 0.087 ± 0.008 cm.

#### Destructive Measurements

At the end of the experiment (Day 31–35) roots were extracted from the soil by washing using sieves of 3 and 2 mm mesh size successively. In the "Split" experiment, compartments were analyzed separately. In the "NoSplit" experiment, the roots grown into the lower 7 cm of the cylinder that could not be imaged were harvested separately. Roots were stored in Rotisol and subsequently scanned on a flatbed scanner (EPSON Perfection V700 PHOTO). The images were analyzed with WinRHIZO 2009b (Regent Instruments, Inc., Quebec, Canada) to obtain total root lengths.

#### Modeling of RWU

For the simulation of RWU we used the numerical model R-SWMS, which solves the water flow equation in the root network and in the soil (Javaux et al., 2008). The numerical solution of the Richards equation (Equation 3, Richards, 1931) with a sink term S based on SWMS\_3D (Simunek et al., 1995).

The water flow equation for the root network is solved based on the radial and axial flow equations (**Equations 4 and 5**) and the mass balance at each root node, resulting in a system of linear equations for ψx, the xylem water potential (Doussan et al., 1998). The system is solved with a biconjugated gradient method.

The root and the soil water flow equations are coupled through the definition of the sink term of the Richards equation and of the water potential at the soil-root interface for the Doussan equation. The sink term of the Richards equation is defined as the sum of the radial root flow into all root segments, k, located within a soil voxel (cuboid), i, divided by the cuboid volume (Equation 6). The soil-root interface water potential at each root node is defined as the distance weighted average of the water potential at the soil voxel nodes.

$$\frac{\partial \theta}{\partial t} = \nabla \cdot \left[ K\left(\psi\right) \nabla\left(\psi\right) \right] + \frac{\partial K\left(\psi\right)}{\partial z} + S(x, y, z, t) \tag{3}$$

$$J\_r = K\_r^\* A\_r \left(\psi\_{s,int} - \psi\_x\right) \tag{4}$$

$$J\_{\chi} = -K\_{\chi}^{\*} A\_{\chi} \left( \frac{d\psi\_{\chi}}{dl} + \frac{dz}{dl} \right) \tag{5}$$

$$S\_i = \frac{\sum\_{k=1}^{n\_k} J\_r^k}{V\_j} \tag{6}$$

where θ [L<sup>3</sup> L −3 ] is the volumetric water content of the soil, K [L T−<sup>1</sup> ] the soil hydraulic conductivity, ψ [P] the soil matric potential, and z [P] the gravitational potential. S [L<sup>3</sup> T −1 ] is the sink term, J<sup>r</sup> [L<sup>3</sup> T −1 ] the radial flow into the roots, J<sup>x</sup> [L<sup>3</sup> T −1 ] the axial flow in the root xylem, K ∗r [L T−<sup>1</sup> P −1 ] is the radial conductivity, K ∗ x [L<sup>2</sup> T <sup>−</sup><sup>1</sup> P −1 ] the axial conductivity, ψs,int [P] is the water potential at the root-soil interface and ψ<sup>x</sup> [P] the xylem water potential, A<sup>r</sup> and A<sup>x</sup> [L<sup>2</sup> ] are the lateral surface and the cross sectional areas of a root segment, l [L] is the length of a root segment. The axial conductance, K<sup>x</sup> = K ∗ <sup>x</sup>A<sup>x</sup> [L<sup>4</sup> T <sup>−</sup><sup>1</sup> P −1 ]. The indices i and k stand for discrete soil voxels and root segments, respectively. V<sup>j</sup> [L<sup>3</sup> ] is the volume of a single soil voxel.

The equivalent hydraulic conductivity of the root system, Kroot [L<sup>3</sup> P <sup>−</sup><sup>1</sup> T −1 ], is defined by the relation between actual transpiration, Tact [L<sup>3</sup> T −1 ] and the difference between the effective soil water potential and the root collar potential (Javaux et al., 2013).

$$T\_{act} = K\_{root} \left(\psi\_{s,eff} - \psi\_{collar}\right) \tag{7}$$

$$
\psi\_{s,eff} = \sum\_{j} \text{SUF}\_{j} \,\psi\_{s,int} \tag{8}
$$

FIGURE 2 | (A) Three dimensional rendered view of the segmented CT images at different scan times. White arrows indicate misclassified objects: NoSplit 2, Day 8: plastic bead Split 1, Day 10: tensiometer, Day 22: paraffin layer, Day 30 soil crack. Split 3, Day 34: container wall. White boxes at Day 8 or 6 show the scaling of the root system: the distance between two ticks

equals 100 pixels, which equals 2.45 cm for NoSplit2 and 2.77 cm for the Split setups. (B) VR reconstructions of root system architectures at the end of each experiment within their respective soil Root systems are colored according to root age and the soil according to the simulated soil water potential.

where ψs,eff [P] is the effective soil water potential, which is weighted by the standard uptake fraction, SUF<sup>j</sup> [-]. SUF<sup>j</sup> represents the relative water uptake by a root segment j in a soil profile with a uniform soil water water potential and can be derived by solving the Doussan equations. A more detailed explanation can be found in Couvreur et al. (2012).

The R-SWMS code and a manual as well as the reconstructed root architectural files are available upon request from the authors.

#### Model Setup

The samples NoSplit 2 from "NoSplit" experiment and Split 1 and Split 3 from "Split" experiment, with fully reconstructed root architectures, were used for the setup of virtual experiments in R-SWMS. In the following when referring to modeling data names of samples will be written in italics.

#### **Soil domain**

We defined rectangular domains with a discretization of 0.5 × 0.5 ×0.25 cm<sup>3</sup> . The domain size was 14 × 14 × 21.5 cm<sup>3</sup> for the "NoSplit" experiment. The domains of the "Split" experiment differed in the z-direction (z = 20 cm for Split 1; z = 20.25 cm for Split 3, **Figure 2B**). The cylindrical geometry of the soil columns was approximated using Pythagoras' Theorem with a cylinder radius of 7 cm. Voxels belonging to this cylinder were defined as soil material; voxels on the outside were defined as wall material. The water retention characteristic was described by a bimodal Mualem - van Genuchten expression (Van Genuchten, 1980; Durner, 1994). The soil hydraulic parameters in **Table 2** were derived from separate HyProp measurements (Peters and Durner, 2008), except the saturated hydraulic conductivity, K<sup>s</sup> , which was predicted using the Rosetta tool (Schaap et al., 2001). Paraffin layers were defined as 0.5 cm thick layers within the cylinder. The modeled layer thickness is thus 10 times larger than the thickness of the split layer in the experiment. However, to achieve a reasonable simulation speed, we had to settle for this trade-off. The split layer material was defined equal to the wall material. However, as a certain leakiness of the split layers became obvious during the time course of the experiment and later on during the modeling, we decided to simulate the leakage by assigning a small hydraulic conductivity to the layers of concern. All soil boundary conditions were defined as zero flux. Initial conditions were defined according to the initial water content at the start of the drying period in the experiments. In the "NoSplit" setup soil matric potential was at hydrostatic equilibrium and in the Split setup, soil water content was equal in each compartment.

#### **Root architecture**

The root architectures for the simulations were obtained from the manually reconstructed CT images. Root hydraulic properties were based on an age dependent parameter set by Doussan et al. (2006) for Lupinus angustifolius (**Figure 3**, bold lines). Radial conductivity of roots was given a constant value of 8.64 × 10−<sup>4</sup> cm d−1hPa−<sup>1</sup> . The axial conductances increased stepwise with segment age. In Doussan et al. (2006) axial conductance (i.e., xylem conductance) of lateral roots increased with age, whereas taproot axial conductance increased with distance to the tip. Thus, for the taproot we had to convert our age information to distance information. For this we divided the given distances by the mean measured elongation rate of the taproot (0.7 cm d−<sup>1</sup> ) to translate the given distances to the according ages.

At a given simulation time only the root segments with an origination time smaller than the actual simulation time were taken into account. The root system was updated at each further run-time step thus enabling predefined root growth over time. We converted the measured daily transpiration rates of each sample to a periodic step function with zero flow during the night and so defined the root flow boundary conditions in the model at the root collar.

#### Scenarios

Each of the three samples was exposed to two or three scenarios to analyze the effect of paraffin layers on RWU. In the first scenario



Saturated and residual water content, θ <sup>s</sup>and θr, respectively; van Genuchten shape parameters, α and n; pore connectivity parameter λ; and saturated hydraulic conductivity, Ks. For the soil, a bimodal θ(ψ) relation (Durner, 1994) was used. Asterisk indicates the saturated hydraulic conductivity of paraffin for the scenario SC.

(CD), a continuous soil domain without any split layers was used. In the second scenario (NC), we defined three non-conductive paraffin layers. Finally, the third scenario (SC), aimed to achieve best agreement to measured data for the "Split" experiment by considering leaking paraffin layers and assigning a low hydraulic conductivity of 0.001 cm d−<sup>1</sup> (**Table 2**) to the split layers. Sample Split 1 was simulated with three slightly conductive layers, and Split 3 with a non-conductive layer at −5 cm and two remaining slightly conductive layers.

A sensitivity analysis was performed to evaluate the uncertainties in the modeling approach due to uncertain age dependent root hydraulic conductivities. We focus on predawn water potentials, ψd, since simulated soil water potentials could be compared with measurements and transpiration rates were used as boundary conditions. Equation 7 shows that in case of zero transpiration, e.g., during night, ψs,eff = ψcollar. Thus, predawn water potential is independent of Kroot and SUF can be used as an indicator for the impact of different root hydraulic conductivities on ψpd. Since SUF represents the water uptake by a root segment, relative to the total of the uptake of the root system, SUF does not depend on the absolute (radial and axial) conductivities of the root segment but on the ratios between the conductivities of one segment to other segments.

The variability of SUF induced by different age dependencies of the hydraulic parameters was examined by comparing different combinations of age dependent and constant axial and radial conductivities for the different reconstructed root architectures (NoSplit2, Split1, Split3) at the end of the growth period. The constant value for K<sup>x</sup> was defined as the arithmetic mean of the age dependent K<sup>x</sup> values and age-dependent K ∗ r values were modified from Doussan et al. (1998) who defined age-dependent K ∗r values for Zea mays L. (**Figure 3**). An overview of the parameterization is given in **Table 3**.

# Results

#### Experimental Results

As expected, plant performance differed markedly between the two experiments (**Figure 4**). In the "NoSplit" experiment plants were bigger and had a larger leaf area (**Figure 4A**). Leaf growth was initially the same in both experiments, but after Day 15 leaf area increased more in the "NoSplit" experiment. A similar pattern could be observed for total root lengths obtained from CT images over time (**Figure 4C**). Root elongation was similar for both, "Split" and "NoSplit" experiment until Day 10. Afterwards elongation rate was higher for "NoSplit." Root length estimations

TABLE 3 | Perturbations of root hydraulic conductivities from Figure 3 for the sensitivity analysis.


from destructively harvested roots using WinRHIZO were on average higher than estimations from CT (**Table 4**).

The vertical root length distribution in the "Split" experiment differed between Split 1 and the remaining samples. Compartment I in Split 1 contained about 3/4 of the total root length, while the distribution for the other replications of the "Split" experiment was more even (**Table 4**). In the "NoSplit" experiment root density increased with depth.

In both experiments transpiration rate initially increased with leaf area (**Figure 4B**). In "NoSplit" a sharp decrease in transpiration rate was seen at Days 23, 25, and 28, respectively for the different samples. Transpiration reduction occurred earliest in NoSplit 3, which was also the largest plant with the highest transpiration rate up to that day. In the "Split" experiment, transpiration reduction could be observed earlier, although the reduction in transpiration was not as strong as in the "NoSplit" experiments. The lower leaf areas and smaller transpiration rates in the "Split" experiment were accompanied by lower stomatal conductance of the youngest unfolded leaves in comparison to the "NoSplit" experiments (**Figure 4D**). Stomatal conductance decreased already from the first measurement, i.e., Day 10, in the "Split" experiment. In the "NoSplit" experiment the variability of stomatal conductance in the different samples was very high, but low values were not measured until Days 23 or 24, respectively.

The addition of paraffin layers ("Split" experiment) also had a pronounced effect on the temporal development of the soil matric potentials in the different soil compartments (**Figures 5A–C**). For the sake of brevity we only present the results of the samples that were later used for modeling (the remaining samples behaved similarly, see **Supplementary Figure 1**). In NoSplit 2, soil matric potential remained high during a long period (approximately until 25 days after the start of the experiment) and there were only small differences between the matric potentials at different depths. After 25 days, the time at which the transpiration in the no-split experiment started to decrease (**Figure 4B**), the matric potentials decreased strongly and more or less simultaneously at different depths in the column. For the "Split" experiments, the matric potentials started to decrease much earlier (from Day 10 onwards) and sequentially from the top toward the bottom compartments. Except for the upper compartment in Split 3, the decrease of matric potential was more gradual and less abrupt than in the "NoSplit" experiments. The tensiometer readings for the "Split" experiment showed a pronounced day-night cycle in the upper and a more damped diurnal signal in the lower compartments.

Water depletion from each compartment was calculated from measured tensiometer values assuming a uniform matric potential within a layer and using the substrate specific water retention curve (**Table 2**). These data were compared to total water loss derived from weighing cells (**Figure 6**). When air bubbles started to form in the tensiometers no further water content change could be calculated. The calculated water content at this point was between 9.5 and 10.6% (ψ<sup>m</sup> = −745 to −431 hPa). In the "NoSplit" setup (**Figure 6A**) there were no true compartments, we therefore assumed that the tensiometers represented the matric potential for the surrounding volume closest to the tensiometer. While the difference between

calculated and measured cumulative water depletion for the "Split" setup (**Figures 6B,C**) converged to below 10% (+9% Split 3, −5% Split 1) at the end of the experiment, it was much higher (17%) in the "NoSplit" setup. Comparison of the slopes over time indicates a poor fit of the dynamics. Calculated water depletion was clearly overestimated at the beginning and underestimated toward the end of the experiments, especially in Split 3.

The arrival of roots in Compartments III and IV in Split 1 was at Day 12 and 18, respectively, nonetheless there was significant (even if overestimated) water depletion from both compartments before these dates.

#### Simulation Results

The three samples (NoSplit 2, Split 1, and Split 3) representing different RSA were subjected to three different scenarios: (CD), a continuous, unrestricted soil domain, (NC) a soil domain with non-conductive split layers, and (SC) with semi-conductive split layers. Mean simulated soil matric potentials in four layers were compared to the measured tensiometer values (**Figure 5**).

#### Choice of Scenario

In scenario (CD) (continuous soil domain) (**Figures 5D–F**), the simulated matric potentials in the different soil layers started declining strongly and nearly simultaneously only toward the end of the simulation period. The simulated decline occurred the earliest and was the strongest in the "NoSplit" experiment reflecting the larger cumulative transpiration from this experiment.

For the "NoSplit" experiment, the simulated matric potentials for scenario (CD) showed a similar behavior as the measurements (**Figure 5D**). The timing and the slope of decrease fitted the experimental data well. The lowest tensiometer (−16.5 cm) was an exception, probably due to the fact that the deep roots could not be detected in the CT and were missing in the model.

For both samples of the "Split" experiment (**Figures 5E,F**), the measured matric potentials of the upper two tensiometers started decreasing much earlier than the simulated matric potentials for scenario (CD). This illustrates the effect of the paraffin layers on the soil water distribution in the "Split" experiment which is ignored in scenario (CD).

Scenario (NC) with non-conductive paraffin layers was simulated only for the "Split" experiments (**Figures 5G,H**). The simulated matric potentials at the tensiometer depths decreased sequentially from top to bottom and the time lag between these decreases was much larger than in scenario (CD) for the same samples. The simulated water potentials started to decrease shortly after roots arrived in a compartment. In Split 3 (**Figure 5H**), simulated average water potential in Compartment I decreased to about −2000 hPa until Day 15 and remained at this level thereafter only showing pronounced diurnal fluctuations until the end of the simulation run. In both samples of the "Split" experiment (**Figures 5G,H**) for scenario (NC) the simulated changes in water potential in Compartment IV were very small due to the small fraction of roots in this compartment.

With Scenario (NC) we were not able to reproduce the measured dynamics of soil matric potentials of the "Split"


TABLE 4 | Root length estimations from CT images and from destructive measurements at the end of each experiment.

\*Value for Compartments II and III combined.

samples. Measured matric potentials did not show a sequential stepwise decrease but a more gradual decrease that started earlier than the simulated decrease and sometimes even earlier than the root arrival time in a compartment. One exception was the matric potential in Compartment I of the Split 3 sample. Scenario (NC) produced large water potential differences between the different compartments, which were not in agreement with the measurements.

The previously described results indicate that paraffin layers were not perfectly isolating, but that there must have been water redistribution between neighboring compartments, albeit at a lower rate than in completely unrestricted soil. Thus, scenario (SC) was applied.

For Sample Split 1 in scenario (SC) (**Figure 5J**), the simulated matric potentials of Compartment I showed a slower decrease than those obtained with scenario (NC) or(CD). At the same time scenario (SC) resulted in an earlier decrease of matric potential in the lowest compartment compared to scenario (NC). The pronounced measured diurnal pattern of soil matric potential in Compartment I was successfully reproduced in scenario (SC).

Likewise, for Sample Split 3 simulated matric potentials of scenario (SC) showed the best agreement with measured tensiometer data. Here the assumption that all layers except the top layer were leaking was important for obtaining the good agreement.

As expected, for the "NoSplit" experiment (**Figure 5I**), agreement between measured soil matric potentials and those simulated with scenario (SC) was very poor. However, it is interesting to note the influence of, albeit leaking, hydraulic barriers to soil water potentials.

In contrast to experimental approaches, which can only detect changes in soil matric potential, the simulation results allow disentangling the different fluxes which contribute to local changes in matric potential and soil water content. The evaluation of fluxes was restricted to those simulations which showed the best agreement between measured matric potentials and simulated once, i.e., scenario (CD) for sample NoSplit 2, scenario (SC) for samples Split 1 and Split 3.

#### Simulated Flow Dynamics

The water balances of the single soil compartments are depicted in **Figure 7**. In case of impermeable split layers, the storage change within one soil compartment should equal root water uptake. However, if the split layers are leaking, which is the case for most of the layers, only adding the net flow through the split layers to the storage change equals root water uptake.

For the NoSplit 2 (**Figure 7A**) simulation RWU was largest in the upper compartment, where it started to decrease from Day 25 onward. The 5–10 cm layer only started to significantly contribute to RWU from Day 17 onward and the 10–20 cm layer only after Day 20, which is related to root arrival time.

It is interesting to note that "early morning values" of RWU in the 0–5 cm layer remained higher than those in the other layers even after 25 days, i.e., during a period where overall contribution of the lower layers to RWU had increased and total transpiration rate was reduced in the experiment.

Simulations showed soil hydraulic redistribution of water from the lower layers to the top 0–5 cm. At 5–10 cm depth inflows from the deepest soil layer and outflows to the 0–5 cm layer were almost of the same magnitude, so the resulting net flow oscillated around zero. Soil hydraulic redistribution started to decrease after Day 25 and seized after Day 31.

Since RWU from a layer corresponds to the sum of the net water flow into and the decrease of the water storage in a soil compartment, it is evident that RWU in a soil layer cannot be derived from water storage changes in that layer. RWU in the 0– 5 cm layer is considerably larger than the changes in water storage whereas the opposite is true for the 10–15 cm layer. It is clearly visible that RWU and storage change did not correspond to each other as long as there was significant soil hydraulic redistribution.

Substantial soil hydraulic redistribution occurred also in the samples Split 1 (SC) and Split 3 (SC), although K<sup>s</sup> values of paraffin layers were only 0.001 cm d−<sup>1</sup> (**Figures 7B,C**). In both simulations RWU did not correspond to water storage change with the exception of Compartment I in Split 3, which was assumed to be separated by a non-conductive split layer. RWU from Compartment IV was very small in both Split 1 (SC) and Split 3 (SC) while the change in soil water content was substantially higher due to flow across the split layer. The same pattern was observed in Compartment III, but net outflow of

water started earlier and was eventually compensated by inflow from Compartment IV. Compartment II showed a contrasting behavior between the two samples of the "Split" experiment. In Split 3 the non-conductive layer at the top prevented water movement in the soil to Compartment I, and the fraction of RWU from compartment II was considerably higher in Split 3 than in Split 1.

In both simulations of the "Split" experiment, there was significant hydraulic redistribution via deep roots into Compartment I. Root hydraulic redistribution was much more pronounced in Split 3. According to the simulations the redistribution via the roots occurred during night and the water was taken up by the roots during the next day.

The comparison of cumulative root water uptake from the different compartments with cumulative water depletion at the end of the simulations highlights the importance of including soil hydraulic redistribution when analyzing the pattern of RWU (**Table 5**). This is most obvious in the unrestricted sample NoSplit 2, where 69% of RWU occurred in the 0–5 cm layer, while the water depletion in this layer was only 16% of total water depletion. But even in Compartment I of Split 3, which was assumed to be perfectly isolated, RWU and water depletion are slightly different, which is probably due to the discretization of simulation outputs and rounding errors.

Further, the development of the root system architecture (**Figure 2**) can be compared to the water flows within the soil and root system (**Figure 7**). Due to the semipermeable split layers in Split 1, most of the RWU takes place in the upmost compartment, the location where also most of the roots are found. In Split 3, where the top compartment is hydraulically isolated, the roots take up most of the water from this layer within the first 15 days, while afterwards the uptake shifts to the lower compartments. This pattern is reflected in the RSA development. The NoSplit setup shows a more or less smooth shift of roots as well as RWU downward in the domain.

#### Sensitivity Analysis

Following Equation (8), the effective soil water potential, in case that transpiration is zero, is equal to the water potential at the root soil interface weighed by the standard uptake fraction, SUF. The SUF was calculated for four different parameterizations of root hydraulic conductivity. **Figure 8A** shows the sum of SUF for the NoSplit setup within given soil depth increments. With agedependent radial conductivity the SUF becomes more uniform over depth. For both Split setups the variability with the different parameterizations is not as large (see **Supplementary Figure 2**).

The SUF, which shows the hydraulic architecture of the root systems, are compared for the three different plants (**Figure 8B**). In contrast to the root system architecture, only small differences can be observed. The differences in predawn water potentials between the different plants were thus mainly due to the soil water distribution and less to RSA.

#### Pre-dawn Water Potential at the Root Collar

Simulated pre-dawn water potential at the root collar (ψpd) was used as an indicator for plant water status (**Figure 9**). ψpd is

time compared to cumulative transpiration from Day 8 for NoSplit (A) and Day 11 for Split 1 (B) and Split 3 (C) until the end of the experiment. Filled areas represent cumulative water content change in the different compartments calculated from tensiometer measurements. Gray line and circles represent cumulative transpiration measured with balances. White asterisks denote the point, when the tensiometer in the compartment showed air bubbles.

independent of actual transpiration rates and can therefore be used to compare different samples. ψpd is generally thought to be in equilibrium with the soil water potential provided that night induced interruption of transpiration is long enough and flow rates in soil root systems are high enough to reach this equilibrium (Donovan et al., 2003). However, the soil matric potentials, simulated in this study were clearly not in equilibrium, especially for the two split samples.

In sample NoSplit 2 (CD), simulated predawn ψpd decreased only slowly until Day 25 and was in equilibrium with soil matric potential in the topsoil (−1.5 cm depth). Due to the homogeneous soil water distribution it was also closely related to the matric potential in the wettest soil accessible to the plant, i.e., the soil at maximum rooting depth at each time step. From Day 25 onwards there was a strong decrease of soil matric potential in the whole column and an according decrease of ψpd. After Day 30, ψpd was more negative than the topsoil matric potential. The disequilibrium increased until the end of the experiment. In both split samples ψpd was more negative than the matric potential at maximum rooting depth but less negative than the topsoil matric potential, indicating that the system did not reach equilibrium at the end of the night. ψpd in Split 1 (SC) was closer to the matric potential in the topsoil, reflecting the higher redistribution through the split layers in Split 1 (SC).

To illustrate the impact of the split layers on soil and thus plant water status, predawn soil water potentials of the different scenarios with and without paraffin layers (SC vs. CD) for each sample were compared. The difference of absolute soil water potentials for the two contrasting soil environments was calculated (1/ψpd/ = /ψpd/SC − /ψpd/CD) (**Figure 10**, bold lines). As expected, soil water potential was constantly more negative in scenario SC than in scenario NC. 1ψpd in Split 1 and in NoSplit 2 were of the same magnitude, while in Split 3, where the upper paraffin layer was assumed to be non-conductive, it increased more rapidly and stronger, indicating an effect of the higher degree of hydraulic isolation of the different soil layers.

When using the previously calculated SUF to determine the impact of parameterization of root hydraulic conductivities on effective soil water potentials, the variability of soil water potentials compared to the plant variability is very small (**Figure 10**, thin lines).

#### Discussion

#### Influence of Paraffin Layers on Plant Growth

CT measurements gave insight into the changes of growth behavior caused by the addition of wax layers. However, the causes for these changes are not revealed by the CT measurements. By using a simulation model CT measured RSA and the low (zero) hydraulic conductivity of the wax layer could be linked to internal plant water potentials. This enables an interpretation of plant water stress and its implications for shoot and root growth.

Experimental results as well as simulations suggested strongly that most of the paraffin layers were not perfectly hydraulically isolating. Tomographic images and visual inspection after destructive harvest showed, however, no evidence of cracks or holes in the wax layers. It is possible that there were cracks at the container walls that were formed due to shearing of the paraffin caused by the weight of the soil in the upper compartments. The only paraffin layer that was evidently tight was consequently the uppermost layer in the sample Split 3. Drew (1975) suggested the use of layers as thin as 0.2 mm, which is even thinner than the layers that were used in this study. Another possible source of leakage is linked to diurnal shrinking and swelling of roots (Huck et al., 1970), which could lead to cavities in the paraffin where it is penetrated by roots. This could not be excluded as CT images were recorded during night.

Roots easily penetrated the paraffin and grew into the lower compartments. Taproots and vertically oriented laterals were not affected by paraffin layer. However, a few roots

continued to grow horizontally within the soft paraffin layers (see **Supplementary Figure 3**).

The plants in the "Split" experiment were overall smaller with lower root densities. Inserting split layers generated a substantial resistance to vertical water flow within the soil and hence water redistribution in the soil column. A restriction of this redistribution led to lower simulated predawn root and collar water potentials, which were related to lower measured stomatal conductance. The lower predawn water potentials pointed at plant stress that resulted in a restriction of root and shoot growth. Even though the root-shoot ratio was shown to increase in Vicia faba in drier environments (El Nadi et al., 1969), this could not be observed in this experiment. A possible explanation for this is the higher bulk density in the split experiment. Slight increases of soil strength can lead to a substantial reduction of root penetration rate (Taylor et al., 1966). We cannot exclude a possible effect of oxygen depletion on plant performance caused by the addition of paraffin layers, as no oxygen concentrations were measured. However, we feel that hypoxia is highly unlikely: The soil was initially not water saturated and the fact that paraffin layers were permeable to water means that soil air could move as well. The rhizon samplers were kept in the soil during the experiment as possible pathways for air. The soil mixture was an artificial mixture without added organic matter, so microbial respiration should be minimal. Experiments with the same quartz substrate showed that even close to saturation redox potentials only decreased after adding significant amounts of organic material (Ackermann et al., 2008).

#### Relation between Measured Water Loss and RWU

The simulations showed the discrepancy between change in soil water content and the location of root water uptake for individual soil compartments, which was caused by soil hydraulic



redistribution. Even a small conductivity of the hydraulic barriers led to considerable redistribution of soil water. The direct calculation of soil water content, and in extension RWU, from measured soil matric potentials was further complicated by the non-linear relation between water potential and water content, which precludes the extrapolation of a single tensiometer reading to the total soil compartment without knowing the gradients. The development of gradients around active roots is shown in **Supplementary Figure 4**.

Even if the vertical soil flow is completely restricted, hydraulic redistribution through the roots might still be a substantial amount of water that is exchanged between the roots and the soil in the drier regions of the root zone. In this case, however, the net water content change should correspond to net root water flow. The share of root hydraulic redistribution was higher when soil water redistribution was restricted by barriers, allowing the formation of a sufficing water potential gradient to drive flow. This may in part explain the controversy in literature as to the ecological relevance of root hydraulic redistribution. Its magnitude spans almost two orders of magnitude and is affected by numerous factors, such as root and water distribution, soil texture, and root-soil hydraulic conductance (Neumann and Cardon, 2012).

#### Predawn Collar Potential

Simulation results suggest that predawn collar water potential (ψpd) cannot be related to the water potential in the wettest part of the root zone, as was previously reported in literature (Hinckley and Bruckerhoff, 1975) (**Figure 9**). When gradients in soil water potential increase ψpd is closer to the driest part of the root zone as water redistribution in the soil is restricted by low unsaturated hydraulic conductivity. Disequilibrium between plant and soil water potentials was caused by the heterogeneity of soil water potential, as previously experimentally shown by Améglio et al. (1999) and Donovan et al. (2003). Root hydraulic redistribution can contribute to the disequilibrium as the nocturnal water loss prevents the recovery of plant water potential (Donovan et al., 2003). This leads potentially to the

equilibration of the system but is ultimately limited by the soilroot resistance to water flow. The largest redistribution in the model, however, takes place through the leaking split layers (**Figure 7**). For this reason, in Split 1 (SC), where the leakage caused the deeper layers to dry earlier, ψpd was very close to the potential of the dry topsoil, while in Split 3 (SC), with Compartment I being perfectly hydraulically isolated, ψpd was between the potentials of the topsoil and the soil at maximum rooting depth.

#### Determination of RSA with CT

Comparison of destructive WinRhizo scans and CT imaging showed a discrepancy of up to 27% for total root length between both methods (**Table 4**). Underestimation of root length with CT imaging had several reasons: (i) 3.5% of total root length

had diameters <0.5 mm. As a diameter of twice the resolution (voxel side length 245 and 277µm, respectively) is required for a safe detection, these roots were possibly missed by CT imaging (Koebernick et al., 2014). (ii) Roots that grow along the cylinder walls are often lost in the course of data processing, when edges of the domain have to be removed. (iii) In the "Split" setup, roots sometimes remained within the soft paraffin layers. These were eventually undetectable with X-ray CT as there is no density contrast between paraffin and roots. (iv) A possible effect of the changing soil moisture content on the segmentation cannot be excluded, since destructive measurements were only available for dry conditions at the end of the experiment. Especially at high soil moisture contents the segmentation of roots can be increasingly difficult (Flavel et al., 2012; Zappala et al., 2013). Conversely, Lontoc-Roy et al. (2006) had more difficulties segmenting maize roots from loamy sand under dry than under water saturated conditions. Our temporally repeated X-ray CT scans suggests that, for the relatively coarse roots of Vicia faba (mean diameter = 1.06 mm), water content did not strongly affect the segmentation results until the end of the experiment, when soil cracks started to form in the upper compartment of Split 1, which prevented the successful segmentation of nearby roots (**Figure 2A**).

#### Parameterization of Root Hydraulic Conductivity

Information on root hydraulic conductivities is very sparse. The use of the xylem pressure probe to determine axial and radial root hydraulic conductivities is technically very demanding, particularly for soil grown plants. Most applications refer to solution culture studies. The root hydraulic parameters for this study were derived from literature data based on experiments with lupin plants (Doussan et al., 2006) and could not be validated by direct measurements or simulation results. Thus, these parameters are a major source of uncertainty. So far, three major uncertainties could be identified:

(1) The absolute value of the conductance of the root system, Kroot, and how it differs between plants

This would affect the absolute value of simulated collar water potentials when transpiration takes place but it does not affect the predawn water potential. Thus, the conclusion that split layers reduce the collar pre-dawn water potential compared to a case where there are no split layers is not affected. The distribution of the water uptake when the soil water potential is non-uniform in the soil profile is affected by uncertainty in the absolute conductance of the root system. However, the relatively good agreement between simulated and measured soil water potentials indicates that the distribution of the root water uptake was simulated satisfactorily using the chosen root conductivities.

(2) The ratio between Kr/K<sup>x</sup>

Previous simulation studies have shown that this ratio affects the location of root water uptake (Couvreur et al., 2014). When Kr/K<sup>x</sup> is small, root water uptake occurs more uniformly along the root profile, whereas for higher Kr/K<sup>x</sup> root water uptake occurs closer to the root collar. In this study we have additional root growth, which affects the location of water uptake. Again, the relatively good predictions of the soil water potentials indicate that the root water uptake distribution was simulated quite accurately.

(3) The change of K<sup>∗</sup> r and K<sup>x</sup> over root segment age

A sensitivity analysis showed that uncertainty about the age-dependency of the root hydraulic parameters has only a small influence on the predawn water potential. However, the age dependency affects the development of the hydraulic conductivity of the total root system and hence also the xylem water potential during transpiration. Further, the root hydraulic properties used in the model could be validated and/or optimized by additional measurements of water potential in the collar or the leaves. The most reliable measurement of leaf water potential (pressure chamber, Scholander et al., 1965) is destructive and hence not suitable for measurement of changes over time. Lately developed sensors for leaf turgor (ZIM-probes, Zimmermann et al., 2013) have the potential to overcome this problem. However, for given root architecture and transpiration rates, the ranking of the collar water potentials that were simulated for our experiments will remain the same if the hydraulic properties of root segments and their dependency on age are assumed to be the same for all plants.

# Conclusion and Outlook

The initial goal was to disentangle root water uptake dynamics in a soil environment with strong water potential gradients. We addressed this question using a novel approach combining experiments, CT scanning and a simulation model. Notwithstanding the uncertainties that arise due to parameterization of the model we demonstrated the synergisms that emerge from combining split root experiments with model simulations and came to the following conclusions:


By knowing the distribution of soil and root water potentials, the combined method presented here would allow to study the direct relation between water use and root or plant growth, as was recently shown by Bao et al. (2014). Nevertheless, this is the first study in which 3-D simulations of water flow in coupled soil-plant studies were performed based on real data of the root architecture and validated against measurements of soil water potential. We did not focus on how to setup an experiment so that root properties and their uncertainty could be derived from such a setup but we rather consider the study as a proof-of-concept. In future studies, inverse modeling could be carried out to determine the root parameters and their uncertainty.

# Author Contributions

NK acquired and analyzed the experimental data, KH did the computational modeling. Both authors collaborated on the writing of the manuscript equally. EK did the setup of the computational modeling and the initial runs. JV, MJ, DV, and HV advised and commented on the manuscript.

# Acknowledgments

This work is a contribution of the Transregio Collaborative Research Center 32, Patterns in Soil-Vegetation-Atmosphere Systems: Monitoring, Modeling and Data Assimilation, which is funded by the German research association, DFG. NK is funded by Helmholtz Impulse and Networking Fund through Helmholtz Interdisciplinary Graduate School for Environmental Research (HIGRADE). The authors would like to thank Patrick Hunger for his help in conducting the experiments and Félicien Meunier for programming support in the sensitivity analysis. Further, the

# References


authors thank the four reviewers for their constructive comments that helped to improve the initial manuscript.

# Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015. 00370/abstract

Supplementary Figure 1 | Development of soil matric potential ψ<sup>m</sup> over time of the samples not used for modeling. Different colors represent measurements in different depths/compartments.

Supplementary Figure 2 | Sums of the standard uptake fraction over soil depth increments of 0.25 cm for (A) the Split 1 root system at t = 30 days and (B) the Split 3 root system at t = 34 days solved for different parameterizations of radial and axial root hydraulic conductivities.

Supplementary Figure 3 | Influence of paraffin layer on root growth: roots grow either unimpeded (left), but can also be deflected within the soft paraffin and later re-penetrate the soil. Split 1, Day 12, Layer at -5 cm, Height of image section: 13.5 mm.

Supplementary Figure 4 | Split 1 scenario SC: (A) line shows the mean and shaded areas the range (min - max) of soil water potential within each of the four soil compartments, (B) single slice at z = −12 cm showing gradients of soil water potential around the roots (black circles).


using neutron radiography and deuterated water. Vadose Zone J. 11. doi: 10.2136/vzj2011.0196


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Koebernick, Huber, Kerkhofs, Vanderborght, Javaux, Vereecken and Vetterlein. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Concepts and Analyses in the CT Scanning of Root Systems and Leaf Canopies: A Timely Summary

Jonathan A. Lafond<sup>1</sup> , Liwen Han<sup>2</sup> and Pierre Dutilleul<sup>2</sup> \*

<sup>1</sup> Département des Sols et de Génie Agroalimentaire, Université Laval, Québec, QC, Canada, <sup>2</sup> Environmetrics Laboratory, Department of Plant Science, McGill University, Montréal, QC, Canada

Non-medical applications of computed tomography (CT) scanning have flourished in recent years, including in Plant Science. This Perspective article on CT scanning of root systems and leaf canopies is intended to be of interest to three categories of readers: those who have not yet tried plant CT scanning, and should find inspiration for new research objectives; readers who are on the learning curve with applications—here is helpful advice for them; and researchers with greater experience—the field is evolving quickly and it is easy to miss aspects. Our conclusion is that CT scanning of roots and canopies is highly demanding in terms of technology, multidisciplinarity and big-data analysis, to name a few areas of expertise, but eventually, the reward for researchers is directly proportional!

#### Edited by:

Andreia M. Smith-Moritz, University of California, Davis – eLutions Integrated Systems, Inc., USA

#### Reviewed by:

Lars H. Wegner, Karlsruhe Institute of Technology, Germany Ralf Metzner, Forschungszentrum Jülich, Germany

> \*Correspondence: Pierre Dutilleul pierre.dutilleul@mcgill.ca

#### Specialty section:

This article was submitted to Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science

Received: 13 August 2015 Accepted: 24 November 2015 Published: 24 December 2015

#### Citation:

Lafond JA, Han L and Dutilleul P (2015) Concepts and Analyses in the CT Scanning of Root Systems and Leaf Canopies: A Timely Summary. Front. Plant Sci. 6:1111. doi: 10.3389/fpls.2015.01111 Keywords: computed tomography scanning, root systems and leaf canopies, scale of observation vs. scale of resolution, CT image processing and CT number analysis, structural complexity and fractal geometry, repeated plant CT scanning and statistical aspects, multidisciplinary applications

# INTRODUCTION

Many discoveries and much research progress have been made in the plant and soil sciences thanks to computed tomography (CT) scanning, from Tollner et al. (1987) and Aylmore (1993) to the 2015 Frontiers in Plant Science Research Topic "Branching and Rooting Out with a CT Scanner" (Plant Biophysics and Modeling; http://journal.frontiersin.org/researchtopic/2132/branching-androoting-out-with-a-ct-scanner-the-why-the-how-and-the-outcomes-present-and-possibly-fu), with several other pioneering studies and more modern publications in the interval (e.g., Aylmore, 1994; Tollner et al., 1994; Dutilleul et al., 2008; Flavel et al., 2012; Mooney et al., 2012). The timing was thus perfect to compile a summary that highlights the main factors, axes and issues involved in a research field which is multidisciplinary by nature and challenging by necessity.

The promise that CT scanning (originally called "computer-assisted tomography" or "CAT scanning") was able to see inside soil columns and monitor processes underground in a continuous, non-destructive manner led Aylmore (1994) to announce that the technology, designed for diagnostic purposes in the medical world, "has the potential to resolve the major controversies in soil physics and soil-plant-water relationships." The accessible challenges mentioned at the time included: understanding soil structure development and water movement in soils and its availability for plant development; documenting root system growth in 3-D space over time; and observing the behavior of soil organisms in situ. The application of CT scanning to study the branching pattern as aerial plant structure of interest and its complexity in relation to light interception by the leaf canopy is more recent (Dutilleul et al., 2005).

This perspective on the CT scanning of root systems and leaf canopies is meant to be not technically driven. The importance of a surrounding medium with similar vs. very different density, relative to the plant structure of interest, is discussed first. Then, the distinction between scale of observation and scale of resolution is emphasized because their ratio has implications for (i) the recommended equipment and its settings, and (ii) the results of graphical and quantitative analyses. When a CT scanning dataset has been collected for the root system of a crop plant or the leaf canopy of a small-size tree, the plant structure must be isolated from the surrounding medium (e.g., soil for roots) or be separated from other plant materials (e.g., seed, stem for roots; leaves for branches). Two main approaches and the companion analytical procedures will be described and illustrated with a root system example. Important mathematics and statistics questions [fractal dimension (FD) estimation, analysis of temporal repeated measures] are treated separately. To close our perspective, recent bridging experiments (alias "combo studies"), in which CT scanning plays a key role but is not the only technology applied to plants, are commented upon, under the umbrella of plant phenotyping.

# ROOT SYSTEMS, LEAF CANOPIES, AND X-RAY DOSES

In CT scanning technology, material density is essential, as it defines the CT number (CTN) value for a "voxel" (3-D extension of a pixel in 2-D space). With X-ray CT scanning, a CTN is expressed as 1000 times the ratio of the difference between the X-ray linear attenuation coefficients for the voxel and water (numerator) to the difference between the X-ray linear attenuation coefficients for water and air (denominator; the coefficient for air is, in fact, 0); its unit is HU (Hounsfield unit; Kalender, 2000). Positive and negative CTN values represent densities higher and lower than that of water, i.e., the expected amounts of X-rays absorbed are, respectively, greater and smaller than for water. Accordingly, floating plant materials have negative CTN values, between —1000 (air) and 0 (water) HU. Furthermore, the air medium surrounding a leaf canopy is lighter than the plant material, whereas the contrary is true for a root system growing in a mineral soil, for which CTN values are much greater than 0 and CTN values of any non-floating plant material (Han et al., 2008, Figure 3). The case of root systems in organic soils is the most challenging because root and soil voxels then have overlapping CTN values (Mooney et al., 2012). Depending on the type of soil, adjusting soil moisture content may be helpful; Lontoc-Roy et al. (2006) thus obtained better results for corn root systems, by CT scanning them in dry homogeneous sand and water-saturated loamy sand in a 2 × 2 factorial design.

Consequently, the amount of X-rays required to penetrate the root system and soil contained in a pot with volume V is expected to be greater than the amount required for a leaf canopy with equal volume; "penetration" means that a strictly positive portion of the X-rays emitted by the source is recorded by detectors on the opposite side in the gantry. The question of X-ray dose in plant CT scanning (especially when temporally repeated) has been the subject of a debate between Stuppy et al. (2003) and Dutilleul et al. (2005), and a welcome, informative update was recently provided by Zappala et al. (2013). It is the energy spectrum (also called "tube voltage," kV) that defines the penetrability of X-rays and their expected relative attenuation while passing through materials; higher-energy X-rays penetrate more effectively but are less sensitive to changing density than lower-energy X-rays, so a compromise must be found (Ketcham and Carlson, 2001). Tube voltage values for medical CT scanners remain in the range of 70–150 kV. On modern industrial CT scanners, settings may allow lower values (60 kV for micro-CT) and values as high as 420 kV. Radiation output increases strongly with tube voltage, but a combination of factors actually represents the radiation level delivered. Those factors also include the tube current (mA), distance from source (cm), total scan time (s), filter nature and thickness (mm). The helical scan option (i.e., when several CT images are constructed from CT scanning data acquired in one rotation) reduces the X-ray exposure time. Using the online software of McGinnis (2002– 2009) with 120 kV, 100 mA, 35 cm, 440 s, and a 3-mm Al filter (Régent et al., 2013), we calculated an X-ray dose of 6.86 Gy (6860 R) for the potato seedling of our root system example. According to Zappala et al. (2013), such X-ray dose (<30 Gy) is not damageable for plant growth or soil microbial populations.

# SCALE OF OBSERVATION VS. SCALE OF RESOLUTION

Effective work in plant CT scanning crucially depends on the scale of observation (different from the statistical concept with same name) and the scale of resolution (Ketcham and Carlson, 2001; Mooney et al., 2012). The size of the object to be CT scanned defines the scale of observation; its volume ranges from m3 to mm<sup>3</sup> , with dm<sup>3</sup> and cm<sup>3</sup> in-between. Accordingly, the scale of resolution, which is defined by the size of voxels for which CTN data are acquired, ranges from mm<sup>3</sup> to µm<sup>3</sup> , with (100 µm)<sup>3</sup> and (10 µm)<sup>3</sup> in-between. The combination of the two scales defines the type of CT: from conventional to micro-, with high-resolution and ultra-high-resolution in-between.

Interestingly, some scanning systems are said to be micro-CT, but their configuration for the reported experiments provided CT scanning data at the ultra-high resolution at best. The smallest resolution (µm<sup>3</sup> ) can be reached with synchrotron X-ray sources, and true microtomography may then allow the detection of hair roots with a diameter of a few µm. That is for a very small part (<1 cm<sup>3</sup> ) of the root system, though, and with excessively large CTN datasets to analyze if expanded (Mairhofer et al., 2012). Accordingly, synchrotron-based CT scanning has very originally been established as a technology to assess the filling status of xylem vessels and detect embolisms (Lee and Kim, 2008; Choat et al., 2015), and to unravel anatomical features of the vascular system (Kim et al., 2014). Concerning root systems and leaf canopies, medical CT scanners can provide a finer spatial resolution than anticipated in the X–Y plane of CT scanning (perpendicular to the couch), thanks to a zoom factor option (see Dutilleul et al., 2015, for leaf canopies). Even in this case, it is important that voxels be as cubic as possible, to approach the

FIGURE 1 | Photographs of root systems of (A) three 8-day-old mung bean (Vigna radiata) seedlings (left) and three 3-week-old wheat (Triticum aestivum) seedlings (right), grown hydroponically, and (B) a 5-week corn seedling (Subramanian et al., 2015) and an 11-week potato seedling (Han et al., 2009), after growing in homogeneous sand in pots, digging and washing, as examples of plant structures to be studied at scales of observation of (A) mm to cm and (B) cm to dm.

isotropy condition; the tendency for "isotropic voxels" is being generalized.

On medical CT scanning systems, a number of "fields of view" (diameters, cm, in the X–Y plane) pre-define the scale of observation (before application of any zoom factor): e.g., SS (18); S (24); M (32); and L (40). Given a 512 × 512 size of CT image and a zoom factor value (≥1.0), it is easy to calculate the X–Y dimensions of a voxel and see that high resolution is reachable with the 18-cm field of view, for a sufficiently small Z-depth. On the non-medical side, there appears to be a greater flexibility for smaller objects and finer resolution; see, e.g., the ∼2.5-cm scale of observation and ultra-high resolution in Mairhofer et al. (2012, Table 1, first Wheat column). **Figure 1A,B** here shows root systems that would be CT scanned at ultra-high resolution vs. high resolution.

# ANALYTICAL PROCEDURES FOR CT IMAGES AND ASSOCIATED CT NUMBER DATASETS

Below, the emphasis is on the isolation of roots from plant-soil CT scanning data and the subsequent construction of 3-D root system images; the cases of hydroponically-grown roots, leaf canopies and branching patterns are easier without being straightforward, the search for interfaces between leaves and branches having its own challenges (Dutilleul et al., 2015). The combination of very fast scans and automated analytical procedures may be the ultimate goal, so that CT scanning technology realizes its potential as a high-throughput technique for the quantification of roots in soils (Flavel et al., 2012), much has been done but there is still much to do regarding analytical procedures. Classically, the root

isolation problem in plant CT scanning is presented in a way that allows two antagonist approaches: "top-down"—a first set of root voxels is isolated and neighboring root voxels are joined based on some criterion; "bottom-up"—an initial group of voxels containing candidate root voxels is successively refined to remove all non-root voxels (Mairhofer et al., 2012). That problem can be posed in another way, as we explain hereafter.

Alternatively, two approaches can be developed and followed, depending on whether researchers have access to (i) the CT images only or (ii) the CTN datasets mapped into CT images. Under (i), a graphical and semi-quantitative approach based on the grayscale values in CT images (e.g., 256 tones) is followed, and the open-source package ImageJ (The U.S. National Institutes of Health; Rasband, 1997–2014) can be used. Under (ii), a root system skeleton is first traced manually through the CT images (graphical phase) and then, the skeletal roots are expanded in a neighborhood analysis involving the CTNs (quantitative phase); built-in and customized programs written in the technical computing language MATLAB (The MathWorks Inc.) have been used for this and other related purposes (Subramanian et al., 2015; see also, e.g., Koebernick et al., 2015; Paya et al., 2015). A potato root system example is presented in **Figure 2**, and we refer to the legend for more technical details. Obviously, the MATLAB-type of analysis is less automated than the ImageJ-type, even though the expansion from root system skeleton to volume is performed with a customized MATLAB program. Accordingly, the results are fragmented (ImageJ) vs. continuous (MATLAB) in **Figure 2B** (lower panel) vs. 2C (upper panel), and approach (i) can be said "bottom-up" whereas approach (ii) would be "topdown" in classical terms. On the commercial side, the general graphical features of VGStudio Max (Volume Graphics GmbH) were found useful by Flavel et al. (2012), Bao et al. (2014), and Metzner et al. (2015).

At least 100 segmentation methods are documented (see Sezgin and Sankur, 2004; Tuller et al., 2013 in the porous media and soil sciences). They include segmentation by global or local thresholding, when boundaries in grayscale or CTN values are applied to select voxels over a larger or smaller extent. And automated procedures may be aimed to avoid operator bias (Baveye et al., 2010), but are not a substitute to operator intelligence, so semi-automatic procedures may represent an acceptable compromise.

# FROM CT SCANS TO FRACTALS AND REPEATED MEASURES ANOVA

After visualization of the plant materials in the CT scanning data and isolation of the structured part, the resulting 3-D image of the plant structure of interest (root system, branching pattern) provides a basis to estimate its complexity; see **Figures 2B,C** for our root system example. However, to avoid any bias due to the thickness of roots or branches, that 3-D image must be further processed and submitted to a "skeletonization" (reduction of thickness to one voxel). Under the fractal geometry assumption (i.e., within the limits of the spatial resolution of CT scanning data, the structure repeats itself at decreasingly smaller scales; Mandelbrot, 1983), the complexity of the 3-D plant structure can be quantified by its FD, estimated with a cube-counting procedure; the higher (lower) the FD value, the higher (lower) the complexity.To ourknowledge,fractal geometry elementswere first applied by Lontoc-Roy et al. (2005)for the 3-D analysis of root system images constructed from CT scanning data; comparison to 2-D results obtained from washed roots (destructive sampling) was made in Lontoc-Roy et al. (2006). Recent, detailed plant applications of the cube-counting procedure for FD estimation can be found in Subramanian et al. (2015) for root systems and in Dutilleul et al. (2015) for branching patterns. The distinct information provided by FD, compared with root, leaf and branch traits (e.g., lengths, areas, volumes) that are not structural complexity measures, makes it a key complement to include in the quantitative analysis of plant CT scanning datasets.

If one CT image is the graphical representation of a 512 × 512 matrix of CTN values (∼250,000 entries), then 500 CT images constructed for a root system represent ∼125 million data. From the skeletonized 3-D image of a root system, however, only one estimated FD value should be retained, and from the nonskeletonized 3-D image, only one estimated root volume and one estimated total root length. That means three sample sizes of 1! Although first-order root traits can be measured individually, sample sizes are not as large as they may first seem in plant CT scanning experiments.

Zappala et al. (2013), following Dutilleul et al. (2005), established conditions in which temporally repeated CT scanning is possible with plants. Statistically speaking, temporal repeated measures on the same plant, or "subject" in general terms, do not mean increased sample size. Instead, it implies the use of "subsamples" for each subject. Since these are not random, the statistical analysis of temporal repeated measures requires an adjustment for autocorrelation and heteroscedasticity (heterogeneity of the variance) when the classical ANOVA (analysis of variance) F-tests are invalid; they suffer from an inflated rate of rejection of the null hypothesis when true. Modified ANOVA F-tests are then performed for within-subject effects (time-related effects in the ANOVA model); classical ANOVA F-tests remain valid for between-subject effects such as treatment main effects (Crowder and Hand, 1990; Dutilleul, 2011). When sample sizes are small, a mixed-model analysis is not recommended in general because of asymptotics in the estimation and testing.

# PLANT CT SCANNING COMBINED WITH OTHER TECHNOLOGIES AND METHODS OF RESEARCH

In the 1980s and 1990s, researchers had to adjust to the continuous development of a medical technology, and improve as much as possible the graphical analyses of CT scanning data in non-medical applications; Ketcham and Carlson (2001) provide an excellent review of the efforts made on the problem till then in the geosciences, including corrections to reduce beam-hardening effects in CT images of dense materials. The search for appropriate processing and analysis of CT images in an extended range of fields is fertile ground for computer-science and engineering joint research; see, e.g., the

preliminary transformation into new image files with 256 gray tones for (B).

automated graphical procedure involving multiple processing stages, proposed by Entacher et al. (2007) to generate treering profiles from the CT image showing a cross-section of the trunk.

Plant-soil CT scanning is an ideal platform to build bridges between fields, and provide supplementary information to researchers in them. Prior to the advanced applications in phytopathology by Han et al. (2008, 2009), with the common scab-inducing pathogen Streptomyces scabies and potato as the experimental crop, a medical CT system was used in soil science by Grose et al. (1996) to measure moisture content in bulk soil and the soil around roots, in order to predict suitable growth conditions for plant pathogenic fungi Rhizoctonia solani and Gaeumannomyces graminis. Because CT scanning technology is density-based, it cannot resolve solely all the research objectives in some plant studies, and other means or methods of data collection are then required. Thus, the joint use of CT scanning and phenotypic/genetic analyses allowed Bao et al. (2014) to identify a mechanism that plant roots might follow to grow toward available water. Combining shade tolerance indices from the ecological literature with leaf canopy and branching pattern traits measured from the CT scanning data collected for miniature conifers, Dutilleul et al. (2015) found differences in mean values of traits and correlations between traits depending on the leaf type, scalelike or needlelike. In a soil-water-root hydrodynamic study, Koebernick et al. (2015) included root architectures reconstructed from CT scans in a simulation model for water potentials in soil and roots in 3D and water uptake by growing roots at different depths. In a recent phytopathological application, Sturrock et al. (2015) quantified, thanks to CT scanning, the damping-off effects caused by Rhizoctonia solani on roots of wheat and oil seed rape, and related their visual assessment of the disease to pathogen DNA quantification in soil using real-time PCR.

Continued advances in CT scanning data collection and CT image analysis algorithms, for root systems under ground level and leaf canopies at interfaces with branches, will make more high-throughput applications and complete plant phenotyping possible, e.g., in greenhouse growing conditions. That is, in a breeding plan, plant structures associated with greater water and nutrient uptake from soil media and higher interception of sunlight will be revealed, but also "after the fact," new plant varieties resulting from a genomic or biotechnological improvement will see their structures characterized exhaustively. In all of this, a spatio-temporal approach based on careful use of repeated CT scanning is possible, and represents an undeniable

## REFERENCES


advantage. In closing, this is only the beginning of plant CT scanning "combos" and a large number of exciting, bridging experiments may be expected in the next years.

### ACKNOWLEDGMENTS

We are grateful to Dr. Sowmyalakshmi Subramanian and Dr. Selvakumari Arunachalam from Prof. Donald L. Smith's laboratory (Department of Plant Science, McGill University), for the preparation of the plant materials used for **Figure 1A**. We sincerely thank the two reviewers for their constructive comments on our manuscript.

the seed and peripheral organs of potato during growth using computed tomography scanning data. Trans. Am. Soc. Agric. Biol. Eng. 52, 305–311. doi: 10.13031/2013.25924


Mandelbrot, B. B. (1983). The Fractal Geometry of Nature. San Francisco: Freeman. McGinnis, R. (2002–2009). Rad Pro X-ray Device Dose-Rate Calculator Online. Cambridge, MA: Massachusetts Institute of Technology.



**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Lafond, Han and Dutilleul. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.