- © 2015 American Society of Plant Biologists. All Rights Reserved.

## Abstract

Leaves are arguably the most complex and important physicobiological systems in the ecosphere. Yet, water transport outside the leaf xylem remains poorly understood, despite its impacts on stomatal function and photosynthesis. We applied anatomical measurements from 14 diverse species to a novel model of water flow in an areole (the smallest region bounded by minor veins) to predict the impact of anatomical variation across species on outside-xylem hydraulic conductance (*K*_{ox}). Several predictions verified previous correlational studies: (1) vein length per unit area is the strongest anatomical determinant of *K*_{ox}, due to effects on hydraulic pathlength and bundle sheath (BS) surface area; (2) palisade mesophyll remains well hydrated in hypostomatous species, which may benefit photosynthesis, (3) BS extensions enhance *K*_{ox}; and (4) the upper and lower epidermis are hydraulically sequestered from one another despite their proximity. Our findings also provided novel insights: (5) the BS contributes a minority of outside-xylem resistance; (6) vapor transport contributes up to two-thirds of *K*_{ox}; (7) *K*_{ox} is strongly enhanced by the proximity of veins to lower epidermis; and (8) *K*_{ox} is strongly influenced by spongy mesophyll anatomy, decreasing with protoplast size and increasing with airspace fraction and cell wall thickness. Correlations between anatomy and *K*_{ox} across species sometimes diverged from predicted causal effects, demonstrating the need for integrative models to resolve causation. For example, (9) *K*_{ox} was enhanced far more in heterobaric species than predicted by their having BS extensions. Our approach provides detailed insights into the role of anatomical variation in leaf function.

Leaf hydraulic conductance (*K*_{leaf}) varies widely among species (Brodribb et al., 2005; Sack and Holbrook, 2006; Sack and Scoffoni, 2013). Because the resistances inside and outside the leaf xylem (*R*_{ox}) also vary widely and are, on average across species, of a similar order of magnitude (Sack and Holbrook, 2006), both vein traits and mesophyll anatomy have potentially strong influences on *K*_{leaf}. This variation has important implications for the ecological consequences of leaf anatomy, for the coordination of water status and water flow across scales in plants, and for stomatal regulation, which may be influenced by microscale variations in leaf water potential (Buckley, 2005; Mott, 2007). However, the mechanistic basis of variation in the hydraulic conductance outside the xylem (i.e. across the bundle sheath (BS) to the sites of evaporation; *K*_{ox} [=1/*R*_{ox}]), is poorly understood (for a list of parameters and symbols used in this study, see Table I).

A strong empirical correlate of *K*_{leaf} is vein length per unit of leaf area (VLA; Sack and Frole, 2006; Brodribb et al., 2007), which is predicted to increase both leaf xylem hydraulic conductance (*K*_{x}) and *K*_{ox}: the former by providing additional parallel flow paths through the vein system, and the latter by decreasing horizontal path length for water transport from the minor veins to the sites of evaporation. High VLA may also be associated with shorter vertical path length if VLA is negatively correlated with leaf thickness, as is observed within certain species sets and lineages but not others (Noblin et al., 2008; Sack et al., 2013, 2014; Zwieniecki and Boyce, 2014). However, *K*_{ox} might be correlated with VLA due to the influence of other traits that are structurally associated with veins and are positively correlated with *K*_{leaf}, such as the size and hydraulic permeability of BS cells and the presence and size of bundle sheath extensions (BSEs). Mesophyll tissue thickness and the ratio of spongy to palisade mesophyll tissue thickness are also both correlated with *K*_{leaf} (for a comprehensive review of anatomical determinants of *K*_{leaf}, see Sack et al., 2015). Additionally, across species, mesophyll anatomy, venation architecture, stomatal conductance, and *K*_{leaf} tend to be intercorrelated (Sack et al., 2003; Aasamaa et al., 2005; Brodribb and Jordan, 2008; Carins Murphy et al., 2012, 2014; Brodribb et al., 2013; Feild and Brodribb, 2013). Thus, many of the key anatomical traits that may influence *K*_{ox} tend to be highly correlated across species (John et al., 2013), making it difficult to infer causal relationships.

Clarity on these issues requires the application of detailed anatomical data to a model that links leaf anatomy to the physics of water transport, allowing testable predictions about *K*_{ox} to be generated from alternative hypotheses about water movement beyond the xylem. Earlier models demonstrated that leaf anatomy can play a critical role in determining the sites of evaporation and major resistances within the leaf and the consequences of these features for stomatal regulation (Meidner, 1976; Tyree and Yianoulis, 1980). More recent work has led to new insights, as well as new questions, about the nature and role of vapor-phase water transport within the leaf, highlighting the need to better represent the anatomical structure of the mesophyll and surrounding airspaces in models (Rockwell et al., 2014; Buckley, 2015). The latter study made steps toward a more anatomically explicit model of leaf water flow and presented an analysis of the effects of epidermal and mesophyll anatomy on partitioning of flow among apoplastic, symplastic, and gas-phase transport modes. However, that analysis did not include several key tissues (the BS and BSEs), and it did not attempt to integrate across tissues, transport modes, and directions of flow to estimate values of *K*_{ox} comparable to experimental data. A new approach was needed to refine and test hypotheses for the influence of anatomy on water flow outside the xylem.

The objective of this study was to test hypothesized relationships between leaf anatomy and outside-xylem water transport by extending the framework of Buckley (2015) to create a new, spatially explicit model of outside-xylem water transport, MOFLO (mesophyll and outside-xylem flow), that includes all leaf tissues, including BS and BSEs. MOFLO computes *K*_{ox} and its BS and outside-BS components (*K*_{b} and *K*_{ob}, respectively) by simulating steady-state water transport outside the xylem in an areole (the smallest region of a leaf bounded by minor veins). We estimated 34 anatomical parameters from light micrographs of transverse leaf sections from 14 species diverse in phylogeny, leaf structure, and ecology and assessed the mechanistic influence of these parameters on *K*_{ox} by varying each parameter in isolation in the model while holding the others constant. We performed a range of alternative simulations to address uncertainty in parameters that could not be confidently estimated by light microscopy (LM). We used these simulations to address five interrelated questions. (1) Where are the major resistances located outside the xylem (i.e. in which tissues, and in which type of flow pathways), and particularly, how much resistance is contributed by the BS? (2) How do BSEs affect *K*_{ox}? (3) How do other cell and tissue anatomical traits influence *K*_{ox} and *K*_{leaf}? (4) Can these influences explain previously described correlations of anatomical traits, and particularly VLA, with *K*_{leaf}? (5) What are the roles of gas-phase transport, temperature, and vertical temperature gradients in determining *K*_{ox}?

## RESULTS

### Comparison of Simulated Values of *K*_{ox} across Species with Measured Values

Observed *K*_{ox} ranged from 3.5 to 54.3 mmol m^{−2} s^{−1} MPa^{−1} across the eight species for which measurements were available (Table II; Fig. 1). The mean and median simulated *K*_{ox} across those eight species (16.8 and 13.6 mmol m^{−2} s^{−1} MPa^{−1}, respectively) were greater than, but of similar order of magnitude to, the mean and median observed *K*_{ox} (11.9 and 5.4 mmol m^{−2} s^{−1} MPa^{−1}, respectively; Table II; Fig. 1). For seven of the eight species measured, the observed values fell between the low and high simulated values from simulation set 1, which used a wide span of values for each of the six unknown parameters of leaf design (Table III). The exception was *S. canariensis*, for which measured *K*_{ox} exceeded the high simulated value. The measured and modeled values of *K*_{ox} were uncorrelated across species (*P* > 0.05; data not shown), which was to be expected, given that our modeled estimates of *K*_{ox} are based on assumed values for several parameters whose true values are unknown and may differ across species.

### Modeling the Water Potential Drawdown outside the Xylem

Figure 2 shows an example of the simulated distribution of water potential drawdown outside the xylem (δψ) in a transverse section of a radially symmetrical areole for one species, *C. diversifolia*, using default values for all parameters (Tables III–IV). The drawdown increases from the BS (left-hand edge, in rows 18–22) to the lower (abaxial) epidermis at the center of the areole (bottom right corner). Although the drawdown exceeds –2.2 MPa, the volume-weighted average drawdown is only –0.6 MPa (or –0.63 MPa excluding the BS itself). One reason for this difference is that much of the leaf’s water is in palisade mesophyll, which is outside of the main pathways for water flow from the xylem to the transpiring epidermis and consequently experiences little drawdown. In this example, simulated *K*_{ox} was 7.9 mmol m^{−2} s^{−1} MPa^{−1}, *K*_{b} was 19.1 mmol m^{−2} s^{−1} MPa^{−1}, and *K*_{ob} was 13.4 mmol m^{−2} s^{−1} MPa^{−1}.

### Partitioning Hydraulic Resistance outside the Xylem

Across all 14 species, simulated *K*_{ox} ranged from 4.0 to 28.6 mmol m^{−2} s^{−1} MPa^{−1}, with a median of 9.0 and mean of 11.8 (Table V). Simulated *K*_{ob} varied from 4.8 to 47.4 (median, 13.2), and *K*_{b} varied from 7.1 to 136 (median, 41.6). On average, for default parameter values, most outside-xylem resistance occurred outside the BS: although the BS contribution ranged from 12% to 71%, the median was 18%.

### The Importance of Tissue Types and Transport Modes in outside-Xylem Water Transport

Tissue types and transport modes varied widely in their contributions to outside-xylem water transport. On average across species, the bulk conductivity [flow per unit of water potential gradient per unit of bulk tissue area; mol s^{−1} m^{−2} (MPa m^{−1})^{−1}] was greatest in the lower epidermis and lowest in the palisade mesophyll (for horizontal transport), followed closely by spongy mesophyll transport (Fig. 3). Bulk conductivity in BSEs and across the BS itself were more than double that of the spongy mesophyll (Fig. 3). Apoplastic pathways provided most transport in all tissues, although transmembrane and (isothermal) gas-phase transport modes together contributed nearly half of the bulk conductivity in the spongy mesophyll. (The roles of anisothermal vertical gas-phase transport driven by temperature gradients, and of temperature itself, are discussed further below.)

### Effects of Changes in Six Unknown Parameters: Apoplastic Pore Diameter, Cell Membrane Permeability, BS Suberization, Palisade Connectivity, Cell Wall Thickness, and Vertical Temperature Gradient

*K*_{ox} was highly sensitive to the values of parameters that could not be estimated confidently, which we refer to here as unknown parameters (listed in Table III). Under default values for other parameters, *K*_{ox} increased by 668% when the *R*_{a}) increased from 3 to 10 nm and decreased by 71% when *R*_{a} decreased from 3 to 0 nm (Fig. 4). However, *K*_{ox} was less sensitive to the *P*_{m} under default values for other parameters, increasing only 52% when *P*_{m} was increased 4-fold from 40 to 160 μm s^{−1} and decreasing just 18% when *P*_{m} was reduced from 40 to 0 μm s^{−1} (Fig. 4). However, if the BS apoplast was assumed to be suberized, *R*_{a} and *P*_{m} had similar influences on *K*_{ox} (Fig. 4).

The fraction of horizontal palisade surface area in contact with adjacent palisade cells (*f*_{cph}) had little effect on *K*_{ox}, which increased only 45% when *f*_{cph} increased from 0% and 100% of the apparent value measured by LM (i.e. when the ρ_{fcph} increased from 0 to 1); furthermore, most of this increase occurred below ρ_{fcph} = 0.2 (Fig. 5). Cell wall thickness was far more important in determining *K*_{ox}: *K*_{ox} increased by 400% when cell wall thicknesses used in simulations were increased from 20% to 100% of the values determined by LM (i.e. when the ρ_{ta} was increased from 0.2 to 1; Fig. 5).

Mean *K*_{ox} across species was strongly enhanced by the presence of a vertical temperature gradient within the leaf: doubling the gradient from its default value of 0.1°C increased *K*_{ox} by 75%, and eliminating the gradient reduced *K*_{ox} by 27% (Fig. 6; note that 0.1°C was the average temperature drop from the point of maximum temperature to the lower epidermis across species; in practice, we used the same gradient [4.6 × 10^{−4}°C μm^{−1}] for all species, so that the absolute temperature drop varied across species in relation to leaf thickness). Comparing these simulations with another that excluded gas-phase transport altogether, we calculated that the average gas-phase contribution to *K*_{ox} increased from 16% to 65% as the temperature gradient increased from 0°C to 0.2°C.

We also assessed the effect of temperature itself, as distinct from temperature gradients. *K*_{ox}, *K*_{b}, and *K*_{ob} all increased strongly with temperature (Fig. 7), but the relative increases in *K*_{ox} and *K*_{ob} were far greater than that for *K*_{b}: *K*_{ox} and *K*_{ob} increased by 286% and 378%, respectively, as temperature increased from 10°C to 40°C, whereas *K*_{b} only increased by 81% over the same temperature range (note that the effect of temperature on *K*_{b} results only from changes in the diffusivity of liquid water in water, because our model did not include any gas-phase water transport across the BS due to the lack of airspaces in the BS).

Changes in the six unknown parameters did not, in most cases, result in substantial changes in the partitioning of hydraulic resistance outside the xylem, which proved robust across most simulations: less than 25% of outside-xylem resistance was contributed by the BS under any tested combination of values for *R*_{a} and *P*_{m}, provided the BS apoplast was not assumed to be suberized (Fig. 8). When a BS Casparian strip was included in the simulations (thus preventing apoplastic transport across the BS), the BS accounted for nearly 40% of total outside-xylem resistance under default values for other parameters and up to 75% for high *R*_{a} (10 nm) and low *P*_{m} (20 μm s^{−1}; Fig. 8). However, changes in ρ_{ta} had little effect on the percentage of outside-xylem resistance in the BS, which decreased from 25.2% to 23.2% as ρ_{ta} increased from 0.2 to 1 (data not shown).

### Functional Consequences of Known Anatomical Traits on *K*_{ox}: VLA, Vein Positioning, Leaf Thickness, BSEs, and Leaf Airspace Fraction

Table VI lists standardized slopes for linear regressions between each anatomical parameter and modeled *K*_{ox}. By far, the strongest influence of leaf anatomy on *K*_{ox} was that of VLA: *K*_{ox} increased 121% with a doubling of VLA (Fig. 9), due in part to the effect of VLA on BS surface area per unit of leaf area (which affects *K*_{b}; Fig. 9) and in part to the fact that VLA reduces the horizontal pathlength for water transport to the transpiring epidermis (which affects *K*_{ob}; Fig. 9). The pathlength effect was stronger than the BS area effect (increasing VLA increased the proportion of outside-xylem resistance in the BS; data not shown). The VLA effect was over 3 times stronger than the next strongest anatomical effects: the increase in *K*_{ox} resulting from greater relative proximity of the vascular bundle to the abaxial epidermis (represented here as the ratio of distances between the BS and the upper versus lower epidermis; Fig. 10) and the decrease in *K*_{ox} caused by increasing spongy mesophyll cell radius (Fig. 11). (The spongy cell radius effect arises because of the dominance of apoplastic transport: if cell radius increases without a concomitant increase in wall thickness, the apoplastic fraction of available transport area declines.) For both of the latter effects, *K*_{ox} changed by approximately one-third with a doubling of the parameter value (Table VI).

Across species, *K*_{ox} was uncorrelated with leaf thickness, and leaf thickness had a smaller mechanistic influence on *K*_{ox} (doubling thickness reduced *K*_{ox} by 18%; Fig. 10; Table VI) than the relative proximity of the vascular bundle to the lower epidermis. The lack of a cross-species correlation between leaf thickness and modeled *K*_{ox} may partly reflect a positive correlation between leaf and cell wall thicknesses in our species (data not shown), which would tend to counteract the effect on *K*_{ox} of increased vertical pathlength in thicker leaves.

Eight of our 14 species were heterobaric (they possessed BSEs), and six were homobaric. We assessed the mechanistic effect of BSEs on *K*_{ox} by comparing standard simulations with another set of simulations in which BSEs were replaced with mesophyll tissue in the model. These simulations found that BSEs directly increased *K*_{ox} by 10% on average across the eight heterobaric species (Fig. 12). However, *K*_{ox} was 34% greater in heterobaric than homobaric species (Fig. 12), which suggests that the enhancement of *K*_{ox} in heterobaric species is mostly due to factors other than the BSEs themselves.

### Correlations of Anatomy across Species with *K*_{ox}: Divergence from Mechanistic Relationships

In each of the cases described above, the correlation between each parameter and the simulated values of *K*_{ox} across species was in the same direction as the mechanistic effect. The opposite was true for several other parameters, however. For example, the mechanistic effect of the fraction of spongy mesophyll cell area in contact with adjacent cells (*f*_{cs}) was positive (simulated *K*_{ox} increased 24% with a doubling of *f*_{cs}), whereas the correlation across species was strongly negative (Fig. 11; Table VI). The converse was true for the ratio of palisade to spongy mesophyll thickness: the mechanistic effect of this ratio on *K*_{ox} was weakly negative, but the correlation across species was strongly positive (Fig. 11; Table VI). Spongy mesophyll airspace fraction (*p*_{s}) also had a positive mechanistic influence on *K*_{ox} (Fig. 11), with *K*_{ox} increasing 35% as *p*_{s} increased from 0.1 to 0.6, whereas these variables were uncorrelated across species (Table VI).

## DISCUSSION

We elucidated and addressed key hypotheses for the anatomical basis of *K*_{ox} by applying measured variations in leaf anatomy across a set of very diverse species (Table VII) to a novel computational model, MOFLO. Our analysis led to several predictions consistent with previous work, but equally, to a number of surprising novel predictions. We addressed several questions, discussed below.

### Where Are the Major Resistances outside the Xylem?

Our simulations converged in showing that most resistance beyond the xylem occurs in the spongy mesophyll and that the BS contributes a minority of outside-xylem resistance. The spongy mesophyll is intrinsically more resistive than other tissues, because its airspace fraction is high (averaging 37% across species, nearly twice that of the palisade) and its cell-to-cell connectivity is low (an average of 26% of spongy cell surface is in contact with other cells), both of which reduce the area effectively available for liquid-phase flow. Our calculations suggest that the epidermis is over three times as conductive on a bulk area basis than spongy mesophyll, on average, across our 14 study species. Only horizontal transport in the palisade has a lower bulk conductivity than the spongy mesophyll, but this has little impact on *K*_{ox} because most water flows through the spongy mesophyll in hypostomatous species (12 of the 14 species in this study).

The true contribution of the BS to outside-xylem resistance remains somewhat ambiguous due to uncertainty about the occurrence of a suberized layer (Casparian strip) in BS cell walls. Such a strip would greatly reduce apoplastic conductivity across the BS, rendering the BS analogous to the root endodermis, and its presence is one of the major outstanding questions in leaf design. Previous studies have suggested a BS Casparian strip in certain grass species, *Plantago* spp., and at least several other taxa (Lersten, 1997; Mertz and Brutnell, 2014), and the expression of similar genes during development in BS and root endodermis suggests functional similarities (Slewinski et al., 2012). In any case, even when the model was modified to include a BS Casparian strip, the average BS contribution to outside-xylem resistance only increased from 10% to 37% under default values for other parameters. Thus, we tentatively conclude that the BS contributes a significant but minority share of outside-xylem resistance.

### Does Liquid Flow outside the Xylem Follow Apoplastic and/or Transmembrane Routes?

Previous studies using staining or conceptual modeling have reached differing conclusions about the relative importance of transport across living cells or around them in the apoplast. Apoplastic tracer studies (Canny, 1986) and the discovery of aquaporins (Agre et al., 1993; Chrispeels and Agre 1994) have promoted the view in recent years that transmembrane flow may dominate outside-xylem transport (Tyree et al., 1981, 1999; Sack and Tyree, 2005), at least in the light, when aquaporins may be activated (Cochard et al., 2007). However, a theoretical study by Buckley (2015) that used membrane permeability values from published studies carried out on illuminated leaves concluded that apoplastic transport should dominate. MOFLO extends upon that study and similarly predicted that that apoplastic bulk flow contributes the majority of *K*_{ox} (68% on average across species), thus dominating both transmembrane and gas-phase pathways under most conditions. This is due to the intrinsically greater efficiency of apoplastic bulk flow than either liquid- or gas-phase diffusion. Although our LM-based measurements of cell wall thickness (which strongly determine apoplastic conductance) were much greater than most published estimates for other species, this does not explain the model’s predictions concerning apoplastic transport, because by default we reduced our LM-based estimates of cell wall thicknesses by 80% before applying them to the model (ρ_{ta} = 0.2). Transmembrane pathways contributed only 19% of *K*_{ox} on average, and this fraction was smaller still (6%) if LM-based cell wall thicknesses were used. (The contribution of gas-phase pathways is discussed below.)

These conclusions assume that bulk flow in the apoplast can be modeled using Poiseuille’s law, which is derived from the Navier-Stokes equations of continuum fluid mechanics. Continuum hydrodynamics is valid provided the flow channels are large relative to the chemical species. The relevant size measure for liquid water molecules in this context is the lattice spacing, which is approximately 0.31 nm. Eijkel and Van Den Berg (2005) note that “friction is seen to increase from the macroscopic [continuum-derived] value when the separation between two surfaces becomes less than, roughly, ten molecular layers,” or approximately 3 nm in this case. This is identical to the low end of the range estimated by Buckley (2015) for the diameter of channels for water flow created by spaces between adjacent microfibrils or bundles of microfibrils in the apoplast (3–20 nm), based on published measurements of cell wall microstructure (McCann et al., 1990; Fleischer et al., 1999; Fahlén and Salmén, 2005; Kennedy et al., 2007), which suggests that the continuum approximation is probably reasonable for apoplastic transport.

The framework developed by Buckley (2015) included a term for diffusive resistance across the cellular interior (transcellular resistance) in series with transmembrane resistance. Further thought and discussions with colleagues led us to conclude that any water transport across the cellular interior probably occurs mostly by bulk flow, provided that the flow area consists of channels much greater than the 0.31-nm lattice spacing of water. Even if those channels had a typical radius similar to those in the adjacent cell walls, transcellular resistance would be on the same order of magnitude as apoplastic resistance (and thus, far smaller than transmembrane resistance) if the transcellular area available for water flow were similar to the apoplastic flow area. Regardless, if this is incorrect and transcellular resistance is large, that would only strengthen our conclusion that apoplastic transport dominates outside-xylem water transport.

### The Effect of BSEs

Previous studies that inferred the effect of BSEs on *K*_{leaf} from anatomy, simpler hydraulic models, *K*_{leaf} responses to light, and stomatal responses to evaporative demand in heterobaric versus homobaric species have hypothesized that BSEs are a major route for water flow from the veins to the epidermis to the stomata (Wylie, 1952; Scoffoni et al., 2008; Buckley et al., 2011; Sommerville et al., 2012; Zsögön et al., 2015). MOFLO allowed us to directly quantify the effect of BSEs on *K*_{ox} by replacing BSEs with mesophyll tissue in the model. The results suggested that BSEs enhance *K*_{ox} by an average of 10% across the eight heterobaric species in this study. However, simulated *K*_{ox} was 34% greater in these species than in the six homobaric species. This finding suggested that the presence of BSEs is correlated with one or more other traits that also enhance *K*_{ox}. The only anatomical parameter that differed significantly between heterobaric and homobaric species in our data set was spongy mesophyll cell radius (*r*_{s}; *P* < 0.05, two-tailed Student’s *t* test with unequal variances): *r*_{s} was greater in homobaric species (21 ± 3 versus 12 ± 2 μm). This is consistent with our mechanistic trait analysis, which predicted that *K*_{ox} should decrease by 30% for a doubling of *r*_{s} (Table VI).

### Effects of Cellular Dimensions on *K*_{ox}

Most individual anatomical traits affected *K*_{ox} only weakly. The major exceptions involved spongy mesophyll anatomy, which had much larger influences than palisade anatomy because most of our study species (12 of 14) were hypostomatous, so little water transport occurs through the upper half of the leaf. The apparent effect of *r*_{s} in our trait analysis arose because, when all other parameters are held constant, increasing *r*_{s} increases the transmembrane fraction of the total cross-sectional area available for flow, which decreases the apoplastic fraction, in turn decreasing *K*_{ox}. However, *r*_{s} is often correlated with spongy mesophyll cell wall thickness across species (John et al., 2013), which would tend to reduce the direct effect of *r*_{s}. Another explanation for the similarity between the correlative and mechanistic relationships that we found between *r*_{s} and *K*_{ox} (Fig. 10B) is that *r*_{s} was negatively correlated with VLA and with the relative proximity of vascular bundles to the lower epidermis (*r*^{2} = 0.25 and 0.61, respectively; *P* < 0.0001 for both), both of which had positive mechanistic effects on *K*_{ox}, as discussed below. A similar negative correlation between VLA and the sizes of mesophyll and epidermal cells was reported previously to hold across species of Proteaceae by Brodribb et al. (2013).

### Effects of VLA, Leaf Thickness, and Distance from Vascular Bundles to Epidermis

The specific role of VLA in increasing outside-xylem flow has been a topic for debate. Sack and Frole (2006) suggested that higher VLA led to shorter horizontal flow distances, increasing *K*_{leaf}. This was also found by Brodribb et al. (2007), who additionally hypothesized that a shorter vertical distance between vein and epidermis would also increase *K*_{leaf}. Indeed, because high-VLA leaves are often thinner as well, a correlation that has been hypothesized to be optimal for water transport based on modeling using artificial leaf assemblies (Noblin et al., 2008), a high VLA would also correspond to such shorter vertical distance. Brodribb et al. (2007) combined the hypothesized effects of horizontal and vertical distances in their variable *D*_{m}, representing a diagonal distance from veins to epidermal evaporating sites, and reported a strong correlation between this diagonal distance and leaf hydraulic resistance, which was mostly driven by VLA. However, Sack et al. (2013) suggested that greater leaf thickness should contribute to higher *K*_{ox}, given the greater number of parallel pathways for horizontal transport to the sites of evaporation, provided those sites are distributed throughout the leaf. MOFLO allowed us to test these putative mechanisms. We found that increasing VLA, reducing total leaf thickness, and reducing the relative distance of vascular bundles from the lower epidermis all increased *K*_{ox} in the model because of their effects on reducing flow pathlengths, although the effect of VLA was by far the strongest and that of leaf thickness the weakest of the three. The model found that VLA affects *K*_{ox} in two ways: by increasing BS surface area per unit of leaf area, which affects *K*_{b}, and by decreasing horizontal pathlength, which affects *K*_{ob}. Although both effects were quite strong, the pathlength effect was stronger (*K*_{ob} and *K*_{b} increased 113% and 94%, respectively, with a doubling of VLA; Fig. 9).

The model also found a negative mechanistic effect of vertical pathlength (as influenced by either total leaf thickness or relative vein-to-epidermis distance), but these effects were only one-sixth and one-third as strong, respectively, as the horizontal distance effect of VLA (Table VI). The main reason for the smaller effect of changes in vertical pathlength (i.e. of leaf thickness) than horizontal pathlength (i.e. VLA) on *K*_{ox} is that adding vertical layers simultaneously also reduces the horizontal resistance by providing additional parallel pathways for horizontal transport (Sack et al., 2013). In contrast to its mechanistic effect, we found that leaf thickness was not significantly correlated with simulated *K*_{ox} across our species, due to compensating effects of other parameters that covaried with leaf thickness. For example, leaf thickness was strongly and positively correlated with cell wall thickness in each tissue type (*r*^{2} between 0.33 and 0.69, *P* < 0.0001 in all cases; data not shown), all of which had strongly positive mechanistic effects on *K*_{ox} (Table VI). These results verify that the often-observed correlation between *K*_{leaf} and VLA is mechanistic in origin (Sack and Frole, 2006; Brodribb et al., 2007; Brodribb and Jordan, 2008; Carins Murphy et al., 2012, 2014; Feild and Brodribb, 2013), and they further suggest that the horizontal pathlength component of the VLA effect is more important than the vertical component.

### The Role of Gas-Phase Transport and Vertical Temperature Gradients

Recent work has raised the possibility that gas-phase water transport contributes a substantial fraction of the total conductance for water movement through the mesophyll, perhaps comparable in magnitude to that provided by liquid-phase pathways, particularly for vertical transport in the presence of large vertical temperature gradients (Rockwell et al., 2014; Buckley, 2015). Our analysis extended that work by providing, to our knowledge for the first time, an integrated measure of *K*_{ox} that includes both horizontal and vertical components of gas-phase transport, all in the same leaf area-based hydraulic conductance units. The model found that gas-phase transport contributed an average of 39% of *K*_{ox} across species under default conditions (which include a baseline temperature of 25°C and a vertical temperature gradient of 0.1°C). This rose to 65% for a gradient of 0.2°C and fell to 16% for zero gradient. Thus, we conclude that the contribution of vapor transport within the leaf to the apparent conductance for water transport can be quite substantial.

This has several implications for interpreting leaf function and gas exchange. First, it implies that the generation of vertical temperature gradients by preferential absorption of light near the upper leaf surface can enhance *K*_{ox} greatly: by over 20% for 0.1°C gradients or 40% for 0.2°C gradients. This corresponds to average 16% and 31% enhancements of *K*_{leaf}, respectively, across the eight species in our data set for which we measured *K*_{ox}. These effects could contribute to the observed effects of light on *K*_{leaf}, in addition to other mechanisms such as increased aquaporin activity (Cochard et al., 2007; Scoffoni et al., 2008; Voicu et al., 2009).

Second, a major role for vapor transport implies that a great deal of water may evaporate from cells deep within the leaf. This contrasts with some earlier conclusions (Tyree and Yianoulis, 1980) that the great majority of evaporation occurs from cells very close to the stomatal pore, but it is consistent with the conclusions of Boyer (1985) based on measurements of vapor diffusion pathlength by Farquhar and Raschke (1978). The question of where evaporation occurs within the leaf has remained one of the most challenging and critically important in plant water transport for decades (Meidner, 1983; Barbour and Farquhar, 2004) and demands further discussion here. In the context of water transport, evaporation represents a shift of water from a liquid pathway to a gas-phase pathway. Water flow will distribute itself across pathways so as to minimize total resistance; therefore, some water will switch from a liquid to a gas-phase pathway whenever the gas-phase conductance increases relative to the liquid-phase conductance (Buckley, 2015). Thus, evaporation should occur wherever the gas-phase fraction of total conductance increases along a trajectory of flow (a pathway normal to isoclines of water potential). That fraction increases substantially in three areas: (1) at the outer margin of the BS (where the fraction rises from 0 to some positive value when water first encounters airspaces in the leaf); (2) at the boundary between palisade and spongy mesophyll (where the gas-phase fraction increases due to increasing tissue airspace fraction and decreasing vertical liquid-phase conductance); and (3) at open stomatal pores, where the gas-phase fraction approaches 100% (because all water exits the leaf as vapor). This suggests that evaporation is clustered in three locations in hypostomatous leaves: the BS, the upper spongy mesophyll, and surfaces immediately adjacent to open stomata. A similar argument would apply to amphistomatous species with spongy mesophyll in the center of the leaf, except that the prevailing direction of water flow would be from spongy into palisade mesophyll, implying that condensation rather than evaporation would occur at the spongy/palisade transitions. The liquid-phase share of transport from those regions to the transpiring epidermes would thus be greater in amphistomatous species than in hypostomatous species (due to the greater liquid conductivity and smaller porosity of palisade as compared with spongy mesophyll), which in turn implies that a greater share of evaporation would occur from surfaces very close to the stomata in amphistomatous species.

### The Role of Temperature Itself

The direct effect of temperature on *K*_{ox} (independent of temperature gradients) was also substantial in the model: under otherwise default parameter values, *K*_{ox} increased 25% as leaf temperature increased from 25°C to 30°C and increased 233% for an increase from 25°C to 40°C. This effect arises partly from the temperature dependence of liquid-phase conductivities (chiefly due to decreasing dynamic viscosity) but more so from increasing gas-phase conductivities (due to strong increases in both the molecular diffusivity of water vapor in air and the saturation vapor pressure). These direct temperature effects could further contribute to the light responses of *K*_{leaf} in nature, where temperature usually increases with the absorption of sunlight. A direct increase in *K*_{ox} with temperature could also help to sustain turgor when water loss increases as a result of leaf warming rather than drying of the air; such an effect may also help to explain positive correlations reported between *K*_{leaf} and transpiration rate (Simonin et al., 2015) in cases where changes in transpiration are temperature driven.

### Implications for Stomatal Sensing of Leaf Water Status

Our model suggested that large water potential gradients could occur between the xylem and the most distal epidermal tissues: in the example shown in Figure 2 (for *C. diversifolia*), the drawdown from the xylem to the lower epidermis at the center of the areole was 3.7 times greater than the average drawdown outside the xylem. This ratio varied across species, reaching 6.3 in *M. grandiflora*, and it was substantial at 2.2 even in the amphistomatous species *H. annuus*. These results support the hypothesis that a transpiring epidermis (and the stomatal guard cells embedded therein) may experience far greater swings in water potential in response to changes in transpiration rate than one would infer from changes in bulk leaf water potential. This may help to reconcile isohydric behavior (near homeostasis in leaf water potential [ψ_{leaf}]) with a mechanism for stomatal responses based on a feedback response to changes in water potential somewhere in the leaf (Sperry, 2000; Buckley, 2005). The large drawdowns predicted by the model also suggest that the upper and lower epidermes are in effect hydraulically sequestered from one another, which may help to explain the observation that stomata at one surface appear only minimally responsive to changes in transpiration rate at the other surface (Mott, 2007). We tested this idea directly in MOFLO by tripling the transpiration rate at the upper surface of a simulated *H. annuus* leaf while holding transpiration constant at the other surface; the resulting change in water potential at the center of the areole in the upper epidermis was 4.3 times greater than in the lower epidermis (Fig. 13).

## CONCLUSION

Our novel analyses provide, to our knowledge for the first time, quantitative integration of the effects of leaf anatomy on water flow outside the xylem, in terms directly comparable to experimental data. Our model confirmed some earlier predictions about the relation of *K*_{ox} to leaf anatomy, including that VLA is the strongest anatomical determinant of *K*_{ox} and that BSEs and thermally driven vapor transport through spongy mesophyll can enhance *K*_{ox}, but also provided novel insights, including that the BS probably contributes a minority of outside-xylem resistance, that higher *K*_{ox} in heterobaric species is mostly due to parameters other than BSEs, that vapor transport may constitute a majority of *K*_{ox} when large vertical temperature gradients exist in the leaf, and that many cross-species correlations between *K*_{ox} and leaf traits are not mechanistic in origin. Our model provides strong insights into the coordinated function of the living leaf, a tool to explore the implications of variation in leaf anatomy, and a baseline for future trait analyses.

## MATERIALS AND METHODS

### Empirical Measurements of *K*_{ox}

We determined *K*_{ox} from measured *K*_{leaf} and *K*_{x} for eight of our 14 study species (Table II). *K*_{leaf} was obtained from whole-leaf hydraulic vulnerability curves previously published using the evaporative flux method (Scoffoni et al., 2012, 2015). Because *K*_{leaf} declines as water potentials become more negative, we calculated for each species the average *K*_{leaf} for the interval of leaf water potential near full hydration (we used 0 to −0.3, 0 to −0.5, or 0 to −1 MPa, depending on species, to capture the interval before a strong decline in *K*_{leaf}: *n* = 5–12). *K*_{x} was obtained as described previously (Scoffoni et al., 2015) using the vacuum pump method. Briefly, minor veins of fully hydrated leaves were cut under water over a light bench to ensure that no major veins were severed. Cuts were made between about 95% of tertiary veins, yielding 5 to 33 cuts mm^{−2} depending on sample size (larger leaves have their major veins spaced farther apart, so that fewer but longer cuts were made; Sack et al., 2012; Scoffoni and Sack, 2015). These cuts were enough for water to move directly out of minor veins and not through outside-xylem pathways (Sack et al., 2004; Nardini et al., 2005). After minor veins were cut, leaves were connected by tubing to a water source on a balance and placed in a vacuum chamber. A steady flow rate was determined for five levels of partial vacuum (0.06, 0.05, 0.04, 0.03, and 0.02 MPa). *K*_{x} was calculated as the slope of the flow rate against pressure, corrected for leaf temperature, normalized for leaf area, and averaged (*n* = 5–11). *K*_{ox} was calculated using Equation 1, and se values were obtained from propagation of error:(1)We note that estimates of *K*_{ox} thus depend on the accuracy of *K*_{leaf} values. In particular, the evaporative flux method requires steady-state transpiration and stable leaf water potential to enable the determination of *K*_{leaf}. We followed the procedure tested and established for a wide range of species in previous work (Scoffoni et al., 2008, 2012; Pasquet-Kok et al., 2010; Guyot et al., 2012). In measuring *K*_{leaf}, 30 min was chosen as a minimum to ensure that leaves had acclimated to high irradiance and stomatal conductance had stabilized. Previous studies found these criteria to be sufficient for the stabilization of transpiration rate, water potential, and *K*_{leaf}. Tests for any change in transpiration rate, leaf water potential, and *K*_{leaf} with measurement time (after stable flow was established) across leaves of a given species for seven species with a wide range of leaf capacitance showed no relationship of *K*_{leaf} to measurement time (Scoffoni et al., 2008; Pasquet-Kok et al., 2010).

### Measurement of Leaf Anatomical Traits

We used measurements of 34 leaf anatomical traits (Table IV) across 14 species as described by John et al. (2013), based on light micrographs of fully hydrated leaves fixed in formalin-acetic acid, embedded in LR White, cut in transverse 1-µm sections using glass knives in a microtome, and imaged using a 20× or 40× objective.

### Outline of the Modeling Approach

We created a model that uses anatomical measurements to calculate the hydraulic conductances of pathways outside the xylem in leaves. The model is adapted from the framework developed by Buckley (2015), which calculates horizontal and vertical components of hydraulic conductance in each of three tissue types distal to the BS (epidermis, palisade mesophyll, and spongy mesophyll) and in each of three transport modes (apoplastic, transmembrane, and gas phase). We extended the framework to include the BS itself and the BSEs, applied it to a spatially explicit grid representing a single areole to compute the distribution of water potential across the areole, and used that distribution to compute total *K*_{ox} and its *K*_{b} and *K*_{ob} components.

The original framework of Buckley (2015) included a term for hydraulic resistance due to diffusion across the interior of each cell in series with the transmembrane resistance. Discussions with colleagues led us to recognize that water movement across the cellular interior may occur by bulk flow rather than by diffusion and that the resulting transcellular bulk flow resistance would be negligible relative to the transmembrane resistance. We thus omitted the transcellular resistance from MOFLO. This is discussed further in “Discussion.” We also assumed that the quantitative contribution of plasmodesmatal flow to transpired water movement is negligible, consistent with its narrow circular slit (of width 1–2 nm) available for water flow between the membrane at its perimeter and the interior desmotubule of the endoplasmic reticulum (Doelger et al., 2014).

### The Areole Grid

We simulated a transverse section through a circular areole (the smallest region of a leaf bounded by minor veins) as a grid. Therefore, our results apply to regions of the leaf bounded only by minor veins and not by the lower order (major) veins; although our model does not directly account for free-ending veinlets, the values of VLA used to estimate areole dimensions did include veinlets. This grid had 744 nodes: 24 horizontal (parallel to the epidermis) and 31 vertical (Fig. 14). The aspect ratio of 24:31 was based on the average ratio of areole radius to leaf thickness across species (0.77 ± 0.07; mean ± se). Each node represents a band of tissue delimited by outer and inner radii (horizontal distances from the areole center) and upper and lower depths (vertical distances from the upper leaf surface; Fig. 14). Representing circular bands of tissue as single nodes is equivalent to assuming that the areole is radially symmetrical. Areole radius was computed from VLA following previous models that considered the vein system as a square grid with unit edge length *x*; this implies that each areole is uniquely associated with a vein length of 2*x* and an area of *x*^{2}, so VLA = 2*x*/*x*^{2} = 2/*x* (Cochard et al., 2004; Sack et al., 2004). Equating this area with that of a circle of radius, *r*_{areole} (π × *r*^{2}_{areole} = *x*^{2}), gives *r*_{areole} = *x*/π^{0.5} = 2/(VLA × π^{0.5}).

Each tissue band (node) in the grid was identified with a tissue type (BS, upper or lower BSE, upper or lower epidermis, or palisade or spongy mesophyll). All bands in the top and bottom rows of the grid were identified as upper and lower epidermis, respectively, and all bands in the left-most column (which corresponds to the outer margin of the areole, aligned with the nearest minor vein) were identified as BS or either BSE (in heterobaric species) or mesophyll (in homobaric species). All other tissue bands were identified as either spongy or palisade mesophyll based on measured anatomical proportions (Fig. 14). Formulas for tissue identity at each band are given in Supplemental Text S1.

The heights of the top- and bottom-most rows were taken as the measured thicknesses of the upper and lower epidermis, respectively; the height of each of the remaining 29 rows was set as 1/29th of the remaining leaf thickness. For homobaric species, which lack BSEs, all column widths were set at 1/24th of the areole radius. For heterobaric species, which possess BSEs, the width of the outermost (left-hand) column was set equal to one-half of the measured BSE width (the other half of the BSE width would be associated with the next areole to the left), and the widths of all other columns were set at 1/23rd of the remainder of areole radius. The resulting differences in tissue band dimensions among columns and rows were taken into account when computing the cross-sectional areas and flow pathlengths for connections between adjacent nodes; calculations involving BS nodes were further modified to account for the mapping of the elliptical cross section of the BS onto a rectangular column of nodes (for details, see Supplemental Text S1).

### Computing Flows and Water Potentials in the Grid

We computed the steady-state distribution of water potential across the grid on the basis of mass conservation. For each node *i*, an expression for mass balance can be written as a linear function of the water potentials of all nodes, in which the coefficients are hydraulic conductances between adjacent nodes. For example, the sum of all flows into node *i* from adjacent nodes must equal the net flow out of node *i* through stomatal transpiration:(2)where ψ_{j} is the water potential at node *j*, *K*_{ji} is the conductance (mol s^{−1} MPa^{−1}) between nodes *i* and *j*, *E*_{i} is any loss of water from node *i* by stomatal transpiration (mol s^{−1}), and the sums are taken over all nodes in the grid (for nodes that are not directly connected to node *i*, the conductance *K*_{ji} will be 0). Water enters the grid from the xylem, which is treated as a reference node with a water potential of 0. This reference node is not part of the grid, but its existence and location are implicitly incorporated by including a term for xylem-to-BS hydraulic conductance (*K*_{xb}) in the equation for each BS node:(3)where the subscript b denotes a BS node. In the presence of vertical temperature gradients within the mesophyll, the conductances for vertical connections between mesophyll nodes will include both an anisothermal gas-phase component (*K*_{aniso,ji}), which depends on the temperature difference between the two nodes, and an isothermal component (*K*_{iso},_{ji}), which does not. Rewriting Equation 2 to separate these components gives:(4)*K*_{aniso},_{ji} is given by Equation 5 (which is based on Equation 15 in Buckley, 2015):(5)where *D*_{wa} is the molecular diffusivity of water vapor in air, *v*_{w} is the molar volume of liquid water, *p*_{sat} and *T* are the saturation vapor pressure and absolute temperature, respectively (at node *i* or *j* as indicated by subscripts), *R*_{gas} is the gas constant, *a*_{ji} and *l*_{ji} are the area and pathlength for the connection between nodes *j* and *i*, and γ_{ji} and β_{ji} are unitless corrections that convert simple areas and pathlengths, respectively, to those actually experienced by moving water (for details, see “Calculating the Conductance Matrix” below). The quantity *v*_{w}ψ_{i}/*R*_{gas}*T*_{j} on the right-hand side is <<1 for typical leaf water potentials (e.g. this quantity is 0.0145 for ψ_{i} = −2 MPa and *T* = 298 K), so it can be omitted with negligible error. Omitting that quantity from Equation 5 and multiplying both sides by ψ_{j} − ψ_{i} gives the anisothermal component of the vapor flow (mol s^{−1}) from node *j* to node *i*, *F*_{aniso,ji}, as:(6)Note that *F*_{aniso,ji} is identical to the term inside the second sum on the left-hand side of Equation 4. Because *F*_{aniso,ji} does not directly depend on water potentials, it can be moved to the right-hand side and combined with the stomatal transpiration flux to give a single term on the right-hand side of each linear equation:(7)Equation 7 represents a system of linear equations that can be expressed more compactly in matrix form, as the product of a square matrix of conductance coefficients (**K**) whose elements are the *K*_{iso,ji}, and a vector (**δψ**), whose elements are the water potentials at each node, expressed relative to xylem water potential (i.e. the steady-state water potential drawdowns from the xylem to each node), with a vector **e** comprising the *E*_{i} on the right-hand side:(8)This system can be solved for **δψ** by multiplying the inverse of **K** by the vector **e**:(9)We generated the vector of transpiration rates (the components *E*_{i} of the vector **e**) by multiplying a fixed and arbitrary transpiration rate per unit of leaf area (*E*_{leaf} = 0.005 mol m^{−2} s^{−1}) by the projected leaf area corresponding to each node at each transpiring leaf surface. For amphistomatous species (*Helianthus annuus* and *Romneya coulterii*), we partitioned total transpiration rate between the upper and lower leaf surfaces using the ratio of maximum stomatal conductances at each surface (estimated as the ratio of the products of mean stomatal density and mean inner pore length for each surface). We measured stomatal density by counting stomata in each of three 400× fields of view in three leaves per surface, per species, and measured pore lengths for four stomata in each field of view using ImageJ software. We thus estimated that 58.2% and 43.6% of transpiration occurred from the lower surfaces of *H. annuus* and *R. coulterii*, respectively. All other species were hypostomatous, so we assumed that all transpiration occurred from the lower surface.

### Calculating the Conductance Matrix

We generated the conductance matrix (**K**) as follows. First, we computed a set of intrinsic hydraulic conductivities, κ (molar flow rates [mol s^{−1}] per unit of water potential gradient [MPa m^{−1}] per unit of area [m^{2}]), for each transport mode (apoplastic, transmembrane, and gas phase). The spatial dimensions in these conductivities represent the actual pathlengths and actual flow areas experienced by water moving in a particular tissue. Those pathlengths and areas often differ from the simple or bulk values that one would infer from bulk tissue geometry (e.g. the apoplastic pathlength around a cylindrical cell is longer than the simple distance across that cell, and the area available for gas-phase flow is smaller than the total cross-sectional area). The second step, therefore, was to compute correction factors for pathlength and area in each tissue type and flow direction. The area correction was the ratio of actual flow area to simple (bulk) flow area (γ), and the pathlength correction was the ratio of actual flow pathlength to simple (direct) flow pathlength (β). Third, for each transport mode in a given tissue and flow direction, we multiplied κ by γ and divided it by β to give the corresponding bulk conductivity, *k*:(10)Fourth, we summed these bulk conductivities across transport modes for each tissue type and flow direction. Finally, for each connection between a pair of nodes (*j* and *i*), we converted the appropriate total bulk conductivity to a conductance (*K*_{ji}, flow per unit of water potential difference [mol s^{−1} MPa^{−1}]) by multiplying it by the bulk flow area (*a*_{ji}) and dividing it by the direct flow pathlength (*l*_{ji}) appropriate to the connection between those nodes:(11)For connections between different tissue types (with bulk conductivities *k*_{1} and *k*_{2}, say), we computed the total conductivity as (0.5/*k*_{1} + 0.5/*k*_{2})^{−1}. The *K*_{ji} comprise the elements of the conductance matrix **K** (denoted as *K*_{iso,ji} in Eq. 7). We derive the expressions for κ in the following section. Expressions for γ, β, *a*, and *l* are derived in Supplemental Text S1, and tabulated results for γ and β are given in Supplemental Tables S1 and S2, respectively.

### Calculating κ

We derived intrinsic conductivities from expressions given by Buckley (2015). Note that the term conductivity in that article referred to flow per unit of area per unit of water potential difference (mol s^{−1} m^{−2} MPa^{−1}), whereas here, we use the term conductivity to describe a flow per unit of area per unit of water potential gradient (mol s^{−1} m^{−2} [MPa m^{−1}]^{−1} = mol s^{−1} m^{−1} MPa^{−1}). Thus, in this article, conductances are computed by multiplying conductivities by flow areas and dividing them by flow pathlengths, as described earlier.

For diffusion across a single membrane, the flow per unit of area per unit of water potential difference is *P*_{m}/*R*_{gas}*T* (compare with Equation 1 in Buckley, 2015), where *P*_{m} is the osmotic water permeability of the membrane (m s^{−1}), *R*_{gas} is the gas constant (J mol^{−1} K^{−1} = Pa m^{3} mol^{−1} K^{−1}), and *T* is the absolute temperature (K). To convert this to an intrinsic conductivity, it must be multiplied by one-half of the transcellular pathlength, *L*_{c} (m; because two membranes are encountered for every bulk distance *L*_{c} traveled, the value of *L*_{c} differs among tissue types and flow directions). Thus, the intrinsic conductivity for transmembrane pathways is:(12)The intrinsic conductivity for free diffusion of water, other than across membranes, is:(13)where *D*_{ww} is the molecular diffusivity for water in liquid water (m^{2} s^{−1}). The intrinsic conductivity for bulk flow of water through cell walls with nanopores having an effective Poiseuille radius of *R*_{a} is:(14)where η is the dynamic viscosity of water and *v*_{w} is the molar volume of liquid water. (Note that the analogous expression in Buckley [2015], in Equation 8, also contains factors that appear in our area and pathlength correction factors [γ and β], which are derived in Supplemental Text S1.) For gas-phase transport (water vapor diffusion), the intrinsic conductivity (κ_{gas}) contains an isothermal term that does not depend explicitly on vertical temperature gradients in the leaf and an anisothermal term that does depend on such gradients. The isothermal term is:(15)where *D*_{wa} is the molecular diffusivity of water vapor in air and *p*_{sat} is the saturation vapor pressure. The anisothermal term is:(16)where the subscripts j and i refer to values at the nodes above and below the internodal connection for which κ_{gas,aniso} is to be calculated. Equation 16 requires the vertical distribution of temperature to be specified. We assumed that temperature varied parabolically with depth in the leaf, relative to a maximum temperature in leaf (*T*_{max}) at a relative depth in leaf at which temperature is greatest (*z*_{max}), such that the temperature drop from the maximum value to the lower surface was equal to an input parameter, Δ*T*. Thus:(17)where *z* is relative depth (*z* = 0 and 1 at the upper and lower leaf surfaces, respectively). For each species, we set Δ*T* proportional to leaf thickness, such that its default value was 0.1°C for the mean leaf thickness of 292.5 μm. We assumed *z*_{max} = 0.25, based on Rockwell et al. (2014).

### Computing Integrated Leaf-Level Hydraulic Conductances

In experimental studies, *K*_{ox} is typically calculated from *K*_{leaf} and *K*_{x} using Equation 1, and *K*_{leaf} in our study was determined using the evaporative flux method described above. We note that there are a number of other methods in use for determining *K*_{leaf}, such as the high-pressure flow method (Yang and Tyree, 1994), the rehydration kinetics method (Brodribb and Holbrook, 2003), and the vacuum pump method (Martre et al., 2001; for review of methods and their contrasting assumptions and difference in simulated flow pathways, see Sack and Tyree, 2005; Flexas et al., 2013). Although several studies have shown that the different methods tend to yield similar maximum *K*_{leaf} values (Sack et al., 2002; Scoffoni et al., 2008), we highly recommend the use of the evaporative flux method for the most accurate representation of outside-xylem hydraulic pathways, since water movement in this method would most closely resemble that of a naturally transpiring leaf. In the evaporative flux method, *K*_{leaf} is defined as the ratio of *E*_{leaf} to the difference between stem water potential (ψ) and bulk leaf ψ (ψ_{leaf}) of a leaf bagged during transpiration and then equilibrated. Generally, the equilibrated ψ_{leaf} is assumed to represent the volume-weighted average over the mesophyll cells in the transpiring leaf. This assumes that negligible water is taken up from the xylem to the mesophyll during equilibration, which would be the case if the open conduits in the petiole of the excised leaf contained negligible volume, an assumption that requires testing, given that many species have open vessels of several centimeters extending from the petiole into higher vein orders (Tyree and Cochard, 2003; Chatelet et al., 2011; Scoffoni and Sack, 2015). Even accepting this typical assumption, additional ambiguity in the partitioning of *K*_{ox} into *K*_{b} and *K*_{ob} components arises when *K*_{ox} is calculated using the bulk water potential of the entire symplast. As we show in Supplemental Text S1, this can lead to spurious differences in *K*_{ob} between leaves even when those leaves have identical flow properties outside the BS. These artifacts can be traced to the fact that the bulk water potential used to compute *K*_{ox} includes some tissues that are proximal to the transport pathways that *K*_{ob} is meant to represent.

To allow simulated values of *K*_{ox}, *K*_{b}, and *K*_{ob} to be interpreted as independent measures of outside-xylem, across-BS, and outside-BS hydraulic conductances, respectively, we defined these conductances, for modeling purposes, in terms of a water potential gradient whose end point is distal to the BS. Specifically, for modeling purposes, we defined *K*_{ox} as:(18)where δψ_{ob} is the volume-weighted average water potential drawdown from the end of the xylem to all tissues distal to the BS, given by Equation 19:(19)where *i* is an index representing all non-BS nodes in the grid, *v*_{i} is the volume of liquid water in the tissue band represented by node *i*, and δψ_{i} is the component of **δψ** for node *i* (calculation of *v*_{i} for each node in the grid is described in Supplemental Text S1). We defined *K*_{b} as:(20)where δψ_{bn} is the volume-weighted average water potential of all nodes immediately adjacent to (and distal to) the BS (where bn stands for BS neighbors). Finally, we defined *K*_{ob} as:(21)To allow direct comparison between measured values of *K*_{ox} (defined by Eq. 1) and modeled values (computed by Eq. 18), we also computed alternative modeled values of *K*_{ox} based on the volume-weighted average water potential of all tissues distal to the xylem:(22)where δψ_{ox} is computed in the same fashion as δψ_{ob} but extended to include the BS itself. Modeled *K*_{ox} values from Equation 18 are given in most cases in “Results”; values from Equation 22 are used only when being compared directly with measured values (Table II; Fig. 1).

### Unknown Parameters

MOFLO contains six parameters that could not be estimated with the same confidence as other anatomical parameters (Table III). These are: (1) the percentage suppression of BS apoplastic transport by a suberized layer in BS cell walls; (2) Δ*T*; (3) *R*_{a}; (4) *P*_{m}; (5) ρ_{ta} (discussed further below); and (6) ρ_{fcph}. For the percentage suppression of BS apoplastic transport, we explored the full range of possible values (from 0%–100%); we set the default value at 0% because there is little evidence of BS suberization in leaves of most species (Lersten, 1997). Because measured *f*_{cph} is most likely an overestimate (light micrographs typically cannot distinguish true horizontal connections between palisade cells and the illusion of connections created by overlap of cells in the depth plane), we set the default value for ρ_{fcph} at 0 and explored a range from 0 to 1. We used a default value of 0.1°C and a range from 0°C to 0.2°C for Δ*T*, which spans the range of values in simulations by Rockwell et al. (2014). (Note that 0.1°C was the average temperature drop from the point of maximum temperature to the lower epidermis across species; in practice, we used the same gradient for all species, so that the absolute drop varied across species in relation to leaf thickness, as discussed earlier below Eq. 17.) We explored values of *R*_{a} from 0 to 10 nm and values of *P*_{m} from 0 to 160 μm s^{−1}, with default values of 3 nm and 40 μm s^{−1}, respectively, based on Buckley (2015).

Our anatomical measurements (John et al., 2013) suggest that cell walls in our species range from 0.5 to 2.9 μm in thickness, averaging 1.4 μm across tissue types and species. These values are about 5 times greater than published measurements made in other species based on transmission electron microscopy (Evans et al., 1994; Moghaddam and Wilman, 1998; Hanba et al., 2002; Scafaro et al., 2011). LM measurements of cell wall thickness might be affected by optical artifacts (blurring near the limit of optical resolution might increase apparent wall thickness) or sampling artifacts (e.g. if a cell’s perimeter is oblique to the sectioning plane, say at an angle of 45°, then the perimeter will appear at least 0.7 μm thick in a 1-μm section [approximately 1/tan(45°)] regardless of true cell wall thickness). On the other hand, fixation for transmission electron microscopy requires strong dehydration that may cause cell wall shrinkage. Accurate measurement of cell wall thickness is a future research direction; for this study, we assumed by default that LM measurements were overestimates by a factor of 5, so the default value of ρ_{ta} was 0.2 and we explored a range from 0.2 to 1.

### Simulations to Determine the Effects of Parameters on *K*_{ox}

MOFLO contains three classes of parameters: eight known biophysical parameters such as molecular diffusivity and dynamic viscosity (Table I), six unknown parameters, discussed above, that were either ambiguous in light micrographs or could not be estimated visually (Table III), and 34 known parameters that were confidently estimated from light micrographs of transverse leaf sections (Table IV). We performed three sets of simulations to explore the effects of these 48 parameters on *K*_{ox}.

#### Simulation Set 1

To bound the range of possible *K*_{ox} values consistent with the model, we varied the six unknown parameters simultaneously in two simulations: one with values chosen to minimize *K*_{ox} and another with values chosen to maximize *K*_{ox}. The low-*K*_{ox} values were *R*_{a}, *P*_{m}, and ρ_{ta} = 50% of their respective default values, Δ*T* = 0, ρ_{fcph} = 0, and 100% of BS apoplastic transport blocked by a Casparian strip. The high-*K*_{ox} values were *R*_{a}, *P*_{m}, and ρ_{ta} = 150% of their respective default values, Δ*T* = 0.20°C, ρ_{fcph} = 1, and no BS Casparian strip.

#### Simulation Set 2

To determine the mechanistic effect of each known parameter on *K*_{ox}, *K*_{b}, and *K*_{ob}, and to distinguish between these mechanistic effects and the across-species correlations between each parameter and *K*_{ox}, we varied each of the 34 known parameters across the range of measured values across species, while holding the other 33 known parameters at their all-species means and holding the six unknown parameters at default values. These parameter ranges, means, and default values are given in Tables III and IV.

#### Simulation Set 3

To explore the importance of uncertainty in the unknown parameters for conclusions drawn from the other simulations, we varied each of the six unknown parameters across a range of five likely values, while holding the other five unknown parameters constant at their default values and holding the 34 known parameters at their measured values for each species; these simulations were repeated for each species. Embedded in these 420 simulations (five values × six unknown parameters × 14 species) are 14 (one per species) in which each unknown parameter was at its default value; these 14 simulations give the default predictions for each species.

### Supplemental Data

The following supplemental materials are available.

**Supplemental Table S1.**Formulas for ratios of actual flow area to bulk flow area used to calculate bulk conductivities from intrinsic conductivities.**Supplemental Table S2.**Formulas for ratios of actual flow pathlengths to direct flow pathlengths used to calculate bulk conductivities from intrinsic conductivities.**Supplemental Text S1.**Derivation of expressions for flow area, pathlength corrections, grid areas, grid pathlengths, and volume averaging basis adjustment.

## Acknowledgments

We thank Helen Bramley, Antonio Diaz-Espejo, John Evans, and Margaret Barbour for helpful discussions.

## Footnotes

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Thomas N. Buckley (t.buckley{at}sydney.edu.au).

T.N.B. and L.S. conceived the project; T.N.B. designed, created, and operated the model and drafted the article; G.P.J. performed the leaf anatomy measurements; C.S. performed the hydraulic conductance measurements; all authors contributed to the writing and editing of the article.

↵1 This work was supported by the U.S. National Science Foundation (grant no. 1146514), the Australian Research Council (grant nos. DP150103863 and LP130101183 to T.N.B.), the Bushfire and Natural Hazards Cooperative Research Centre, and the Grains Research and Development Corporation.

↵[OPEN] Articles can be viewed without a subscription.

### Glossary

*K*_{leaf}- leaf hydraulic conductance
*R*_{ox}- resistance outside the leaf xylem
*K*_{ox}- outside-xylem hydraulic conductance
- VLA
- vein length per unit of leaf area
*K*_{x}- leaf xylem hydraulic conductance
- BS
- bundle sheath
- BSE
- bundle sheath extension
- MOFLO
- mesophyll and outside-xylem flow
*K*_{b}- bundle sheath conductance
*K*_{ob}- outside-bundle sheath conductance
- δψ
- water potential drawdown outside the xylem
*R*_{a}- Poiseuille radius of apoplastic nanopores
*P*_{m}- osmotic water permeability of cell membranes
*f*_{cph}- fraction of horizontal palisade surface area in contact with adjacent palisade cells
- ρ
_{fcph} - ratio of true palisade horizontal connectivity to value from microscopy
- ρ
_{ta} - ratio of true cell wall thicknesses to values from microscopy
*f*_{cs}- spongy mesophyll connectivity
*p*_{s}- spongy mesophyll airspace fraction
- LM
- light microscopy
*r*_{s}- spongy mesophyll cell radius
- κ
- intrinsic hydraulic conductivities
- Δ
*T* - vertical temperature gradient within leaves

- Received May 20, 2015.
- Accepted June 15, 2015.
- Published June 17, 2015.