Adaptive trait syndromes along multiple economic spectra define cold and warm adapted ecotypes in a widely distributed foundation tree species

1. The coordination of traits from individual organs to whole plants is under strong selection because of environmental constraints on resource acquisition and use. However, the tight coordination of traits may provide underlying mechanisms of how locally adapted plant populations can become maladapted because of climate change. 2. To better understand local adaptation in intraspecific trait coordination, we studied trait variability in the widely distributed foundation tree species, Populus fremontii using a common garden near the mid-elevational point of this species distribution. We examined 28 traits encompassing four spectra: phenology, leaf economic spectrum (LES), whole-tree architecture (Corner's Rule), and wood economic spectrum (WES). 3. Based on adaptive syndrome theory, we hypothesized that trait expression would be coordinated among and within trait spectra, reflecting local adaptation to either exposure to freeze-thaw conditions in genotypes sourced from high-elevation populations or exposure to extreme thermal stress in genotypes sourced from low-elevation populations. 4. High-elevation genotypes expressed traits within the phenology and WES that limit frost exposure and tissue damage. Specifically, genotypes sourced from high elevations had later mean budburst, earlier mean budset, higher wood densities, higher bark fractions, and smaller xylem vessels than their low-elevation counterparts. Conversely, genotypes sourced from low elevations expressed traits within the LES that prioritized hydraulic efficiency and canopy thermal regulation to cope with extreme heat exposure, including 40% smaller leaf areas, 67% higher stomatal densities, and 34% higher mean theoretical maximum stomatal conductance. Low-elevation genotypes also expressed a lower stomatal control over leaf water potentials that subsequently dropped to pressures that could induce hydraulic failure. 5. Synthesis. Our results suggest that P. fremontii expresses a high degree of coordination across multiple trait spectra to adapt to local climate constraints on photosynthetic gas exchange, growth, and survival. These results, therefore, increase our mechanistic understanding of local adaptation and the potential effects of climate change that in turn, improves our capacity to identify genotypes that are best suited for future restoration efforts.This database contains: Data from 28 traits representing the phenology, leaf, whole-tree architecture, and wood traits spectra from 48 trees grouped into 8 populations and two ecotypes.LES Traits Stomatal anatomy - In 2016, fully expanded leaves were collected from outer leaves on the south facing side of mid-canopy height to assess stomatal density, length, width, and area. Following the nail polish impression method (Hilu and Randall, 1984), 160 impressions on both the abaxial and adaxial sides of the leaves (n = 640 images) were obtained to be observed under an Olympus CX41light microscope and images were taken with a Moticam Pro 282A camera (Motic, Richmond, BC, Canada). Stomatal density was estimated as the number of stomata in eighty 0.59 mm2 digital images at 10× magnification. Stomatal sizes (length × width) were estimated on 800 stomata from digital images at 40× magnification (n= 100 per population) using an open-source imaging program, ImageJ (https://imagej.nih.gov/ij/). Maximum theoretical stomatal conductance to water vapor (mmol m-2 s-1) was calculated from Franks and Farquhar (2001): Gsmax = (1) Where dw is the diffusivity of water in air (2.43x10-5 m2 s-1), v is the molar volume of air (0.024 m3 mol-1) (Jones, 2013), Ds is the stomatal density, amax is the maximum area of the open stomatal pore, approximated as π(p/2)2, where p is stomatal pore length, assumed to be stomatal length divided by two (Franks & Farquhar, 2007). Leaf traits - Specific leaf area (SLA) was calculated as the one-sided area of a fresh leaf, divided by its oven-dry mass (Wright & Westoby, 2002). SLA was measured in June, July, and September of 2017. A subset of 12-20 collected leaves per tree were scanned with a high-resolution computer scanner, and one-sided leaf area was measured with ImageJ. The scanned leaves were then oven-dried for 72 hours at 75 °C and weighed to calculate SLA (cm2 g-1). Individual leaf area (Ail) was derived from the average individual leaf area from these measurements (Ackerly, Knight, Weiss, Barton, & Starmer, 2002). Petiole traits - In September 2017, three to four leaves located at the mid-canopy level and south-facing side were collected from six trees per population to study petiole traits. Leaf samples were stored in a cooler at approximately 7-10 °C and transported to the Imaging and Histology Core Facility at Northern Arizona University. Individual petioles were cut with a razor blade and their mid-portions were sectioned to fit within embedding cassettes (28.5 x 41.0 x 6.7 mm). The samples were prepped using an automated fixation and paraffin embedding process. Specifically, the samples were fixed with formalin, dehydrated with increasing concentrations of undenatured alcohol, cleared with xylene to then be infiltrated with paraffin, and then embedded into the cassette. Cassette blocks were then sliced into transverse sections approximately 5-10 µm thick with a microtome and molded onto positively charged slides, deparaffinized with xylene and rehydrated with decreasing concentrations of undenatured alcohol until rinsing with only DI water, and then stained with 0.1% toluidine blue. Images were produced on a digital light microscope that were subsequently analyzed using ImageJ. Several petiole characteristics were recorded, including the length, width, and area of the entire petiole (Ap). Additionally, we measured petiole xylem vessel diameter (d), vessel density (Dpv), mean and total lumen area of all of the vessels contained in the petiole (Apl), the hydraulically weighted mean vessel diameter (Hdp) was calculated as ∑d5/∑d4 (Scholz et al. 2013; Sperry and Saliendra 1994). Petiole flatness (Lpf) – a distinct characteristic of Populus species – was quantified from the petiole width to length ratio at the mid-rib of the petiole (Lindtke, González-Martínez, Macaya-Sanz, & Lexer, 2013). The petiole lumen fraction (Fpl) was calculated as total Apl per petiole transverse area (Ap) at the mid-rib. Mean petiole theoretical hydraulic conductivity (Kp; mg m MPa –1 s –1 ) was calculated from total petiole vessel lumen diameter using the Hagen–Poiseuille equation (Eguchi et al. 2008; Nobel 2009; Tyree and Zimmermann 2002): Kp = (2) where di is the diameter of a single vessel (m), ρ is water density at 25 °C (998 kg m – 3 ) and ƞw is the viscosity of water at 25 °C (8.9 ×10–10 MPa s) (Eguchi et al., 2008). Theoretical hydraulic conductivity per unit leaf area (Kl, mg m-1 s-1 MPa-1) (Sack, Cowan, Jaikumar, & Holbrook, 2003; Sack & Frole, 2006), was calculated as: Kl = ΣKpAil (3) where Ail is the area of the leaf (m2) attached to the petiole and Kp is the mean petiole theoretical hydraulic conductivity. Additionally, we estimated differences in water use strategies at the population and ecotype levels by dividing leaf area normalized theoretical hydraulic conductivity of the petiole (Kl, mg m-1 s-1 MPa-1) by Gsmax. 2.4 | CR traits In July 2017, we estimated whole-tree leaf area (Al) and sapwood area (As) using population-specific allometric relationships between stem diameter to leaf area through a branch summation approach (Jones, O'Hara, Battles, & Gersonde, 2015; Kenefic & Seymour, 1999). Thus, whole-tree leaf area was estimated per branch by multiplying the mass of all leaves by their respective SLA. Canopy diameters (4 to 8 measurements per tree) and their respective canopy areas together with whole-tree height (H) were measured five times during the 2016 and 2017 growing seasons with a telescoping measuring pole. Canopy area (Ac) was determined using the ellipse equation, πab, where a is the mean radius of longest canopy axes and b is the radius of two perpendicular canopy axes (Ansley, Mirik, Surber, & Park, 2012). Leaf area index (LAI) was estimated by the equation from (Hultine, Burtch, & Ehleringer, 2013): LAI = AlAc (4) 2.5 | WES Traits In May to June 2018, we collected one-year-old stem samples to measure wood and xylem traits. We collected branch cuttings approximately 1 cm diameter and cut the segment to a length of 30 cm. The segments were placed into a plastic bag with a moist paper towel and kept in a cooler until being transferred to a lab refrigerator kept at 4 °C. Specific wood density (Dstem) was determined using Archimedes' principle of water displacement (Cornelissen et al., 2003; Hacke, Sperry, & Pittermann, 2000; Preston, Cornwell, & DeNoyer, 2006). The outer bark was removed from a 1 cm diameter stem segment. The segments were cut to a length of approximately 15 cm with no obvious side branches and total stem volume with the bark was measured. Specifically, a graduated cylinder with water was tared on a scale with 0.01 g accuracy, and the segment was submerged just below the meniscus; the weight was recorded and converted to volume. After the whole stem volume was initially measured, the bark was removed, and the sapwood volume was measured using the same method. The samples were then over-dried at 70° C for 48 h, and dry weight was recorded. Dstem was determined as the ratio of dry weight to the volume displaced of the sapwood. The whole stem volume was measured to assess the stem bark fraction (Fb). Xylem vessel area was determined from the 1 cm diameter segments used for density by first removing 2 cm length segments to mount in a GSL1 sledge microtome (Gärtner, Lucchinetti, & Schweingruber, 2014). The samples were moistened with a Strasburger solution (Eilmann, Zweifel, Buchmann, Graf Pannatier, & Rigling, 2011) before being shaved into 100-150 µm slices. The fresh-cut slices were transferred to a dye solution of 0.1% toluidine blue for 1-2 minutes (Sridharan & Shankar, 2012) before being dehydrated by increasing concentrations of undenatured alcohol (50, 70, 95, and 99.5% EtOH; KOPTEC 200 proof, VWR) (Buesa & Peshkov, 2009). At the end of this process, we used the synthetic mounting medium Permount (Mayr et al., 2014; Ravikumar, Surekha, & Thavarajah, 2014) to embed the samples. Slides were photographed with a Moticam Pro 282A camera (Motic, Richmond, BC, Canada) mounted to an Olympus CX41light microscope. The region of interest included a subset of the previous year's growth ring to include both the earlywood and latewood. Vessel lumen area was measured using ImageJ software, and individual images were stitched together (Preibisch, Saalfeld, & Tomancak, 2009) to analyze xylem vessels of the entire growth year. The number of vessels measured per genotype ranged from 153 to 941 (mean = 460). From the images, we quantified mean stem lumen area (Asl), mean stem hydraulic diameter (Hds) and stem vessel density (Dsv). Stem lumen fraction (Fsl) was calculated as the total Asl per Dsv. Water potential - Monthly measurements of leaf water potential (Ψ) were taken from June to October of 2017 using the Scholander pressure chamber (Cochard, Forestier, & Améglio, 2001; Scholander, Hammel, Bradstreet, & Hemmingsen, 1965). Pre-dawn and midday Ψ measurements were taken to assess possible differences in water potential gradients at the population and ecotype levels. A single shoot tip from each of the 64 trees was cut with a sharp razor blade to measure water potential at predawn (Ψpd) between 0300 and 0500 h, and at midday (Ψmd) between 1100 and 1300 h. 2.6 | STATISTICAL ANALYSIS All statistical analyses were conducted in R version 3.6.2 (R Development Core Team 2011). Prior to analyzing the data, we examined whether each variable met the assumptions of normality and homogeneity of variance, using a Shapiro and Barlett test. When the data were not normally distributed, they were normalized by log10, square root, or box-cox transformations. Once the basic requirements were met, functional traits were analyzed statistically using linear regression with provenance elevation as the predictor and traits as the responses. We also used an analysis of variance (ANOVA) to analyze differences among populations. When a trait showed significant variation, we used Tukey's HSD posthoc test to detect differences at the elevation level (Sokal & Rohlf, 1995). Trait contrasts at the ecotype level were analyzed using a student's t-test. Ecotype differences in traits measured several times during the growing season (Ail, SLA, Ac, H, and Ψpd , Ψmd) were analyzed by individual mixed-model repeated measures ANOVA (type III with Satterthwaite's method) in the lmer R package (Bates, Mächler, Bolker, & Walker, 2015; Kuznetsova, Brockhoff, & Christensen, 2017). In this test, individual traits were represented as response variables while the ecotype and month were treated as fixed effects with two and three levels, respectively. Individual genotype nested within ecotype was the random effect. We conducted separated principal component analyses (PCAs) on the three trait spectra (leaf, architecture and wood) and all traits pooled together using the factoextra and FactoMineR packages (Kassambara & Mundt, 2017; Lê, Josse, & Husson, 2008). Because Gsmax is autocorrelated with Dstom and Sstom, it was excluded from these analyses. Thus, we simultaneously assessed what traits explained most of the variation within each trait spectrum and among all the remaining 27 P. fremontii traits together. For each PCA analysis, the variables showing the highest loading in each PCA were selected as indicators of local climatic adaptability. We initially determined the number of meaningful PCA axes using the Kaiser criterion. This criterion recommends using axes with eigenvalues above 1.0 exclusively. However, because the Kaiser criterion is not recommended to be used as the only cut-off criterion for estimating the number of factors (Freeman & Jackson, 1992; Grossman, Nickerson, & Freeman, 1991; Jackson, 2016; Peres-Neto, Jackson, & Somers, 2005), we also used the Broken Stick Model in the vegan and biodiversity R package to determine significant components. This model randomly divided a stick of components into the same number of elements found in the PCA axis. Then, these elements were rearranged in decreasing order according to their length to be compared to the eigenvalues. Axes with larger eigenvalues than their corresponding stick of components were considered significant (Borcard, Gillet, & Legendre, 2011). In each principal component graph (biplot), trait representation was based on the magnitude of the correlation (loadings) between traits and the given principal component. Thus, in each biplot, traits were represented as vectors with a length and direction indicating the strength and trend of a given trait's relationship among other traits. Specific location of the vector in the biplot indicates the positive or negative impact that a trait has on each of the two components x-axis, first component (PC1) and y-axis, second component (PC2). Additionally, to assess the relationship between the two ecotypes and the traits distribution in every PCA biplot, we constructed two 95% confidence ellipses based on the PCA scores of each of the two ecotype means. Subsequently, linear regressions between significant principal components and elevation of source populations were constructed. Additionally, we performed t-tests and ANOVA Tukey's HSD tests to assess significant differences in PC axes scores at ecotype and population levels. Because PC scores mainly described dominant traits in each axis, we evaluated population differences in all traits included in each PCA simultaneously by using a permutational multivariate analysis of variance (PERMANOVA; (Anderson, 2001) in the vegan R package. Population differences in all 27 traits PCA were further analyzed with pairwise comparisons using permutation MANOVAs on a Pillai test and distance matrix in the vegan and RVAideMemoire R packages. Redundancy Analysis (RDA) (Borcard et al., 2011)) was used to determine how environmental (e.g. provenance latitude, longitude, mean annual precipitation, MAT transfer distance, and ecotype) and genetic (e.g. ecotype) descriptors influence multi-spectra trait variance in P. fremontii. We also used forward stepwise selection of descriptor variables to determine the significance of each variable to the RDA. Forward selection determines the successive contribution of each descriptor to explaining trait variation and adds only those variables with significant contribution