Materials and methods

Study site and sampling procedure. This study was conducted in the basin of Lake Cuitzeo. It has an area of 4,026 km² and is located at 19° 30’ - 20° 05’ N and 100° 35’ - 101° 30’ W in the Trans-Mexican Volcanic Belt in the northern part of Michoacan state and the southern part of Guanajuato state (Figure 1). Climate in the basin is temperate with a marked rainy season during summer months (June to September). There is a significant climate gradient, with precipitation increasing and temperature decreasing from north to south and with altitude (^{Mendoza et al. 2006}). The topography and soils of the study area are product of the volcanic activity of the Quaternary. The dominant soil groups are vertisols, luvisols, leptosols, acrisols and andosols (^{Mendoza et al. 2006}, ^{Chávez-Vergara et al. 2014}).

Quercus castanea grows throughout the Cuitzeo basin, except in the northeastern portion (^{Aguilar-Romero et al. 2016}). In the study region, populations can be found between 2,000 and 2,800 m (^{Herrera-Arroyo et al. 2013}). Twenty-two populations were chosen for sampling, mostly covering the distribution of the species within the Cuitzeo basin (Table S1, Figure 2). At each population, 10 randomly chosen adult trees separated by at least 20 m from each other were sampled by taking four haphazardly selected sun-exposed branches, with terminal twigs with at least 10 matures leaves with no visible damage, at heights between 2 and 5 m. Only in population 10, fewer trees were collected (6) because of their low availability.

Functional traits. In each individual tree, we measured the following functional traits: leaf chlorophyll concentration (CC), leaf area (LA), leaf mass per area (LMA), leaf thickness (LT) and Huber value (HV, the sapwood cross section area divided by the leaf area distal to the stem) (^{Tyree & Ewers 1991}). Studying these functional traits is relevant because they are adaptive and vary along environmental gradients in many plant species (^{Auger & Shipley 2013}, ^{Rosbakh et al. 2015}, ^{Cochrane et al. 2016}). CC directly influences the photosynthetic capacity of plants (^{Croft et al. 2017}). LA is important in the balance of water and energy of the plant and usually smaller leaves are associated with high radiation, heat-, cold- or drought-stress (^{Cornelissen et al. 2003}). LMA is positively related to the investment in structural defenses and with leaf lifespan and usually higher values of LMA occur in environments with higher resources stress. LT is associated with photosynthetic rates per unit of leaf area and is usually correlated positively with mean temperature and with solar radiation (^{Niinemets 2001}). Finally, HV is higher in more arid sites (^{Carter & White 2009}).

Leaf chlorophyll concentration was measured in the field using the Minolta Soil Plant Analysis Development (SPAD-502) chlorophyll meter (^{Markwell et al. 1995}) in three randomly chosen leaves per individual. For the other traits, fully developed leaves without damage were dried at 55 °C, weighed in an analytical balance, and lamina thickness measured in a portion of the leaf without veins using a Mitutoyo Absolute digital caliper (model 500-172-20) with a 0.01 mm precision. Leaves were then imaged with a flatbed high resolution scanner (Epson V700), and the area of each leaf was determined using the ImageJ software (^{Rasband 2010}). Leaf mass per area was calculated as the ratio between dry mass and leaf area values. These traits were measured in ten randomly chosen leaves per tree. For the Huber value, three twigs per tree were selected. The leaf area was determined by digitalizing all the twig leaves with the high-resolution flatbed scanner and summing their individual areas. The area of the sapwood was determined from the diameter of cross sections of the leafless twigs measured with the digital caliper after removing the bark. This trait was determined for three twigs per individual.

Geographic and environmental variables. Spatial coordinates and altitude of each sampling population were recorded using a global positioning system (GPS) unit (Table 1S). Soil type for each population was extracted from a soil surface (^{Cabrera-González et al. 2010}) using GIS Arc View ver. 3.3 (^{ESRI 1999}). To determine soil water holding capacity, three soil samples were taken in different random directions (keeping a minimum angle of 70° between samples) one meter away from the main trunk of each oak tree sampled using a soil core sampler. In the laboratory, soil samples were dried in an oven at 70 °C for 3 days. Subsequently, 10 g of dry soil from each sample were weighted, then wetted to field capacity and weighted again. Soil water holding capacity (WHC) was calculated as follows:

This variable was determined as an average of values of the three soil samples per oak tree and then an average value per population was obtained.

To characterize the climate at each population, 19 bioclimatic variables derived from monthly precipitation and temperature values (period 1910-2009) were extracted at 30 arc seconds and downscaled using a digital elevation model at 30 m of resolution (^{Cuervo-Robayo et al. 2014}, ^{Correa-Ayram et al. 2017}). Nineteen climatic variables were finally obtained for each population using GIS Arc View v3.3 (^{ESRI 1999}).

Data Analysis. The significance of the differences among populations for each functional trait were evaluated through one-way analyses of variance (ANOVA) using average values of individual trees. To determine how the total variance for each functional trait is partitioned among populations, among trees within populations, and within trees, we estimated the variance components using the restricted maximum likelihood (REML) method. Residual variation was considered to correspond to differences among leaves within trees. In this analysis the full database with individual leaf values was used. A stepwise canonical discriminant analysis (CDA) using individual tree averages for each functional trait was performed to further determine which functional traits have the highest variation among the populations sampled. Finally, to evaluate the degree of association among the different functional traits, a pairwise correlation analysis was performed. These analyses were carried out in JMP 8 (SAS Institute, Cary, North Carolina).

The associations between the population mean values of each functional trait and the environmental variables at each sampling population were evaluated with Spearman’s correlation analyses. Before conducting the correlation analyses, redundancy was reduced among the environmental variables by assessing pairwise correlations and discarding the more specific variable in each pair of highly correlated variables (r ≥ 0.8). The included environmental variables were soil water holding capacity (WHC), annual mean temperature (B1), mean diurnal range (B2), isothermality (B3), minimum temperature of the coldest month (B6), temperature annual range (B7), annual precipitation (B12), precipitation of the driest month (B14), precipitation seasonality (B15) and precipitation of the driest quarter (B17). The JMP 8 program was used in these analyses (SAS Institute, Cary, North Carolina).

To separate the effects of spatial distance and environmental variables on phenotypic variation we performed a Redundancy Analysis (RDA) (^{Van den Wollenberg 1977}), which is a constrained ordination method analogous to linear regression for datasets with multiple dependent and multiple independent variables. Three separate RDAs were conducted: (1) a full model with environmental and geographic variables (latitude and longitude of the populations) as explanatory variables; (2) a partial RDA (pRDA1) where we removed geographic variables (i.e., a pure environmental model) and (3) a partial RDA (pRDA2) where we removed environmental variables (i.e., a pure geographic model). We then used variance partitioning to calculate the proportion of variation explained by the independent contributions of environmental and geographic variables and their joint effect (^{Borcard et al. 1992}, ^{Peres-Neto et al. 2006}). The significance of each RDA model was calculated using a permutation test with 1,000 permutations. We calculated the adjusted coefficient of multiple determination (R ² _adj) for all models (^{Peres-Neto et al. 2006}). Data of the five functional traits were used as the response variables in the three RDAs. The five variables were log₁₀ -transformed to correct for skew, and then centered and standardized before the RDA analysis. The explanatory variables included eight of the environmental variables described above, but we excluded isothermality (B3) because centering and standardizing this variable resulted mostly into zero values. Explanatory variables also included spatial variables defined using principal coordinates of neighborhood matrices (PCNM), which were calculated from the geographic coordinates in decimal degrees (^{Dray et al. 2006}). We retained half of the PCNM variables with positive eigenvalues, as has been suggested by some authors (^{Manel et al. 2012}, ^{Fitzpatrick & Keller 2015}), which were seven in our case. Calculation of PCNM, RDA, test for the significance and adjusted coefficient of multiple determination of the models were performed in R (^{R Core Team 2016}), using the ‘vegan’ 2.3-0 package (^{Oksanen et al. 2015}).

Results

One-way ANOVAs indicated highly significant differences among populations for the five functional traits evaluated (Table 1). However, total variance was partitioned differently depending on the trait. LT showed the largest variation among populations (43.9 %) while for the other four traits variation among trees within populations and leaves within trees accounted for the largest proportion of the variance (Table 2).

The stepwise canonical discriminant analysis indicated that the variable that contributes more importantly to differentiate populations is LT, followed by LMA, LA, CC and HV. The two first canonical discriminant functions allowed highly significant discrimination among populations (Wilks’ lambda = 0.048; P < 0.0001) and explained 71.2 and 16.13 % of the variation, respectively (Table 3). The variables that contributed more strongly to the first canonical function were LT and LMA, while LA, CC and HV contributed to the second canonical function.

We found moderate but significant positive correlations of leaf mass per area with leaf thickness and with the Huber value (Table 4). Chlorophyll content had positive weak correlations with leaf mass per area and the Huber value. In contrast, negative correlations were observed between leaf mass per area and leaf area, between chlorophyll content and leaf thickness and between the Huber value and leaf area (Table 4).

According to the Spearman’s correlation analyses, LMA had a significant negative correlation with annual mean temperature (r _S = -0.5; P = 0.017) and a positive correlation with isothermality (r _S = 0.5; P = 0.016). CC and HV were also positively related to isothermality (r _S = 0.48; P = 0.02 and r _S = 0.43; P = 0.04, respectively), while LT was negatively correlated with precipitation seasonality (r _S= -0.6; P = 0.003).

The full RDA model for the combined effect of geographic and climatic explanatory variables on functional trait variation was significant (P = 0.001). The pRDA1 for the association between environmental variables and functional trait variation while controlling for geographic effects was also significant (P = 0.001, R ² _adj = 13.9 %) (Table 5). In this analysis, the most important variables were annual precipitation (B12) and precipitation of the driest quarter (B17), followed by minimum temperature of the coldest month (B6), annual mean temperature (B1), mean diurnal range (B2), temperature annual range (B7) and precipitation of the driest month (B14) (Table S2). The pRDA2 indicated that geography alone explained a slightly higher proportion of functional trait variation (P = 0.001, R ² _adj = 18.4 %). Finally, joint climate and geography accounted for 71.9 % of the variation (Table 5).

Discussion

Our analysis at a landscape scale showed significant variation among populations of Q. castanea in the functional traits studied. Furthermore, we identified the environmental variables with the largest influence on these functional traits at this spatial scale. It is widely acknowledged that leaf area, leaf mass per area, leaf thickness and Huber value are adaptive and vary according to the environment. These traits are related to water use efficiency and tolerance to temperature and water stress (^{Westoby et al. 2002}, ^{Wright et al. 2004}, ^{Lohbeck et al. 2013}), and leaf mass per area and leaf thickness are functional traits within the leaf economic spectrum (LES) (^{De La Riva et al. 2016}).

Variation among populations. In comparison to other oak species studied across the world, we found that Q. castanea showed a high leaf mass per area (range 134-183 g m^-2), relatively low leaf areas (9.38-18.62 cm²) and intermediate leaf thickness (0.11-0.24 mm) (Table S3). The value of leaf mass per area is situated at the higher end of values reported for oaks, and comparable to values found in evergreen Mediterranean and evergreen arid species (^{Gil-Pelegrín et al. 2017}). According to ^{Flexas et al. (2014)} a leaf mass per area value higher than 120 g m^-2 indicates a true sclerophyllous species. Leaf thickness was within the range of mean values reported for deciduous temperate and evergreen arid oak species. In turn, leaf area was higher than for Mediterranean and arid evergreen species, but lower than for tropical and temperate evergreen and deciduous species (^{Gil-Pelegrín et al. 2017}). These three traits (leaf area, leaf thickness, leaf mass per area) showed highly significant variation among the sampled populations of Q. castanea even though, except for leaf thickness, variation among trees within populations was higher, suggesting the importance of microenvironmental factors (or other unaccounted factors) in determining intraspecific functional trait variation.

Leaf mass per area has been considered a key leaf functional trait for plants (^{Wright et al. 2004}, ^{Poorter et al. 2009}). In Q. castanea populations, leaf mass per area showed a positive correlation with leaf thickness and a negative correlation with leaf area (Table 4). However, leaf thickness was the trait that showed the lower intra-individual and intra-population variation and the higher among-population variation (Table 2), suggesting its importance in the adjustment of Q. castanea to the environment.

In a previous study within the Cuitzeo basin, ^{Aguilar-Romero et al. (2017)} had compared water-use strategies among nine oak species (including Q. castanea) distributed along an aridity gradient. In that study, species from more arid areas had more deciduousness and a higher instantaneous water-use efficiency, while their more humid counterparts had less deciduousness and a xylem that was more resistant to embolisms, suggesting a trade-off between the xylem vulnerability to embolism and deciduousness. In that study, Q. castanea was characterized as occupying an intermediate position along the aridity gradient, with a brevideciduous phenology and also intermediate values (in relation to the other eight species) of wood density, water use efficiency, Huber value and xylem resistance to embolism (^{Aguilar-Romero et al. 2017}). However, the breadth of the environmental niche of each species and its correspondence with the degree of intraspecific variability in functional traits was not considered. As we have seen in this study, such evaluation is critical, since some traits may exhibit even higher intraspecific than interspecific variation. For example, both in greenhouse grown seedlings and in adult trees from the field, ^{Ramírez-Valiente et al. (2015)} found higher differences among populations of the tropical oak Q. oleoides than among species within the same clade (Virentes, the live oaks), for various functional traits in response to drought. Therefore, the extent and ratio of intra- versus interspecific variation, and their relationship with the niche breadth of the species and their coexistence within communities is still an open question (^{Cavender-Bares et al. 2004}).

Association with environmental variables. Even at the relatively fine spatial scale of the Cuitzeo basin, we were able to define several significant associations of the traits measured with environmental variables. The strongest pairwise association was a negative correlation between leaf thickness and precipitation seasonality, meaning that leaves are thinner at populations where precipitation is more seasonal (i.e. there are more marked dry and wet periods). In turn, leaf mass per area had a positive correlation with isothermality and negative correlation with annual mean temperature, and the Huber value had a positive correlation with isothermality. Higher Huber values mean that the tree produces more wood per unit of leaf area, a characteristic that suggests an improvement in its water and nutrient storage capacity (^{Callaway et al. 1994}) and that reduces the vulnerability to embolism of the xylem (^{Tyree & Dixon 1986}).

In turn, the multivariate RDA identified annual precipitation (B12) and precipitation of the driest quarter (B17) as the most important environmental variables influencing functional trait variation of Q. castanea in the Cuitzeo basin, with also a considerable contribution of precipitation of the driest month (B14). Overall, these results indicate that leaves of Q. castanea are more sclerophyllous (i.e. have higher leaf mass per area) in sites that are cooler and less seasonal in precipitation and temperature variation, and less sclerophyllous in more seasonal and warmer sites. Previously, sclerophylly has been considered an adaptation to arid environments across biomes (^{Niinemets 2001}). At the intraspecific level, leaf mass per area has been shown to increase with aridity in Mediterranean oaks such as Q. ilex and Q. coccifera (^{Peguero-Pina et al. 2016}). However, the opposite pattern has been shown in populations of Q. oleoides that tended to have more sclerophyllous leaves in more mesic areas (^{Ramírez-Valiente et al. 2015}, ²⁰¹⁷). The authors suggested that this counterintuitive result is due to phenological differences: populations from mesic populations maintain leaves for longer time during the dry season and therefore maintain function with increasing water stress (^{Cavender-Bares & Ramírez-Valiente 2017}). Therefore, selection has favored an increase in deciduousness (and less investment in tissue construction) with increasing dry season severity in this tropical oak species, contrasting with Mediterranean species that are limited by both cold winters and dry summers (^{Cavender-Bares & Ramírez-Valiente 2017}). Therefore, Q. castanea follows a pattern resembling more a tropical than a Mediterranean species. We hypothesize that phenological patterns (i.e. degree of deciduousness) in this species should be negatively associated with leaf mass per area along populations in the Cuitzeo basin, but this should be evaluated in further studies.

We conclude that, while previous studies have considered functional variation in oak species at a whole-range level, we have shown that significant functional differences exist among populations of Q. castenea separated by a few kilometers in the heterogenous landscape of the Cuitzeo basin. This species shows clearly sclerophyllous leaves but leaf thickness varies to a considerable degree across populations in association with precipitation seasonality, indicating that the temporality of water availability represents a significant environmental pressure for this oak species. Whether the observed functional variation is due to local adaptation, plasticity or a combination of both will be clarified in ongoing common garden and landscape genomics studies.

Supplemental data

Supplemental material for this article can be accessed here: https://doi.org/10.17129/botsci.2449