Linking functional traits to multiscale statistics of leaf venation networks

Leaf venation networks evolved along several functional axes, including resource transport, damage resistance, mechanical strength, and construction cost. Because functions may depend on architectural features at different scales, network architecture may vary across spatial scales to satisfy functional tradeoffs. We develop a framework for quantifying network architecture with multiscale statistics describing elongation ratios, circularity ratios, vein density, and minimum spanning tree ratios. We quantify vein networks for leaves of 260 southeast Asian tree species in samples of up to 2 cm, pairing multiscale statistics with traits representing axes of resource transport, damage resistance, mechanical strength, and cost. We show that these multiscale statistics clearly differentiate species’ architecture and delineate a phenotype space that shifts at larger scales; functional linkages vary with scale and are weak, with vein density, minimum spanning tree ratio, and circularity ratio linked to mechanical strength (measured by force to punch) and elongation ratio and circularity ratio linked to damage resistance (measured by tannins); and phylogenetic conservatism of network architecture is low but scale-dependent. This work provides tools to quantify the function and evolution of venation networks. Future studies including primary and secondary veins may uncover additional insights.


Introduction
Leaves have venation networks with architecture that varies widely, from a single vascular strand (e.g. pines) to purely branching structures (e.g. Ginkgo) to open net patterns (e.g. many ferns) to mostly parallel structures in monocots, to highly reticulate patterns in many angiosperms. Networks vary over multiple spatial scales, with several levels of branching at length scales from 10 −5 m (the radius of a single vein) to 10 0 m (the length of some large leaves) (Roth-Nebelsick et al., 2001;. We define spatial scale as a characteristic feature of the network with a certain extent, for example the area of a vein loop or the radius of a vein. Network architecture may be closely linked to multiple functions (Ronellenfitsch & Katifori, 2019). Selection could act to maximize efficiency of resource transport, resistance to damage, mechanical strength, or to minimize total construction cost (Blonder et al., 2018). For resource transport, veins are implicated in photosynthesis and transpiration, as vein-mediated resource transport comprises a large portion of total leaf conductance in most species (Brodribb et al., 2007), though variation in other tissues (e.g. bundle sheath extensions) may also be important (Ohtsuka et al., 2018). For resistance to damage, the presence of loops in the network (Katifori et al., 2010) could prevent the propagation of embolisms that reduce conductance under low water potential (Brodribb et al., 2016), as well as the propagation of tears or cracks (Vincent, 1982;Niklas, 1999) caused by wind or herbivores. Additionally, defensive secondary compounds (e.g. latex) transported by the network could promote further resistance against herbivory (Agrawal & Konno, 2009), as damaged leaves with redundant flow pathways could still deliver latex or other signaling molecules to other portions of the leaf. For mechanical strength, the network could provide a skeleton that allows leaves to remain upright in wind, ultimately supporting photosynthesis (Givnish, 1979;Givnish et al., 2005). For construction cost, lignified tissue comprising veins is costly to construct relative to other tissues (John et al., 2017), and may also displace leaf volume that could be allocated to photosynthetic functioning. Thus, the realized network architecture may reflect tradeoffs among these different functions.
These functional tradeoffs and network architecture may both contribute to function at a range of spatial scales within a leaf. Thus, network features of different spatial scales may make varying contributions to functioning. For example, minor veins contribute more to conductance and photosynthesis due to their close spacing and disproportionate impact on hydraulic resistance (Brodribb et al., 2007;McKown et al., 2010). Conversely, major veins contribute more to cost (John et al., 2017) and mechanical strength (Hua et al., 2019) due to their disproportionate mass allocation.
From an evolutionary perspective, some aspects of network architecture may be more labile than others. Larger-scale veins could be phylogenetically conserved over deeper time (Ellis et al., 2009). Smaller-scale vein patterns are thought to evolve more quickly and show developmental plasticity .
At a single scale (e.g. of minor veins), network architecture can be described using a range of geometrical or topological statistics (Table 1). Single-scale studies of network architecture and leaf function have sometimes produced clear functional linkages (e.g. between vein density and resource transport, or between vein branching and defense chemistry; Blonder et al., 2011Blonder et al., , 2016Blonder et al., , 2017Hua et al., 2019) but in other contexts have not (e.g. between vein density and traits associated with cost, like leaf mass per area; Li et al., 2015). Recent models linking network architecture to leaf function (Blonder et al., 2011) have neglected this multiscale variation.
If network function varies across spatial scales, then network architecture should also be described across spatial scales, to better understand the diversity of evolved forms and the rules underlying this multiscale architecture. While such a perspective is sometimes implicit (Sack et al., 2012;Hua et al., 2019;Kawai & Okada, 2019), prior studies have used categorical descriptors to characterize features at other spatial scales (Sack & Frole, 2006;Sack et al., 2008;Ellis et al., 2009).
There is a need for multiscale description of venation network architecture, as well as for comprehensive assessment of its implications for multiple leaf functions. Recent efforts have focused on the idea of hierarchical loop decomposition (HLD) (Katifori & Magnasco, 2012;Mileyko et al., 2012;Ronellenfitsch et al., 2015). In this approach, a leaf venation network is considered as comprising a set of loops that are nested within each other. The leaf can be partially summarized via a hierarchical tree describing how smaller loops are nested within larger loops, and the statistical properties of the tree (e.g. nesting ratio, topological length) can be used to characterize variation between species. This approach provides an advance in terms of explicitly considering architecture across spatial scales. However, empirical applications of this framework have been constrained by limited venation data and by the unclear linkage between the hierarchical tree properties and hypotheses about leaf functioning (Ronellenfitsch et al., 2015).

Framework for multiscale venation statistics
Here we present a framework, built on HLD concepts, to define multiscale venation statistics. The goal of the framework is to describe scale transitions in network architecture statistics (Table  1), so that a set of single-scale statistics ω i f g, can be transformed into a set of continuous functions ω i r ð Þ f g, where r is a metric of spatial scale.
HLD analyses are based on 'spatial graph' extractions of leaf venation networks, that is, with all vein segments described in terms of their position, length, radius, and connectivity with other vein segments, and all areoles described in terms of their area and shape. The network extraction step can be carried out using software available in Xu et al. (2020). This information is necessary to determine which areoles are nested within which other areoles. Loops are defined either in terms of regions wholly enclosed by veins, or regions enclosed by veins and leaf boundary; regions not wholly enclosed by veins (e.g. due to cropping of an input image) are excluded from the analysis. Vein segments that constitute boundaries between loops are identified and described in terms of their radius (r). Then, a pair of loops with the minimum radius boundary segment (r min ) are fused together, and the boundary segment is deleted from the network. This process begins with fusing the areoles linked by the vein with the smallest r min , and is iterated until all loops are fused together into one large loop representing the entire leaf, with r min equal to the width of the largest vein segment ( Fig. 1). Only loops that are complete are considered (i.e. excluding those cut off by the image boundary). Large veins are split into multiple segments spanning branching points with other segments.
Single-scale statistics ω i f g are then calculated at each of the network fusion steps, where i indexes the various statistics. If a network has n boundary vein segments, then the hierarchical loop decomposition will yield n values of ω i f g, each at a different value of r. The outcome of this process is a scale-dependent description of the network that remains when considering only veins above a certain size class.
These ω i f g values can then be converted into multiscale statistics, ω i r ð Þ f g, which encapsulate the scale transitions that may occur in network architecture. Importantly, not all leaves will have veins present at all scales. Thus, to obtain scale-dependent functions, the first step is to construct a hurdle function describing the presence or absence of a vein at a given scale, gis constructed describing the statistic of interest, for example VD(r), which is the vein density across scales. The final function is: Estimation proceeds by binning the data into uniform ranges of r, then calculating bin-median values of η i and σ i . Using raw data leads to undesirable sample sizenetwork scale correlations, while fitting functional forms to the curve via regression is not Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust New Phytologist (2020) 228: 1796-1810 www.newphytologist.com necessarily appropriate because the ω i values are piecewise continuous rather than fully continuous (i.e. they include gaps where no features of a given size class exist). Vein absences can arise either because a vein is truly absent at a scale, or because the sampled area was too small to be representative (i.e. undersampling).
If the input network is complete (i.e. comprising a whole leaf) then the absences are real, because the input data are a statistical population rather than a statistical sample. In leaf subsections, however, undersampling biases may be important. In these cases, one can assume that network architecture falls on a scale continuum over the scale range of the data, that is, η i r ð Þ ¼ 1 for all r, and then σ i r ð Þ can be gap-filled using interpolation. Sampling limitations can also increase the frequency of gaps and/or bias estimates. The sampling uncertainty in ω i

Research
New Phytologist necessarily increases with increases in r. This is because the number of network components must decrease as portions of the network are deleted. Therefore, analyses should be restricted to ranges of r where more than a certain fraction of samples have η i r ð Þ ¼ 1. Once the multiscale statistics are estimated, they can be used to describe variation among leaves and species, and to serve as potential predictors of leaf function. To proceed, one can use the values directly, or can reduce them to scalars like the linear slope estimate of ω i r ð Þ vs r. When using the predicted values of ω i r ð Þ f g in regression analyses, it is necessary to account for the nonindependence of values at similar spatial scales, for example via partial least squares (PLS) methods.

Hypotheses about trade-offs between network architecture and function
We focused on traits related to four multiscale statistics: vein density (VD), the mean elongation ratio of loops (ER), loop circularity (CR), and the minimum spanning tree ratio (MST) ( Table 1). We also focused on linking them to four leaf functions (resource transport, damage resistance, mechanical strength, and construction cost) ( Table 2).
Vein density, ER, and MST have been shown to capture leading axes of variation in network architecture across many species at single scales (Blonder et al., 2018), while CR(r) may reflect relative allocation of vein perimeter relative to loop area (Blonder et al., 2011). Based on the availability of trait data in this study, we represented resource transport by maximum photosynthetic capacity, though recognizing that this metric of carbon flux does not fully represent water fluxes (i.e. maximum hydraulic conductance; Brodribb et al., 2007), for which we did not have data. Damage resistance was represented by tannin fraction, as this type of secondary chemistry is commonly used as an herbivory defense (Coley & Barone, 1996). To represent mechanical strength we used force to punch, as this is a direct metric of strength (Pérez-Harguindeguy et al., 2013). Finally, we represented cost by leaf mass per area, as dry mass investment is strongly related to tissue construction cost (Poorter et al., 2009).
Specifically, we asked (1) whether multiscale statistics enable useful measurement of species network architecture, (2) whether network architectures vary systematically across spatial scales, (3) whether expected architecture-function trade-offs are supported across different spatial scales, and (4) whether the phylogenetic conservatism of architecture varies across spatial scales. Then, we hypothesized that: (1) The four statistics of network architecture will represent statistically independent axes of variation, such that a wide range of network architectures are possible.
(2) Different portions of the architectural space will be occupied at different scales. That is, network architecture will shift systematically as scale increases.
(3) Distinct trade-offs exist between network architecture and functioning across spatial scales (Table 2). We predict that resource transport will be positively linked to VD(r) at small values of r, due to the importance of small veins in determining hydraulic conductance (McKown et al., 2010), as well as positively linked to MST(r) across scales, as more tree-like networks distribute resources better (Katifori et al., 2010). Damage resistance will be negatively linked to MST(r) and CR(r) across scales, as high secondary chemistry investment may offset investment in redundant flow pathways. Mechanical strength will be positively linked to VD(r) at large scales, as large veins may contribute disproportionately to stiffness and because lignified vessel sclerenchyma are associated with fracture resistance (Choong et al., 1992) and negatively linked to ER(r) at small scales, as more parallel-veined laminas may be easier to tear. Cost will be positively linked to VD(r) at large scales, due to the importance of large veins (which are lignified and dense) in determining total leaf cost (John et al., 2017).
(4) Phylogenetic niche conservatism in all network multiscale statistics will be lowest at small values of r due to their presumed development and evolutionary lability (Trivett & Pigg, 1996;Roth-Nebelsick et al., 2001;Ellis et al., 2009).
We tested these hypotheses using a phylogenetically diverse set of 260 tree species in a forest in Malaysian Borneo. Our analysis was based on leaf samples of c. 2 cm 2 , which contained veins of several orders but typically did not include primary or secondary veins.

Materials and Methods
This study integrates venation network imagery (Blonder et al., 2019) with functional trait data that were collected from the study sites (Table 3) (Both et al., 2018).

Sites and sampling design
Samples were collected in Malaysian Borneo during July-December 2015, within eight 1-ha permanent forest plots comprising mixed dipterocarp lowland forest (Table 3). We pooled data from all sites, as this study focused on functional and evolutionary rather than environmental questions. Vouchers are stored at the Danum Valley Field Centre, Sabah, Malaysia.
Within each plot, all trees ≥ 10 cm in diameter at breast height (DBH) were identified and measured for DBH and height. Branches of 2-4 cm diameter were then collected using rope climbing or pole pruning techniques. Branches were collected from species with highest contribution to biomass, that is, species comprising the top 70% of plot basal area. Additional leaves of rarer taxa were also collected from all trees in three 20 × 20 m subplots within each plot. When possible, both sunlit and shaded branches were collected from each tree, resulting in sampling heights across the range of 2-53 m.

Venation networks
A single mature undamaged leaf was collected from each branch, cleaned with a wet rag, then pressed flat and dried at 60°C for several days. A 1-4 cm 2 section was then cut from each leaf midway from the base to the apex and midway from a primary vein to either another primary vein (for palmate leaves) or to the margin (for pinnate leaves). These dried sections were rehydrated to reverse the effects of sample shrinkage (Blonder et al., 2012) and chemically cleared and stained (Pérez-Harguindeguy et al., 2013;Blonder et al., 2018). Samples were mounted on glass slides and imaged using a compound microscope (Olympus BX43) with a ×2 apochromat objective and a color camera (3840 × 2748 pixel resolution; Olympus SC100; Southend-on-Sea, UK). Approximately 16 image fields were then stitched together to obtain a full image of each sample (resolution 595 px mm −1 ). Image contrast was enhanced using a contrast-limited adaptive histogram equalization on the green channel of each image. Images were then imported into GIMP (GNU Project) image-editing software and the boundaries of the sample (vs background) were manually delineated with a polygon. Next, a c. 700 × 700 px 2 region in each sample was manually segmented using a tracing tablet. The hand-traced vein images were then segmented into binary representations, in which vein pixels were given one value, nonvein pixels another value, and image background an NA value.
These validated data (comprising 686 881 432 manually segmented pixels) were used to train a machine learning algorithm, as described by Xu et al. (2020). An ensemble of six convolutional neural networks (CNNs) was developed, with each network implemented via a U-Net architecture (Ronneberger et al., 2015). To avoid overfitting, input image data were repeatedly rotated and scaled during training. Each CNN was trained on 5/6 of the input dataset and predicted using a sliding window approach on the remaining 1/6. Probability outputs were then averaged across all CNNs in the ensemble and converted to binary predictions using a threshold determined by Receiver Operating Characteristic (ROC) plots based on the prediction relative to the manually delineated validation region in each image. The algorithms sometimes under-segmented large veins that were not included in the validation dataset. Any veins with widths greater than c. 500 µm were also manually segmented throughout the entire sample. These segmentations were then added to the CNN-based segmentations. Segmentations were then masked to the manually delineated boundaries. A total of 32 815 701 653 pixels were algorithmically segmented using this methodology.
The segmented images were then filtered, so that images with sample areas < 90 mm 2 and Fbeta2 scores < 0.8 were discarded. Among the retained images, the segmentation quality was high (Fbeta2 mean = 0.95 AE 0.03 SD). A small number of images were assigned new branch codes different than those in the original image dataset, based on recent revision of handwritten sample labels (Supporting Information Table S1).
Spatial graph representations of the networks were extracted from each segmented image. Next, all unique loops and their boundary vein segments were identified. Unclosed loops (on the boundary of the sample) were removed from the analysis. The four network statistics (Table 1) then were obtained by following the method described by Xu et al. (2020).
After calculating these statistics, the network was iteratively pruned, and loops fused, by sequential removal of the smallest remaining boundary vein segment. After each fusion event, network statistics were re-calculated, along with the minimum vein radius (r min , µm). Multiscale statistics were truncated to the 0-0.2 mm range of r min to minimize undersampling biases (Fig.  S1), though veins of up to 0.58 mm were present in the data.

Trait measurements
Traits were measured from each branch using mature undamaged and cleaned leaves. Different leaves were used for each trait, then data were pooled at branch level. Measurements are described in Table 2 and in a study by Both et al. (2018).

Statistical analysis
We first constructed multiscale statistics for each leaf by binning data into 50 r min bins spanning 0-0.2 mm. As we worked with

Research
New Phytologist leaf subsections, the absence of veins at a certain scale could be caused either by true absences or by undersampling. In cases where η i r ð Þ ¼ 0 occurred, we used two approaches. In the first approach, we assumed undersampling, then set η i r ð Þ ¼ 1 and filled the missing values of σ i r ð Þ using linear interpolation. Extrapolation was handled by using the value at the closest data extreme. In the second approach, we assumed the data were accurate as observed, then separately modeled the η and σ components of ω. For all the main text analyses, except the test of Hypothesis 3, we only used the first approach, as the large number of incomplete cases prevented application of the necessary statistics. Analyses were conducted treating branches as independent and identically distributed random replicates.
To determine the independence of architectural axes across scales (Hypothesis 1), we carried out a principal component analysis across the four axes of VD, ER, CR and MST, using values for each leaf at each scale as replicates. We then assessed variance allocated to each principal component to determine the relative independence of axes.
To evaluate whether different portions of the architectural space are occupied at different scales (Hypothesis 2), we also visualized the principal component scores using convex hulls and 95% confidence ellipses at each scale.
To assess evidence for architecture-function trade-offs across scales (Hypothesis 3), we used two complementary approaches. In the first approach, we carried out partial least squares (PLS) regressions for each of the functional traits. Functional traits were treated as response variables and the multiscale statistics in each bin as predictors. PLS allows for incorporation of nonindependence of the statistics across scales. All four multiscale statistics were included as predictors. Each multiscale statistic was z-transformed before analysis to improve comparability of parameter estimates across multiscale statistics. For all models, most of the variation was explained by the first component, so we restricted analysis to this component. We reported the variance explained by the first component, and also assessed the scale-dependent effect of each venation trait at each scale as the value of the respective loading coefficient. Each hypothesis was then tested by examining plots of loading coefficients vs r min .
In the second approach, we used ordinary least squares regression models, one for each functional trait, multiscale statistic, and r min binned value combination. At each combination, we pooled data for all species with η i r ð Þ ¼ 1 (i.e. only those species with observed veins at that scale). We then conducted the regression using the functional trait values as the response variable and the σ i r ð Þ values as the predictor variable. We reported the overall variance explained by the regression and the slope estimates. This approach does not account for nonindependence of data across spatial scales, or nonindependence of predictors, nor does it allow us to assess statistical significance (due to the high falsediscovery rate). However, it does allow for use of more data in each analysis, as missing cases do not need to be dropped or interpolated.
To determine levels of phylogenetic niche conservatism across scales (Hypothesis 4), we used the binned values of each multiscale statistic. For each spatial scale and each multiscale statistic, we computed the K statistic (Blomberg et al., 2003). Values of K < 1 indicate faster trait evolution than expected under a Brownian model (trait lability), while values of K > 1 indicate slower trait evolution (phylogenetic niche conservatism). The tree was estimated based on two recently published mega-trees (Jin & Qian, 2019).
All statistical analyses were conducted in R software (v.3.5.1), using the V.PHYLOMAKER, APE and PLS packages.

Data availability
Venation data, subsetted trait data, and the phylogeny are available at https://doi.org/10.5287/bodleian:QR11d1PD2. R code to replicate all analyses using these files is available at https://doi. org/10.5287/bodleian:E9JP2gjyP. Venation network image data are available in Blonder et al. (2019). Segmented venation network images (like those shown in Fig. 2) are available in Xu et al. (2020). Algorithms to analyze images and to calculate multiscale statistics are available in Xu et al. (2020). The full trait dataset (including numerous other variables) is available at https://doi. org/10.5281/zenodo.3247631.

Quantification of species network architecture across scales
All four multiscale statistics demonstrated high variation among species, as well as complex patterns across spatial scales (Fig. 3), giving support to our Hypothesis 2. In general, relationships for each multiscale statistic were consistent among leaves within species, that is, intraspecific plasticity for network architecture traits was low, but variation between species was high. Data for focal species are shown in Fig. 3, with data for all species shown in Figs S2-S5. Some species maintained low values of ER (ER ≈ 1) across scales (Figs 3a, S2), reflecting circular loops nested within larger circular loops. Other species showed sharp increases (ER ≈ 5-10) at r min ≈ 0:1 mm, reflecting nesting of circular loops in longer loops (e.g. D. lanceolata). Most species then reduced ER at r min ≈ 0:2 mm, while a small number maintained high or increasing values (e.g. M. wrayi).
Many species had a high value of CR (CR ≈ 0:75) at small r min , then reached a lower value (CR ≈ 0:25) by r min = 0.2 mm, reflecting a transition from loops with more complex boundaries to loops with simpler boundaries (Figs 3b, S3). However, the rate of decrease varied strongly among species, with some showing Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust New Phytologist (2020) 228: 1796-1810 www.newphytologist.com leveling off at intermediate scales (P. pinnata) and others increases at larger scales (D. lanceolata). VD decreased monotonically with r min in all species, reflecting the greater prevalence of small veins relative to large veins in these species at the studied scales (Figs 3c, S4). However, the rates of decrease varied across species. Some species showed rapid linear decreases (H. crassifolia), consistent with vein tapering and a scale continuum, while other species showed abrupt decreases, consistent with discrete transitions between veins of different orders (M. wrayi and D. lanceolata).
Contrarily, MST increased monotonically with r min in all species (Figs 3d, S5), reflecting a tendency for more branching and less looping at larger spatial scales. At small r min , MST values ranged from 0.5-0.8, indicating that a large fraction of all vein segments contributed to loops rather than branches, while at larger r min , MST values approached 1 (i.e. no loops). Some species (H. crassifolia) had much more tree-like networks at all scales. Although MST values must reach 1 ultimately (i.e. when only a single vein remains in the network), it is mathematically possible for MST to both increase and decrease with r min at intermediate scales. However, no decreases occurred in this dataset. Abrupt increases at some scales were present in some species (P. pinnata), indicating discrete transitions in branching architecture.
The principal components analysis showed that the four multiscale statistics can be distinguished in at least three axes of variation (Fig. 4), partially supporting Hypothesis 1. The first axis represented high VD and low MST, and explained 67% of the variation in the data. The second axis represented high ER, and explained 23% of the variation in the data. A third axis represented high CR and explained 8% of the variation in the data. As r min increased, species tended to increase in PC1, decrease in PC2, and increase then decrease in PC3 (Fig. S6), indicating shifts at larger spatial scales from high to low VD, high to low MST, high to low ER, and low to high to low CR. Additionally, there was a pronounced shift from variation occurring primarily along PC1 and PC2 at small scales (r min < 0.05 mm), to variation occurring primarily along PC1 at intermediate scales (r min > 0.05 mm), and then primarily along PC2 at larger scales (r min > 0.1 mm).
Trade-offs in network architecturefunction across spatial scales PLS analysis showed that network architecture was more associated with mechanical strength (in the first PLS component, force to punch explained 10% of model variation) and damage resistance traits (tannins fraction: 8.9%) than to resource transport

Research
New Phytologist (photosynthetic capacity: −2%) and cost (leaf mass per area: 3.5%) (Fig. 5). Inclusion of additional PLS components did not qualitatively improve explained variation or reveal other network architecture linkages (data not shown).
The importance of each multiscale statistic varied across spatial scales. For force to punch (mechanical strength proxy), VD and MST made large contributions at r min < 0.05 mm, while ER and CR made large contributions at r min > 0.05 mm (Fig. 5a). These results are largely not consistent with the hypothesized positive linkage to VD at large scales and negative link to ER at small scales. For tannins (damage resistance proxy), only VD contributed at r min < 0.04 mm, while CR, ER, and MST contributed at r min = 0.05 mm, and only ER contributed at r min > 0.15 mm (Fig. 5d). These results are not consistent with the hypothesized negative linkage to MST and CR at all scales. Results for photosynthetic capacity (resource transport proxy; Fig. 5c) had slightly negative explained variation (i.e. the model was worse than one not including the predictors), which is also not consistent with the hypothesized positive linkage to VD at small scales and positive link to MST at all scales. Results for leaf mass per area (cost proxy) (Fig. 5b) also had low explained variation, which is not consistent with the hypothesized positive linkage to VD at large scales. Therefore, we did not find strong evidence for the predictions of Hypothesis 3.
We also repeated these analyses using single-variable regressions at each scale. Results were qualitatively similar (Fig. S7), indicating that the choice of methodology did not strongly influence our findings.

Phylogenetic conservatism of network architecture
All four multiscale statistics varied extensively across the phylogeny with some visual evidence of clustering at certain spatial scales (Fig. 6). Notably, clustering was apparent at small r min for VD and MST for the families Euphorbiaceae and Dipterocarpaceae, which comprise many species that are ecologically dominant in southeast Asia. By contrast, CR and ER varied extensively across the phylogeny at small r min . At larger r min , more variation was apparent, with low values of CR present in some clades (e.g. the Clusiaceae), and more rapid increases in MST in the Fagaceae and Dipterocarpaceae. ER showed extensive variation but no clear association with phylogeny at any scale. That is, patterns of phylogenetic clustering were highly scale dependent.
The K analysis was consistent with this perspective. All four multiscale statistics varied at least as rapidly (K < 1) as under a Brownian trait evolution model (Fig. 7). In the case of VD and CR, Brownian evolution (K = 1) could be rejected at almost all scales; for MST, K < 1 was supported at r min < 0.12 mm, and for ER, at some scales of r min < 0.08 mm. Contrary to Hypothesis 4, there was a trend for more trait lability (smaller K ) at larger r min .

Discussion
Broadly, we found that (1) the network architecture varied extensively among species and across scales, (2) the phenotype space shifted and decreased at larger scales, (3) the relationships between venation multiscale statistics and leaf function in this dataset were weak, and (4) the phylogenetic niche conservatism in multiscale statistics was low and variable across scales.

Independent axes of variation
Our framework provided a vocabulary for making quantitative and specific statements about network architecture variation. The four axes of VD, CR, ER, and MST provide complementary information about network architecture, with likely three independent axes present (Fig. 4). Other metrics of network architecture not considered here (e.g. vein branching angles; Hickey, 1973) may also be useful, so these four statistics should not be considered a complete set of descriptors.

Variation in architecture across scales
The architectural phenotype space varied across scales (Fig. 3). In particular, we found a reduction in the size of phenotypic space at larger scales, and a shift to more tree-like and less elongated network geometry. The overall shape of this phenotypic space is potentially consistent with the existence of a number of phenotypic extremes, in which intermediate phenotypes occur as tradeoffs among functions represented by the extreme functional archetypes (a 'Pareto front'; Shoval et al., 2012). Consistent with this perspective, theoretical work has identified extreme network phenotypes (e.g. all branching, all looping) corresponding to optimization of resource transport and resistance functions, respectively (Ronellenfitsch & Katifori, 2019). Our work now extends these models to a wider range of phenotypic axes and empirical datasets. However, the weak empirical results for Hypothesis 3 leave unclear whether the extreme phenotypes do indeed represent functional archetypes.
The phenotype space appeared to vary along certain axes differentially across scales, suggestive of multiple drivers affecting patterning at each scale. That is, variation (as seen via confidence ellipses and convex hulls) appeared primarily along PC2 at small and large scales, but along PC1 at intermediate scales. A simple Fig. 4 Principal components analysis of the four-dimensional architectural space comprising elongation ratio (ER), circularity ratio (CR), vein density (VD), and minimum spanning tree ratio (MST) across scales, at each value of minimum vein size (r min ). 95% confidence ellipses enclose the data at r min value to show the central tendency in the data; convex hulls are also shown to highlight species with extreme phenotypes. Parenthetical values indicate variance explained by each axis.

Research
New Phytologist null explanation is that reduced sampling at large scales reduces the available variation (e.g. MST must reach a value of 1 for all leaves once r min reaches the maximum vein size). However, this idea seems insufficient to explain the differential variation along other axes, as well as the contraction of confidence ellipses at intermediate scales.
These results highlight the complexity of the phenotype space, and provide new targets for models meant to capture their evolution (Ronellenfitsch & Katifori, 2017. Broadly, we found that network statistics at one scale are not predictive of statistics at another scale, and that the limits to phenotype space are scaledependent. This high heterogeneity also invites questions about the reasons why species may occupy different subsets of the architectural phenotype space, what factors drive shifts, and whether such variation leads to divergent or convergent functioning.

Functional linkages between network architecture and function
Each model explained < 10% of variation in each functional trait (Fig. 5). Stronger linkages were obtained for tannin fraction (damage resistance) and force to punch (mechanical strength), rather than LMA (cost) and photosynthetic capacity (resource transport). This finding is consistent with prior studies finding vein linkages to structural and defensive traits (Li et al., 2015;Blonder et al., 2018;Hua et al., 2019;Kawai & Okada, 2019) and inconsistent with those finding linkages to hydraulic or resource flux traits (Sack et al., 2008;Brodribb et al., 2016).
Although we found a linkage between force to punch and venation traits, the relationships observed were in the opposite directions of our predictions (Hypothesis 3). For example, we expected mechanical strength to exhibit a positive link to VD at larger scales and to ER at small scales, but we found the opposite pattern (Fig. 5a). The most likely explanation is that spatial scale variation in r min in these data was not enough to capture some of the 'large'-scale functional linkages. The largest veins studied here were 0.2 mm, when midveins can be at least 100× larger; similarly, the largest sample was c. 4 cm 2 in area, while leaves can be nearly 1 m 2 in area for some species (Kattge et al., 2011). These larger features could be important for explaining mechanical strength, and other axes of functioning. Therefore, describing network architecture across all orders of magnitude is a priority for future investigations. However, one cannot easily quantify large image extents at high resolution, at least with brightfield microscopy. One solution is to take a hybrid approach, using higher resolution microscopy for measurements at smaller scales and lower resolution cameras or scanners at larger scales. Because both analyses yield multiscale statistics with the same scale variable (r min ), it should be possible to fuse multiscale statistics from multiple ranges of scales, as shown in Fig. S8.  For damage resistance (tannin), we predicted a negative link to MST and CR at all scales, but found a significant contribution only at r min = 0.05 mm (Fig. 5d). This result suggests that a high investment in secondary chemistry may not be enough to offset investment in redundant flow pathways. However, we recognize that we examined damage resistance primarily through the lens of herbivory defense, but there are other aspects that may also be linked to venation. For example, the ability to resist damage from drought may also be linked to the venation network, as the probability of cavitation varies with vein size (Sack et al., 2008;Brodribb et al., 2016). Another important functional axis not explored in this study was resilience to damage. While resistance is related to the ability to prevent damage, resilience is the capacity to maintain function after damage has occurred. The presence of vein loops may provide redundant pathways that enable continued flow after damage, caused either by abiotic (e.g. drought, frost) or biotic factors (e.g. herbivory) (Sack et al., 2008;Katifori et al., 2010). However, too much redundancy may be detrimental and might actually decrease resilience, by facilitating the spread of embolisms or diseases.
This study assessed functional linkages of network architecture with a limited set of functional traits and functional axes. The underlying trait dataset we used here actually includes a larger number of traits, including secondary chemistry and isotope and elemental stoichiometry data. These other traits may map onto the functional axes we have already discussed. However, analyses using these additional traits as proxies for each functional axis did not qualitatively shift results, with maximum PLS1 R 2 for any trait of 11% (Fig. S9). However, other proxies might have relevant connections. In particular, the resource transport axis could be refined in future studies, separating out functions related to carbon uptake and water New Phytologist (2020) 228: 1796-1810 Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust www.newphytologist.com Research New Phytologist transport. For example, hydraulic conductance (unmeasured) might be more closely linked to water transport functions, whereas photosynthetic capacity would be more linked to carbon uptake functions (measured). The role of venation in predicting hydraulic conductance appears to be complex, even at single scales. Studies examining the linkage between minor veins and hydraulic conductance have obtained divergent results, suggesting that selection for this relationship may be clade-or biogeographically-dependent (Brodribb et al., 2007;Walls, 2011;Gleason et al., 2016). More broadly, various aspects of leaf functioning may be decoupled, if there is independent selection on multiple axes of functioning beyond the 'fast-slow' continuum (Wright et al., 2004;Reich, 2014), for example, as has been found in the decoupling between hydraulic and economic functioning in Chinese (Li et al., 2015) and Australian (Gleason et al., 2016) species. Such patterns might also occur if functioning is achieved through other traits not directly related to venation. For example, leaves may modify traits like their mesophyll density or their thickness in order to meet different structural constraints (John et al., 2017), and similarly, outside-xylem pathways (e.g. bundle sheath extensions) may be important factors in determining hydraulic conductivity (Buckley et al., 2015;Ohtsuka et al., 2018).
An alternate explanation for the weak architectural-functional relationships found in this study is that not all aspects of network architecture have immediate functional linkages. That is, much of the diversity of form seen in networks may reflect evolution in the absence of strong selection. For example, the weak relationships between leaf functional traits and ER (Fig. 5) may indicate that elongation of areoles may occur, but rather that it is a consequence than a cause of other evolutionary forces. For example, leaf aspect ratio and areole elongation are sometimes linked in nonmonocots (Blonder et al., 2016), and are potentially coupled via developmental mechanisms related to leaf elongation (Kang & Dengler, 2004). As such, selection on leaf aspect ratio could indirectly drive variation in ER.

Phylogenetic conservatism of network architecture
The results also showed that all venation network multiscale statistics varied more rapidly than expected under a Brownian evolution model (Blomberg's K < 1). Moreover, there was no increase in K at large scales for any statistic (Fig. 7). That is, the prediction of more niche conservatism at larger scales (Hypothesis 4) was apparently rejected. These results are consistent with prior analyses of fine-scale network architecture variation (Boyce et al., 2009;Blonder et al., 2018), and extend them to a wider range of scales. This finding runs contrary to ideas that large scale network architecture should be highly conserved, as is expected in plant systematics (Hickey, 1973), but is consistent with rapid rates of evolution seen in some adaptive radiations (Dunbar-Co et al., 2009;Blonder et al., 2016). However, these results should be interpreted with some caution as the sampling design was focused on dominant species; thus, patterns may change if rare species were included.
As in the functional analysis (Hypothesis 3), there was likely insufficient scale variation present in the image data to fully address this research question. Thus, our results do not yet call into question the utility of primary and secondary venation in www.newphytologist.com systematics. It seems likely that larger-scale analyses would likely shift the conclusions of this analysis, especially given the extensive prior use of primary and secondary venation characters as aids to species identification and systematics. Our framework provides an approach that, with additional data, will be able to delineate the scales and contexts in which these applications are wellfounded.

Conclusion
We advanced a framework to quantify network architecture across species and spatial scales. By building on ideas for hierarchical loop decomposition (Katifori & Magnasco, 2012), and leveraging recent machine learning algorithms (Xu et al., 2020), we were able to move past topological (Ronellenfitsch et al., 2015) or single scale (Blonder et al., 2018) descriptions of leaves toward multiscale statistics. Datasets with wider scale variation than ours will further clarify understandings of network form-function linkages.

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.         New Phytologist is an electronic (online-only) journal owned by the New Phytologist Foundation, a not-for-profit organization dedicated to the promotion of plant science, facilitating projects from symposia to free access for our Tansley reviews and Tansley insights.
Regular papers, Letters, Research reviews, Rapid reports and both Modelling/Theory and Methods papers are encouraged.
We are committed to rapid processing, from online submission through to publication 'as ready' via Early View -our average time to decision is <26 days. There are no page or colour charges and a PDF version will be provided for each article.
The journal is available online at Wiley Online Library. Visit www.newphytologist.com to search the articles and register for table of contents email alerts.
If you have any questions, do get in touch with Central Office (np-centraloffice@lancaster.ac.uk) or, if it is more convenient, our USA Office (np-usaoffice@lancaster.ac.uk) For submission instructions, subscription and all the latest information visit www.newphytologist.com New Phytologist (2020) 228: 1796-1810 Ó2020 The Authors New Phytologist Ó2020 New Phytologist Trust www.newphytologist.com