Introduction

Water plays an essential role as a solvent for nutrients, minerals, and other biomolecules in plant tissues1,2. To date, the inability to non-invasively image and quantify water transport directly within root tissues has been a key stumbling block for researchers seeking to understand hydrodynamics within living plant cells and tissues. Current techniques developed to monitor water uptake in roots either suffer from being indirect (tracking radiotracers3 or monitoring pressure4) or invasive (e.g., pressure chamber5, root and xylem pressure probes6, and heat pulsing7). The most promising techniques are nuclear magnetic resonance flow imaging8,9 and rapid neutron tomography10, which allow non-invasive measurement of structure and water flow within tissues. However, these techniques lack the spatial and temporal resolution to monitor water uptake kinetics at the cellular scale in root tissues.

Raman microspectroscopy (RMS) represents a non-invasive, non-destructive spectroscopic technique with the ability to investigate biological samples in situ, under physiological conditions11. By detecting inelastically scattered photons from a narrow-band laser source, RMS provides spectroscopic information at the molecular level without the need for sample preparation, fluorescent labels, or fixation. Biomolecular components (e.g., water, lipids, proteins, nucleic acids, and metabolites) appear as Raman lines at different positions inside the spectrum corresponding to their vibrational energy levels. Since the amount of energy lost by the Raman-scattered photons depends on both mass and geometry of the molecular bonds in the sampling volume, the presence of isotopes such as deuterium can be followed by monitoring distinct Raman shifts.

In this work, vibrational contrast is achieved by monitoring the appearance of deuterated water (D2O) within the Raman “silent region” (void of functional groups between 1800 and 2800 cm−1, Supplementary Fig. 1, except for triple-bond vibrations). The spatio-temporal dynamics of water flow is retrieved through an approach combining non-invasive RMS imaging following application of D2O to the root tip, and inverse modeling of D2O advection-diffusion down to the cellular scale. Our study explores whether this approach can provide non-invasive measurements of water fluxes in living plant tissues at a cellular (2 µm step size) spatial resolution and sub-second (0.3 s) temporal resolution.

Results and discussion

Water en-route to the shoot does not re-enter outer root tissues

To demonstrate the feasibility of using RMS to non-invasively image the spatial and temporal dynamics of water isotopes in plant tissues, we initially performed a proof of concept study in root tissues of Arabidopsis thaliana (Fig. 1 and Supplementary Fig. 2). RMS measurements performed as a line scan in the mature root (Fig. 1A, B, D), 11 mm from the tip, detected a pulse of deuterated water within 80 s of exposing the seedling root tip to a source of D2O (phase termed “D2O wash-in”; Fig. 1C, F). Similar temporal dynamics were detected when replacing D2O with H2O (phase termed “D2O wash-out”; Fig. 1E, G). Spatially, the D2O pulse was detected solely in inner root tissues (i.e., the endodermis and the stele, which contains water-transporting xylem vessels) and not in outer root tissues, at any point of a 64-min experiment composed of wash-in and wash-out phases (Fig. 1C, E). Hence, our results demonstrate that RMS is capable of monitoring hydrodynamics within living plant tissues. Furthermore, this result suggests that water taken up by root tips does not diffuse toward outer layers of mature root tissues in measurable quantities. This is likely to reflect the impact of both the inward advection of H2O and the diffusion barriers consisting of (i) the lignin belt—termed the Casparian strip (CS)—encircling endodermal cells within radial cell walls, and (ii) suberin deposition between the endodermal plasma membrane and the cell wall12,13 (details in Supplementary Table 1). It has been argued that the CS is porous to water since the lignin pore diameter is ~1 nm14, which is 10 times larger than a water molecule. Hence, the CS acts as a barrier to most solutes but would leave the door open to a fraction of purely apoplastic radial diffusion. Our RMS results do not support water diffusion from the stele to the cortex as A. thaliana roots take up water from their environment.

Fig. 1: Overview of the water imaging system and its outputs in line-scan mode.
figure 1

A Experimental setup for imaging hydrodynamics, showing A. thaliana under Raman measurement (laser). D2O loaded via the intact root tip is separated from the rest of the plant by a silicone grease barrier. The H2O-containing agar cap placed over the root prevents movement and desiccation. The location of panel B is indicated by the orange box. B Bright-field micrograph of intact root (scale bar 25 µm). The dotted red line indicates the location of line scans in (C–G). Note the visible xylems in this central plane of focus. C Time-course of D2O wash-in during successive line scans across the root center. Timings (in seconds) define the start of each line scan. D Root anatomy, showing the location of the confocal RMS line scan (C, E) across the full diameter, through the slightly tilted stele, with highlighted endodermis (yellow), protoxylem vessels (blue), and immature xylem (white). Spatial resolution is diffraction-limited at a 2 µm step size, temporal resolution is 300 ms/point (detector limit). E Time-course of D2O wash-out during successive line scans. Timings (in seconds) define the start of each line scan. F 3D plot showing spatiotemporal dynamics of D2O pulse arrival at the line under observation. G 3D plot showing spatiotemporal dynamics of D2O washing out of the preparation when H2O is re-introduced via the intact root. Spatial resolution is diffraction-limited at a 2 µm step size. Temporal resolution is 300 ms/point (detector limit). Dotted gray lines in C–G indicate the outer boundary of the root under observation. Source data underlying Fig. 1C, E–G are provided as a Source Data file.

Full size image

Quantified xylem water fluxes concur with stomatal aperture data

Our RMS imaging system offers the opportunity to visualize in situ xylem D2O dynamics in living roots at unique spatial and temporal resolutions, and to translate these measurements into quantitative xylem water velocities using inverse modeling. In a complementary experiment, in order to parameterize this model at an improved temporal resolution, the laser was directed at a single protoxylem vessel (identified as a local maximum of D2O signal between wash-in and wash-out phases), and D2O hydrodynamics was captured at 1 Hz from a diffraction-limited volume.

RMS imaging was used to investigate alterations of D2O dynamics and xylem water velocity in A. thaliana genotypes impaired in the integrity of their endodermal diffusion barriers. Measurements were performed on roots of wild-type plants (WT, Col-0); a WT line expressing the pCASP1::CDEF construct which degrades endodermal suberin (termed CDEF)15; and the sgn3-3 myb36-2 mutant, which displays no CS but normal suberization16. RMS imaging of “D2O wash-out” from individual protoxylem vessels revealed contrasting hydrodynamic behaviors. For example, the CDEF line needed the shortest time to complete D2O wash-out (linear slope of the wash-out curve: −5.2 × 10−3 ± 0.7 × 10−3 s−1, N = 18, significantly steeper than in WT: −3.9 10−3 ± 0.9 10−3 s−1, N = 24, Anova-1 p-value < 10−4), whereas sgn3-3 myb36-2 had almost no effect on the wash-out slope (−4.4 × 10−3 ± 1.1 × 10−3 s−1, N = 9) compared to WT (Fig. 2F, where traces are normalized to start at a unit value).

Fig. 2: Scheme of the inverse modeling loop used to extract quantitative xylem water flow rate data from RMS D2O wash-out traces.
figure 2

A Scheme of the subcellular compartments and pathways of the hydraulic model MECHA. B A. thaliana root anatomy used in the model. C Scheme of D2O wash-in and wash-out conditions. D Snapshot of a simulated D2O map for a mutant with no Casparian strip during the wash-in phase. E Fit of measured and simulated root hydraulic conductivities used to constrain the cell hydraulic parameters (WT: N = 15; CDEF: N = 12; esb1 CDEF: N = 10; sgn3 myb36: N = 10 independent plants; vertical bars show the extent of one standard deviation on both sides of the mean values). F Fit of measured and simulated D2O wash-out traces (WT: N = 24; CDEF: N = 18; sgn3 myb36: N = 9 independent plants; vertical bars and shaded areas show the extent of one standard deviation on both sides of the mean values). G Retrieved xylem water flow rates (WT: N = 24; CDEF: N = 18; sgn3 myb36: N = 9 independent plants; boxplots show 25th, 50th, and 75th percentiles, whiskers the most extreme points excluding outliers “+”). Triple stars indicate Anova-1 p-values below 0.001. Source data underlying Fig. 2E–G are provided as a Source Data file.

Full size image

From the physical point of view, xylem D2O wash-out dynamics is determined by water flow (i.e., advection, both axially and laterally) and by the mixing of H2O and D2O as they diffuse across the root composite pathways (Fig. 2A, D). Translating our wash-out traces into xylem water flow rates, therefore, required a framework capturing these hydraulic processes. We exploited a recently developed model of the root “hydraulic anatomy” termed MECHA17. This model includes explicit tissue geometries, subcellular hydraulic properties (cell walls, membranes, and plasmodesmata), and localization of diffusion barriers (Fig. 2A, B). In the Methods section, we present the implementation of its three-dimensional D2O advection-diffusion equations. MECHA can simulate D2O spatiotemporal dynamics during wash-in and wash-out cycles (Fig. 2C, D) in WT, CDEF, and sgn3-3 myb36-2 hydraulic anatomies under physiological conditions (e.g., snapshots at the laser focal point before wash-out starts in Supplementary Fig. 3, showing D2O leakage in sgn3-3 myb36-2 due to the absence of the CS). Conducting such simulations first required the estimation of hydraulic conductivity and diffusivity parameters values using an inverse modeling scheme based on an iterative search algorithm (loop indicated by curved solid arrows between panels in Fig. 2) that fine-tunes these values until convergence between measured and simulated variables. The convergence between root hydraulic conductivities (termed Lpr) simulated and measured with a pressure chamber5 under control and azide treatments (inactivating aquaporins18) (Fig. 2E, see individual points in Supplementary Fig. 4) drove the search algorithm to optimal values for the three subcellular hydraulic conductivity parameters (Supplementary Table 2). To better constrain the estimation of hydraulic conductivity parameters during this first step, we complemented the Lpr data in WT, CDEF, and sgn3-3 myb36-2 with data from the literature19 of the well-characterized esb1.1 CDEF mutant line which shows ectopic lignification in endodermal radial walls, absent suberization, and downregulated aquaporins. In a second step, diffusion parameters and xylem water flow rates were estimated in another inverse modeling scheme, this time aimed at reproducing measured xylem D2O wash-out traces in each individual replicate plant (overview of distributions in Fig. 2F and individual fits in Supplementary Figs. 5–7).

MECHA simulations suggest that the main accelerator of the D2O wash-out in Fig. 2F is the advective water uptake in the distal region, close to the tip. This seems to be because H2O flushes out the remaining D2O reserves from outer root tissues in the distal root region previously immersed in D2O. In WT, due to the presence of both endodermal barriers close to the RMS observation point (but not in the distal region), most of the uptake is distal and hence evacuates D2O reserves relatively quickly during the wash-out phase. To reach similar or steeper wash-out slopes as observed in WT (Fig. 2F), the endodermal barrier mutants thus required xylem water flow rates significantly higher than in WT (Fig. 2G, Anova-1 p-value < 10−3, see individual values in Supplementary Table 2). This result corroborates consistently higher stomatal aperture in CDEF relative to WT observed independently19, naturally generating higher xylem water flow rates.

In conclusion, this non-destructive hydrodynamic imaging approach produces meaningful quantitative results and parameters, supported by the concurrence of both stomatal and RMS observations. Furthermore, our imaging approach and model predictions suggest water transported via the root xylem does not re-enter outer root tissues or the surrounding soil when en-route to shoot tissues in A. thaliana plants with an intact endodermal diffusion barrier, thus distinguishing ‘two water worlds’.

Methods

Raman microspectroscopy

Raman spectra were recorded with a bespoke Raman microspectrometer optimized for biological samples. The instrument consists of a Raman laser (Tsunami 3960, Spectra-Physics) fitted with a 60 ps GTI and tuned at 725 nm wavelength to minimize sample autofluorescence while maximizing the Raman signal. The laser beam is passed through a clean-up filter (FF01-720/13-25, Laser2000) then directed inside an inverted microscope frame (Nikon Ti-E Eclipse), via the backport, using a dichroic mirror (FF735-DiO1-25 × 36, Semrock). A silver mirror (21010, ChromaTech) inside the microscope filter cassette is used to direct the laser into a water-immersion objective (63×/NA 1.0; Zeiss). Laser power at the sample was approximately 150 mW with a 2.2 µm laser spot size. Cells irradiated by lasers in the near-infrared region (700–1100 nm) remain viable under longer exposures at much higher laser power densities than in our experiment20. Raman photons collected by the microscope objective are passed through a long pass high-Q filter (HQ735LP, ChromaTech) upon exiting the backport of the microscope and focused using a 25 mm focal lens (AC127-025-B-ML, Thorlabs) into an optical fiber (WF50/125A 0.12NA, Ceramoptec). The optical fiber is connected to a spectrometer (77200, Oriel), equipped with an 830 lines/mm grating and a cooled deep-depletion back-illuminated CCD detector (iDus401, Andor). A high-precision automated stepper-motor stage (Prior, Cambridge, UK) fitted to the microscope is used to provide automated XY positioning. The wave-number axis of the spectrometer was calibrated using a NeAr atomic line source (IntelliCal, Princeton Instruments) and the spectral resolution was ~3 cm−1 in the 2000–3000 cm−1 region. Purpose-designed stainless steel microscope slides were fabricated that incorporated a Raman-silent MgF2 window (0.11 mm thick, 15 mm diameter) to enable the acquisition of Raman spectra.

For measurements, a freshly cleaned microscope slide (stainless steel with an MgF2 window insert) was prepared in advance by placing a water-proof barrier layer using a small syringe filled with silicone grease (Techne Ltd., part number 6101351). The barrier was dispensed onto the stainless steel part of the slide very close to the MgF2 window. A few drops of water were added on both sides of the barrier to prevent the dehydration of plants. Plant samples were placed on the slide with the leaves on the MgF2 side of the barrier and the root tip on the stainless steel side of the barrier (Fig. 1A). Then, another layer of the water-proof barrier was added on top of the root, overlapping the first barrier layer to completely seal the root. Finally, an agar cap was placed over the root on the MgF2 section to prevent sample movement and root dehydration during the measurements. For the D2O wash-in experiments, the water (H2O) on the stainless steel side of the grease barrier was replaced with D2O, while for the D2O wash-out experiments the D2O was replaced with H2O. Room temperature was 20 °C, relative humidity 43%, and light intensity 100 µmol m−2 s−1, which were all kept constant across experiments.

Line scanning (Fig. 1). Raman spectra were measured at different positions inside the plant by raster-scanning the root across the diffraction-limited laser focus in 2 µm steps (equivalent to a line of 60 points). This process was repeated 60 times (across the same line) in order to capture the temporal evolution of D2O during the wash-in and wash-out processes. The acquisition time at each position was 300 ms, yielding a total measurement time of approximately 32 min (includes raster scanning overhead). The same sample could be reproducibly scanned through multiple imaging cycles without visible photodamage. Data pre-processing consisted of three steps: cosmic ray removal, background subtraction, and D2O peak area integration. A very small fraction of individual spectra contained cosmic rays (typically < 1%). These are identified and the individual points (contaminated by cosmic rays) inside the Raman spectrum are discarded and replaced with the average value of the neighboring points. The average of the Raman spectra measured at points outside the plant represents the background spectrum (contribution from the agar, MgF2, coverslip, and microscope objective). The Raman spectrum representing each point inside the plant was obtained by algebraic subtraction of the background spectrum. Finally, the area of the D2O peak in the 2200–2800 cm−1 region was estimated using a Gaussian peak fit with a linear baseline. The obtained value is proportional to the D2O content at each point inside the plant.

Time-course (Fig. 2F). Raman spectra were measured at the same position inside the xylem (roughly 6–9 mm away from the root-hypocotyl junction) for approximately 16 min (1000 s), acquiring one measurement per second using kinetic mode acquisition inside the Andor SOLIS software. To confirm sample viability, each sample was pre-loaded with D2O via the root tip and placed on the microscope. During the D2O loading phase (approximately 15 min) a suitable measurement site inside the xylem was identified and the microscope focus was adjusted to maximize the Raman signal. After the loading phase was complete (no further increase in the Raman signal); D2O was replaced with H2O at the root tip and the D2O wash-out measurement was performed. Once complete, H2O was replaced with D2O and the wash-in measurement was performed. Data preprocessing consisted of four steps: cosmic rays removal, laser power fluctuation correction, background subtraction, and D2O peak area integration. Laser power fluctuation correction is performed for each individual Raman spectrum by normalizing every point inside the spectrum to the average value of all Raman points in the 2800–3000 cm−1 region. The background spectrum (contributions from the agar, MgF2 coverslip, microscope objective, and the plant itself) is calculated by averaging the last 20 spectra (981–1000 s) in the case of D2O wash-out, or by averaging the first 10 spectra (0–9 s) for the case of D2O wash-in. The Raman spectrum representative of each point during the time course was obtained by algebraic subtraction of the background spectrum. Finally, the transient D2O profile was obtained by using Gaussian peak fitting similar to Fig. 1.

Plant material, growth conditions, and Lpr measurements

The accession used for all experiments in this study is A. thaliana Columbia-0 (Col-0). Plants were required for three types of experiments: Raman imaging of water transport, root anatomical imaging, and root hydraulic measurements. As the last two methods are invasive, these were conducted on different plants. Only plant material used for root anatomical imaging was fixed. For Raman imaging of water transport, seeds were surface-sterilized with 50% (v/v) sodium hypochlorite and 0.01% (v/v) Triton X-100, washed three times with sterile water, and sown on plates with 0.5× Murashige and Skoog (MS) medium (Sigma) solidified with 1% (w/v) Bactoagar (Difco). After two days at 4 °C in the dark, plates were placed vertically in a growth room at 21 °C with continuous light at 100–150 µmol m−2 s−1. To obtain root sections used as templates for MECHA, seeds were treated as above but grown at 22 °C.

Root hydraulic measurements were performed with pressure chambers on 21-day-old plants cultivated for 10 days in vitro and an additional 11 days in hydroponic solution in a 16 h photoperiod growth chamber at 21/20 °C. Root hydraulic conductivity (Lpr) was determined in freshly detopped roots using a set of pressure chambers filled with a hydroponic culture medium. Excised roots were sealed using dental paste (Coltène/Whaledent s.a.r.l., France) and were subjected to 350 kPa for 10 min to achieve flow stabilization, followed by successive measurements at pressures 320, 160, and 240 kPa. Root hydrostatic conductance (Kr) was calculated by the slope of the flow-to-pressure relationship. The hydraulic conductivity was calculated by dividing Kr by the root dry weight. For sodium azide (NaN3) experiments, Lpr was determined from continuous measurement at 320 kPa. In order to compare measured and simulated Lpr expressed as hydraulic conductances per dry weight and per root surface area, respectively; values were normalized by the mean Lpr of the WT plants under control conditions, hence the arbitrary units in Fig. 2E and Supplementary Fig. 4.

Root anatomical imaging

In order to preserve accurate cell geometry, roots were fixed before acquiring and extracting the anatomy from microscope images. Root samples were fixed for 3 h at 20 °C in a solution containing 4% glutaraldehyde, 4% formaldehyde, and 50 mM sodium phosphate buffer at pH 7.2. Serial ethanol dehydration was then performed (30, 50, 70, 90, and 95% [twice]) at room temperature for 1 h at each step. Samples were embedded in Technovit 7100 resin (Kulzer) according to the manufacturer’s instructions. Sections were cut, dried onto glass slides, and stained for 20 min in an aqueous 0.1% calcofluor solution. Sections were observed under a fluorescence microscope equipped with a UV filter set.

Solution of solute transport equations from the subcellular scale

The version of the micro-hydrological model MECHA17 used in this study solves three-dimensional transient solute transport in finite-difference hydraulic networks with the implicit Crank–Nicolson method21. The hydraulic networks comprise symplastic and apoplastic “nodes”, where solute concentration is defined. Nodes are connected by hydraulic conductances (Fig. 2A) conveying water and solutes through paths including cell walls, plasma membranes, and plasmodesmata. Solutes may move along these paths with water mass flow (advection) or due to random molecular movements from a region of high to a region of low concentration (diffusion). In case the solute is deuterated water (D2O), no substantial metabolism, exclusion, or active transport applies so that advection and diffusion are the only relevant processes to consider when simulating its transport. For more elaborate simulations of the transport of major ions and associated processes in A. thaliana roots, an axisymmetric radial-longitudinal model has been developed by Foster and Miklavcic22,23. In MECHA, primary cell walls, plasma membranes, plasmodesmata, and xylem vessels are attributed to specific diffusion coefficients and hydraulic conductivities. The diffusivity of D2O within the protoplast is considered as non-limiting given its high porosity relative to membranes and plasmodesmata so that D2O concentration is assumed to freely equilibrate and thus to be uniform within each protoplast. In consequence, each protoplast only needs a single node at its center, while the apoplast is discretized in apoplastic blocks, according to cell wall geometry (parts of walls between two adjacent cells form one of these blocks). Apoplastic nodes are located at the center of each apoplastic block and at their junctions (see Fig. 2A, dotted and dashed white circles, respectively) for a total of 517 nodes in the apoplast and symplast within each two-dimensional plane in this study. In addition to connections between nodes within each transverse plane, apoplastic nodes from consecutive planes are axially connected through cell wall hydraulic conductances, while symplastic nodes are axially connected through two conductances in parallel standing for (i) plasmodesmata, and (ii) a cell wall and two plasma membrane layers. Unlike other apoplastic nodes connected by hydraulic conductances representative of primary (permeable) or secondary (hydrophobic) walls, the axial conductances of xylem vessels are calculated with the Poiseuille–Hagen law, based on their diameters determined from root sections24.

In this study, we focus on regions of the root upstream of the observation point which substantially affects xylem water composition. We, therefore, exclude the elongation zone, which does not have functional xylem vessels and radially absorbs part of the water it needs for cell elongation. For this reason, and because roots stayed fully hydrated during experiments, we assumed that cell volumes represented in the three-dimensional modeling framework did not change over time. Consequently, the root hydraulic anatomy has constant volumes for pieces of compartments, with specific diffusivities at their interfaces. An ad-hoc version of the solute advection-diffusion equation accounting for these specificities is solved. It relies on (i) the calculation of solute flow rates between compartments, (ii) the solute mass balance in each compartment, and (iii) solute boundary conditions at nodes in direct contact with the root environment. In the following elaboration, time (e.g., day), length (e.g., meter), and quantity (e.g., mole) units are referred to with the T, L, and N symbols, respectively. For enhanced clarity, scalars are represented by symbols in italics, vectors by symbols in bold italics, and matrices by symbols in bold.

As common plant root tissues do not fractionate H2O and D2O25, solute advective flow rate (Qs,a, N T−1) across each path simply equals the product of water flow rate ((Qw, L3 T−1) by D2O concentration (C, N L−3) at the origin (node i in Eq. (1), node j in Eq. (2))

$${Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{a}}}}}},{i}\to {j}}={Q}_{{{{{{\rm{w}}}}}},i\to j}{C}_{i}$$
(1)

or

$${Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{a}}}}}},{i}\leftarrow {j}}={Q}_{{{{{{\rm{w}}}}}},i\leftarrow j}{C}_{i}$$
(2)

where the subcellular water flow rate between nodes i and j is previously solved as in Couvreur et al.17.

Unlike solute advective flow, solute diffusive flow rate (Qs,d, N T−1) does not depend on the direction of water flow. It is driven by the difference of solute concentration between neighboring compartments:

$${Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{d}}}}}},i,j}={D}_{p}\frac{{A}_{i,j}}{{l}_{i,j}}({C}_{i}-{C}_{j})$$
(3)

where Dp (L2 T−1) is the D2O effective diffusion coefficient in path type p, Ai,j (L2) and li,j (L) are the cross-section and length of the path connecting node i and its neighboring node j, respectively.

In order to track D2O concentration change in time, solute advective and diffusive flow rates are balanced in the following equation for each node

$${V}_{i}\frac{\partial {C}_{i}}{\partial t}=-\mathop{\sum}\limits_{{{{{{\rm{outflow}}}}}}}{Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{a}}}}}},i\to j}+\mathop{\sum}\limits_{{{{{{\rm{inflow}}}}}}}{Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{a}}}}}},i\to j}-\mathop{\sum}\limits_{j}{Q}_{{{{{{\rm{s}}}}}},{{{{{\rm{d}}}}}},i,j}$$
(4)

where Vi (L3) is the fixed volume of the cell solution allocated to node i, outgoing (respectively incoming) solute advective flow rates negatively (respectively positively) contribute to the solute mass in compartment i, and solute diffusive flow rates are summed over all connections with neighboring compartments (j).

A consequence of the Vi factor on the left-hand side of Eq. (4) is the buffering of solute concentration in compartments with high volumes (e.g., cortical protoplast), while concentration may fluctuate faster in compartments with low volumes (e.g., cell walls). Based on Gaff and Carr26, water is assumed to occupy 69% of the primary cell wall volume and 70% of the protoplast volume.

Finally, two types of Neumann boundary conditions are set in Eq. (4) for nodes at the root surface and at the proximal end of xylem vessels. First, a solute advective inflow or outflow rate term equal to the product of water flow rate by the associated solute concentration at the origin (e.g., 0 for root surface walls bathing in H2O; equal to water inflow rate for root surface walls bathing in D2O, then subtracted from the boundary condition vector FBC at the row corresponding to the node, see Eq. (6)). Second, a solute diffusive flow rate term, which is a function of the concentration difference between the chosen boundary concentration and the root surface wall. The latter term may accelerate the propagation of the solute front and tends to zero as the compartment’s concentration equilibrates with the boundary concentration.

Equations (1–4) are solved with the implicit Crank–Nicolson method21, which uses the following approximation to ensure the convergence and accuracy of the solution

$$\frac{{{{{{{\boldsymbol{C}}}}}}}_{t+dt}-{{{{{{\boldsymbol{C}}}}}}}_{t}}{dt}\simeq \frac{1}{2}\frac{\partial {{{{{{\boldsymbol{C}}}}}}}_{t}}{\partial t}+\frac{1}{2}\frac{\partial {{{{{{\boldsymbol{C}}}}}}}_{t+dt}}{\partial t}$$
(5)

where Ct and Ct+dt (Ntot × 1, N L−3) are the vectors of solute concentrations at all root nodes at times t and t + dt, respectively (dt being a relatively small time step). Ntot is the total number of nodes (specific nodes indices are referred to as i or j in Eqs. (1–4)). In order to solve Eq. (5), one must express t as a function of t + dt, then regroup it to the left-hand side. Here, this is possible because the water flow field is at steady-state at least from time t to t + dt (note that D2O concentration does not need to be at steady-state), so the relation between C and its temporal derivative is known, see Eq. (6).

Equations (1–4) can be compacted into the following matrix operation

$${{{{{\rm{diag}}}}}}({{{{{\boldsymbol{V}}}}}})\cdot \frac{\partial {{{{{\boldsymbol{C}}}}}}}{\partial t}={{{{{\mathbf{M}}}}}}\cdot {{{{{\boldsymbol{C}}}}}}-{{{{{{\boldsymbol{F}}}}}}}_{{BC}}$$
(6)

where V (Ntot × 1, L3) is the vector of cell solution volumes associated to each root node (associated to the specific node i in Eq. (4)), V (Ntot × Ntot, L3) is a matrix containing the vector V on its central diagonal and zeros elsewhere, M (Ntot × Ntot, L3 T−1) contains the factors multiplying solute concentrations on the right hand side of Eqs. (1–3), while Neumann solute boundary conditions are concatenated in FBC (Ntot × 1, N T−1).

The matrix M is built as follows: for each node i, the outgoing water flow rates Qw,i→j are subtracted from M[i,j] (first coordinate for row index, second coordinate for column index) and incoming water flow rates Qw,i←j from each neighboring node j added to M[i, j]; for each node i, the “diffusion connectivity” \({D}_{p}\frac{{A}_{i,j}}{{l}_{i,j}}\) with each neighboring node j is subtracted from M[i, i] and added to M[i, j].

Combining Eqs. (5) and (6), Ct and Ct+dt being two distinct versions of the vector C, yields the following matrix operation

$$({{{{{\rm{diag}}}}}}({{{{{\boldsymbol{V}}}}}})-\frac{dt}{2}\cdot {{{{{\mathbf{M}}}}}})\cdot {{{{{{\boldsymbol{C}}}}}}}_{t+dt}({{{{{\rm{diag}}}}}}({{{{{\boldsymbol{V}}}}}})-\frac{dt}{2}\cdot {{{{{\mathbf{M}}}}}})\cdot {{{{{{\boldsymbol{C}}}}}}}_{t}-dt\cdot {{{{{{\boldsymbol{F}}}}}}}_{{BC}}$$
(7)

In this system, after setting initial concentrations in the vector Ct, the remaining unknowns are grouped in Ct+dt. The system is solved with the Python Scipy function “spsolve” for sparse matrices.

Model parametrization

The root hydraulic layout reproduces the anatomy of an A. thaliana primary root in the differentiation zone, with a mature protoxylem (Fig. 2B, Supplementary Fig. 2). The skeleton of the hydraulic anatomy can be digitized with the software CellSet27, as in the current study (version 1.5.1 used here), or simulated with the root anatomy simulator GRANAR28. As axial water fluxes through phloem sieve tubes are relatively small compared to xylem water fluxes, here we assume phloem flow does not affect the xylem D2O content at the RMS measurement point. Hence, for simplicity, phloem elements and their companion cells are treated as stele parenchyma cells in the simulations of D2O advection-diffusion. The cross-section is given a third spatial dimension, assuming 10−4 m long cells29 stacked axially (here we tested the model with up to 280 stacks since the longest distance between the start of the differentiation zone and the laser focal point was about 2.8 cm). As the focus of this study extends to regions of the root that transport water toward the RMS measurement point, MECHA explicitly simulates D2O advection-diffusion from the transverse section located at the maturation point of protoxylem vessels onto the RMS measurement point, about 2 cm closer to the shoot.

Two types of hydrophobic barriers are modeled in the endodermal cell walls of the selected lines. First, the CS—made of lignin— blocking apoplastic flow in radial cell walls of endodermal cells, already at the point of maturation of protoxylem vessels in WT. This barrier is absent up to 3 mm from the protoxylem maturation point in the mutant esb1.1 (and by extension esb1.1 CDEF, whose hydraulic data is used to constrain the model parameters) before an ectopic CS is formed19. The CS is completely absent in sgn3-3 myb36-2, even after the formation of the suberin lamellae, so that apoplastic flow in radial walls of endodermal cells may continue at all stages of maturation16. Second, the suberin lamellae, made of a waxy substance called suberin, is located between plasma membranes and primary walls on all sides of endodermal cells. Therefore, it blocks water flow between the apoplast and the symplast of the endodermis (short dark gray resistances in Fig. 2A). Suberin lamellae fully cover the endodermis about 1 cm after the CS in A. thaliana WT19. The suberin-degrading mutant CDEF (and by extension esb1.1 CDEF) does not display an endodermal suberin lamellae19. To isolate the impact of hydrophobic barrier locations on D2O dynamics, and as no substantial difference in root anatomies were observed between WT and mutants, the same network geometry is used for all simulations. Thus, simulated variations between different lines only stem from cell hydraulic properties adjusted at hydrophobic barrier locations as summarized in Supplementary Table 1.

Water flow and diffusion in the 10−7 m thick A. thaliana cell walls30 are limited by their hydraulic conductivity and diffusion coefficient. Assuming that they are fully hydrophobic, CS and suberin lamellae are attributed null hydraulic conductivities and diffusion coefficients, so they locally block any D2O flux. The primary cell wall hydraulic conductivity and diffusion coefficient that allow the best fits of Lpr (Fig. 2E) and wash-out traces data (Fig. 2F) are 2.5 × 10−12 m2 MPa−1 s−1 and 9.4 × 10−10 m2 s−1, respectively.

Advective water flow through non-suberized cell plasma membranes is known to occur through protein channels called aquaporins31 as well as a small part through the phospholipid bilayer of cell membranes. The phospholipid bilayer hydraulic conductivity is set to 2.6 × 10−8 m MPa−1 s−1, assuming similar physical properties as in Couvreur et al.17, while the contribution of aquaporins to the plasma membrane hydraulic conductivity that allows the best fit of Lpr data (Fig. 2E) is 5.4 × 10−7 m MPa−1 s−1, except in esb1.1 CDEF, whose aquaporin contribution is fivefold lower due to downregulated aquaporins19. D2O diffusion across the phospholipid bilayer of cell plasma membranes is also modeled, with a diffusion coefficient of 3.4 × 10−17 m2 s−1 after multiplication by the partition coefficient of water in the lipid bilayer32,33, whose thickness is only 3.0 × 10−9 m34.

Advection and diffusion within the symplast occur through plasmodesmata, whose open cross-section is set to 7.5 × 10−5 m2 over a length equivalent to twice the thickness of a primary cell wall (2 × 10−7 m)35. Based on Zhu et al.36, plasmodesmata frequencies are set to tissue-specific values reported in Supplementary Table 3.

The hydraulic conductance per individual plasmodesma and associated diffusion coefficient that allow the best fits of the Lpr (Fig. 2E) and wash-out traces data (Fig. 2F) are 1.1 × 10−19 m3 MPa−1 s−1 and 8.9 × 10−10 m2 s−1, respectively.

The last media to have its own hydraulic conductivity and diffusion coefficient is the xylem vessel. The specific hydraulic conductance of xylem poles calculated with Poiseuille–Hagen law is 2.1 × 10−14 m4 MPa−1 s−1 and the diffusion coefficient that allows the best fit of the wash-out traces data (Fig. 2F) is 1.2 × 10−9 m2 s−1, which is relatively low and barely affects the simulated traces in xylem vessels. As a matter of comparison, the self-diffusion coefficient of D2O in free water at 25 °C approaches 2.3 × 10−8 m2 s−1 37,38. Diffusion coefficients with similar values to those we found have been reported for small ionic compounds such as Na+ and K+ in A. thaliana cell walls23, which we interpret as a positive sign in view of the large uncertainties associated with the estimation of diffusion coefficients in plant tissues39. Apart from the hydraulic conductivity of membranes, which matches the range found in maize40, the subcellular hydraulic conductivities of cell walls and plasmodesmata are smaller than those found in maize with the same modeling framework17,41. Reduced diffusion coefficients and hydraulic conductivities in A. thaliana porous media are likely due to lower porosity, constrictivity, and higher tortuosity42.

Though the model parameters are all physical, no direct method allows measuring them easily, simultaneously, or non-invasively. Inverse modeling is an indirect method allowing the retrieval of model parameter values by searching for the values that best reproduce measured variables that can be simulated by the model. In this study, the measured variables are the root hydraulic conductivity (Lpr) and the xylem D2O washout trace. Different parameter values are tested as part of a “search loop” (Fig. 2A–F) until simulations and measurements converge.

With the search algorithm Multistart from the MATLAB (MathWorks, Inc.) global optimization toolbox, we find parameter values of primary cell wall hydraulic conductivity (kw), individual plasmodesma hydraulic conductance (KPD), and aquaporin contribution to plasma membrane hydraulic conductivity (kAQP) that best fit the Lpr data of WT, esb1.1 CDEF, CDEF, and sgn3-3 myb36-2, under control and azide treatments, except for sgn3-3 myb36-2 for which only control Lpr data is available (Fig. 2E, optimal parameter values in Supplementary Table 2). As the measured Lpr data has units of water flow rate per pressure differential per root mass, while Lpr simulated with the 3-D version of MECHA has units of water flow rate per pressure differential per root surface, all Lpr values are normalized by the WT Lpr under control treatment (leaving six other data points to constrain three unknowns). Note that Lpr simulations include the region of the root starting at the point of maturation of protoxylem vessels and consecutive 1.5 cm (150 stacks) shootward.

The same principle is used to find parameter values for the diffusion coefficients of primary cell walls (Dw), plasmodesmata (DPD), xylem vessels (Dx), and for the xylem axial water flow rate associated to each trace (Qxyl). Overall, simulated D2O dynamics are quite sensitive to the type of apoplastic barrier and to Qxyl, which has two consequences: (i) measured D2O traces can be expected to be sensitive to xylem water flow rate, which is a requisite for our enterprise of quantifying such fluxes in living tissues based on D2O traces, and (ii) in plants with different apoplastic barriers, similar D2O traces may have been driven by substantially different values of xylem water flow rates (as confirmed in Fig. 2F, G) and vice versa. Absolute values of D2O content would be hard to capture experimentally due to variations in the instrument transfer function between samples and fluctuations in laser power over the weeks to months during which the experiments took place. Therefore, we focus on reproducing the shape of the normalized wash-out traces (linearly transformed to start at a value of one, and end at a value of zero, 600 s after the start of the wash-out phase, Fig. 2F). Interestingly, unlike absolute fractions of D2O in simulated wash-out traces, the normalized curves do not vary substantially when scaling root length with xylem water flow rate (see Supplementary Fig. 8) while the sensitivity to the type of apoplastic barrier remains (see Supplementary Fig. 9). Thus, in order to save computing time and resources, the inverse modeling loop is conducted on a shorter virtual root of 0.15 cm (15 stacks), with 0.05 cm exposed to D2O during the wash-in phase (replaced by H2O during the wash-out phase), 0.05 cm through the barrier separating distal from proximal pools of D2O/H2O, and 0.05 cm in the proximal H2O pool. The distribution of measured and simulated wash-out traces (average and ±1 standard deviation area) are shown in Fig. 2F, as well as individually in Supplementary Fig. 5–7. Optimal parameter values that allowed reaching these fits are displayed in Supplementary Table 2. Note that as there is experimental data for 51 wash-out traces, and 54 wash-out parameters (the diffusion coefficients being conserved across simulated traces), the average number of parameters per experimental trace is 1.06.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.