Latitudinal variation of leaf stomatal traits from species to community level in forests: linkage with ecosystem productivity

To explore the latitudinal variation of stomatal traits from species to community level and their linkage with net primary productivity (NPP), we investigated leaf stomatal density (SDL) and stomatal length (SLL) across 760 species from nine forest ecosystems in eastern China, and calculated the community-level SD (SDC) and SL (SLC) through species-specific leaf area index (LAI). Our results showed that latitudinal variation in species-level SDL and SLL was minimal, but community-level SDC and SLC decreased clearly with increasing latitude. The relationship between SD and SL was negative across species and different plant functional types (PFTs), but positive at the community level. Furthermore, community-level SDC correlated positively with forest NPP, and explained 51% of the variation in NPP. These findings indicate that the trade-off by regulating SDL and SLL may be an important strategy for plant individuals to adapt to environmental changes, and temperature acts as the main factor influencing community-level stomatal traits through alteration of species composition. Importantly, our findings provide new insight into the relationship between plant traits and ecosystem function.

Changes in the environment can affect ecosystem processes and function directly through influences on abiotic controls, and indirectly through influences on the physiology, morphology, behavior of individual plants, the structure of populations, and the composition of communities [1][2][3] . Recent studies have suggested that the spatial and temporal variation of ecosystem functions is difficult to be explained sufficiently by classical climatological or biogeographical approaches, and plant traits have the potential to advance the understanding in this field 2,4 . Plant traits have been well recognized as important determinants of ecosystem function 2,5 , and it is possible to aggregate functional traits measured on organisms to explain the function of populations, communities, ecosystems, and beyond 2,4-6 . Thus, understanding plant traits and their controlling factors at a large scale is of importance to explore how organisms and environmental conditions collectively shape the variation of ecosystem function in space and time 2 , and to predict ecosystem function in response to environmental changes.
Spatial and temporal variation of ecosystem function, such as net primary productivity (NPP), evapotranspiration, and water use efficiency, has received much attention to date 2,7,8 . Some studies have demonstrated that leaf traits, including leaf size, specific leaf area (SLA), and tissue nitrogen and phosphorus concentration, are strongly correlated with ecosystem productivity 2,3,5 . Abiotic factors (e.g., soil, climate, and disturbance) may exert indirect effects on the variation of NPP through regulating leaf Scientific RepoRts | 5:14454 | DOi: 10.1038/srep14454 for SD L and 9.34-93.76 μ m for SL L , respectively (Fig. 2). Moreover, the variation of SD L was larger than that of SL L (coefficient of variation, CV = 0.70 vs. 0.38, Fig. 2). For nine forest communities, the mean values of community-level SD C and SL C were 1,497.09 stomata mm −2 and 36,032.13 μ m mm -2 (see Supplementary Table S2), with CVs of 0.55 and 0.49, respectively. SD L and SL L varied markedly across different plant functional types (PFTs) ( Fig. 3 and Supplementary  Fig. S1). Compared with shrubs and herbs, the leaves of trees had denser but smaller stomata (F = 234.80, P < 0.001 for SD L ; F = 185, P < 0.001 for SL L ). Among trees, broadleaves had higher SD L and lower SL L than conifers (F = 13.48, P < 0.001 for SD L ; F = 11.73, P = 0.001 for SL L ), and evergreen broadleaves had higher SD L and less SL L than those of deciduous counterparts (F = 24.10, P < 0.001 for SD L ; F = 57.83, P < 0.001 for SL L ).
Relationships between SD and SL at the species and community levels. Strong negative relationship between SD L and SL L was observed across species (Fig. 4A) and different PFTs (Fig. 5), in which SL L decreased linearly with the increase in SD L (after log 10 transformation). The results of standardized major axis (SMA) tests showed that the slopes differed significantly among growth forms (P < 0.05, Fig. 5A and Supplementary Table S3), and the slope for herbs was steeper than those of trees or shrubs. SMA tests for common slopes revealed no significant difference between evergreen and deciduous trees  (P > 0.05), but difference was observed in the y-axis intercept, with higher SL L at any given SD L in deciduous trees than that in evergreen trees (P < 0.05, Fig. 5B and Supplementary Table S3).
In contrast to the results at the species and PFT levels, a significantly positive relationship between SD C and SL C was observed at the community level (R 2 = 0.97, P < 0.001, Fig. 4B).

Latitudinal variation in stomatal traits and their influencing factors.
At the species level, plant species at lower latitudes had higher SD L than those distributed in higher latitudes (R 2 = 0.08, P < 0.001, Fig. 6A), whereas the observed latitudinal trend of SL L was opposite (R 2 = 0.07, P < 0.001, Fig. 6A). When scaling up to the community level, SD C and SL C decreased nonlinearly with increasing latitude (R 2 = 0.83, P < 0.05 for SD C ; R 2 = 0.76, P < 0.05 for SL C , Fig. 6B).
Mixed-effect models were used to explain the effects of environmental factors and PFT on the stomatal traits (Table 1 and Supplementary Table S5). At the species level, the best model for SD L included PFT, mean annual temperature (MAT), soil N, and their interactions as predictors; the best model for SL L only included PFT and MAT (Supplementary Table S5). PFT could explain the largest part of the variation in SD L and SL L (15.3-22.8%) at the species level, and the environmental variables accounted for only a small proportion (4.7-5.4%, Table 1), although the relationships between MAT or soil N and stomatal traits were significant ( Fig. 7A-C). In addition, the interactions among PFT, MAT and soil N accounted for a large proportion of the variation in SD L (16.6-17.9%, Table 1).
To partition the variance of SD L and SL L into among-site and within-site components, the data were further analyzed through nested random-effect analysis of variance (ANOVA), with species nested within  sites. We found that > 50% of the variance in SD L and SL L existed within sites, and 10.3-13% of variance occurred across sites (Supplementary Table S6).
At the community level, both SD C and SL C increased significantly with the increase in MAT (P < 0.05, Fig. 7D,E). The best models of SD C and SL C only included MAT (Supplementary Table S5), which accounted for 35.2-36.7% of total variation (Table 1). Furthermore, the site explained 30.9% and 45.9% of the variation in SD C and SL C , respectively (Table 1).   Relationships between community-level stomatal traits and NPP. Community-level SD C and SL C were closely correlated with forest NPP (Fig. 8). NPP increased linearly with increasing SD C and SL C (P < 0.05, Fig. 8). Moreover, in the general linear model analysis, SD C explained 51% of the total variation in forest NPP at the large scale (Supplementary Table S7), whereas SL C did not have a significant effect on NPP (P > 0.05, Supplementary Table S7), probably because of the high self-correlation between SD C and SL C (R 2 = 0.97, Fig. 4B).

Discussion
This study is the first to explore the latitudinal variation in stomatal traits at the species and community levels from tropical monsoon rainforest to cold-temperate coniferous forest. Furthermore, our findings establish a relationship between stomatal traits and the spatial variation in ecosystem NPP, and thus provide useful biological information for modeling ecosystem carbon and water cycles under global changes.  As expected, community-level stomatal traits were closely related to forest NPP in our study (Fig. 8). In previous studies, tight relationships between aboveground NPP and other community-aggregated traits (such as SLA and leaf N) have been documented 5,22 . Community structure is the basis for scaling up plant traits from species to community level 5 . Furthermore, leaf area index (LAI) has been widely considered as an important parameter controlling the interception of light and water, transpiration, photosynthesis, carbon, and nutrient cycles 3 . In this study, we developed a new methodology to determine community-level stomatal traits through leaf-level measurements weighted by species-specific LAI. Such defined community-level stomatal traits were useful to qualify the total stomatal number and length per unit ground surface area within communities, and could reflect the potential to exchange water vapor and CO 2 between plant community and the atmosphere. The results from Fig. 8 and Supplementary Table  S7 showed that community-level SD C positively correlated with forest NPP and accounted for 51% of spatial variation in NPP along the NSTEC. These findings suggest that the stomatal data measured at the species level may be associated with NPP through community structure, which provides new evidence for the linkage between plant traits and ecosystem function 2 . Moreover, community-level SD C can better characterize the capacity of CO 2 and H 2 O exchange between plant community and atmosphere, and thus may be a new ecological parameter in modeling CO 2 and H 2 O exchange at the community scale.
Because of measurement difficulties and other limitations, few studies have systematically analyzed the biogeographical variation in stomatal traits across wide gradients of environmental factors. In this study, we found that latitudinal variation in species-level SD L and SL L was minimal (R 2 ranging from 0.07 to 0.08, Fig. 6A), and environmental variables weakly modulated stomatal traits (4.7-5.4%, Table 1). The main reason for these weak correlations was that a substantial proportion of the variation in SD L and SL L was observed within sites (> 50%, Supplementary Table S6) and thus resulted in relatively weak spatial patterns of species-level stomatal traits along the large-scale environmental gradients. Similarly, high within-site variation is normal for global-scale analyses of plant morphological, chemical, anatomical, and physiological traits 23,24 . The large variance explained by within-site differences implies that site scale (or community level) is of particular relevance when studying mechanisms of species assembly 25 .
Changes in PFTs were dominant factors influencing the biogeographic distribution of species-level stomatal traits, explaining the largest part of the variation in species-level SD L and SL L ( Table 1). Considerable variation was observed among different PFTs (Fig. 3), reflecting the differences in adaptive strategies of PFTs to their external environment. For example, smaller and denser stomata in trees are associated with higher stomatal conductance and transpiration rates 10,26 , which benefits water and nutrient transmission through longer xylem pathways in woody plants 27 . In contrast, large stomata are critical for herbs to optimize light capture in the light-limited environment 9 . Differences in stomatal traits between broadleaves and conifers reflect the differences in leaf vascular architecture, where broadleaved angiosperms reach high stomatal conductance with numerous small stomata, and conifers obtain lower stomatal conductance with fewer large stomata 28 . Compared with deciduous trees, evergreen broadleaves with high leaf SD L can likely acquire high diffusive conductance for water vapor or CO 2 , thereby meeting the high transpiration requirement and simultaneously providing efficient cooling of the leaf 9,29 . Thus, changes in PFTs along environmental gradients contributed substantially to generating the observed variation in SD L and SL L .
Shifts of PFTs accounted for the major part of the latitudinal variation in SD L and SL L ; however, climate and soil, especially MAT and soil N ( Fig. 7A-C), may directly or indirectly influence stomatal traits by shaping the vegetation biogeography (changes in species composition) 3,30 . In this study, the interactive effect of MAT and PFT accounted for 16.6% of the variation in SD L , and that of PFT and soil N explained another 17.9% of the variation (Table 1).
In this study, the latitudinal trends of SD C and SL C were significant at the community level (P < 0.05, Fig. 6B). Both SD C and SL C decreased with latitude (Fig. 6B), and could be explained largely by MAT (Table 1 and Supplementary Table S5). These results indicated that through regulation of species composition and PFTs within a community, temperature acted as the main environmental factor driving the latitudinal patterns of community-level stomatal traits, even though the influence was moderate at the species level. On the one hand, evergreen broad-leaved forests with higher LAI and SD L also have higher community-level SD C and SL C , allowing them to attain effectively high stomatal conductance to meet the high requirements of transpiration and photosynthesis at low latitudes 3 . On the other hand, coniferous forests dominate at high latitude and cold habitats with low SD L (thus low SD C and SL C ) that allows them to minimize water loss by transpiration and mitigate water-deficit stress resulting from low air and soil temperature 31 .
Stomatal morphological traits and behaviors collectively determine the balance of CO 2 uptake for plant photosynthesis against water loss by transpiration 9 . In addition to the short-term dynamics resulting from altering the width of the stomatal pore apertures, plants can respond to environmental change by adjusting their stomatal morphology via long-term evolution and development. In this study, SL L correlated negatively with SD L across all species and PFTs (Figs 4A and 5), which was consistent with a general relationship reported by previous studies across multiple species and geological time scales [9][10][11] . The trade-off between SD L and SL L can be partially explained by physical and energetic constraints 11 . The physical (space) constraints hypothesis assumes that plants adjust their stomatal size and density to optimize stomatal conductance while satisfying a given stomata-to-pavement cell ratio, and thus the allocation of leaf epidermal space to stomata is limited. The second plausible hypothesis for the trade-off Scientific RepoRts | 5:14454 | DOi: 10.1038/srep14454 between SD L and SL L is the energetic constraint, which emphasizes maximizing the return in terms of stomatal conductance and photosynthesis for a given investment in the construction of stomatal complexes 11 . High stomatal conductance may be related to increased metabolic cost, which may explain why some species revert to producing fewer but larger stomata, particularly for plants grown in low light or low nitrogen availability with a low requirement for stomatal conductance and photosynthesis 11 .
The negative relationship between SD L and SL L governs short-term (plastic) and long-term (evolutionary) adaptations of plant physiological activities to the environment 11 . In contrast, a positive relationship between SD C and SL C was observed at the community level, indicating that the general trade-off between species-level SD L and SL L did not exist between the community-level SD C and SL C . Therefore, we assume that the trade-off mechanism through regulation of SD L and SL L may be an important strategy at the species level for adaptation to the external environment, while the adaptation of community-level stomatal traits is achieved mainly through the modulation of community composition. This mechanism could compensate for the limitation of biological conservatism resulting from the trade-off between leaf SD L and SL L at the species level.
In conclusion, community-level SD C and SL C decreased significantly as latitude increased, although latitudinal variation in species-level SD L and SL L was minimal. The adaptive strategies of stomatal traits to the external environment may be different between the species and community levels. We assume that the trade-off between SD L and SL L represents adaptive mechanism of leaf stomata to respond to environmental changes at the species level. While community-level stomatal traits regulate species composition to maximize ecosystem productivity in a given environment, and to some extent compensate for the limitation of biological conservatism resulting from the trade-off between SD L and SL L . More importantly, SD C could explain 51% of the variation in forest NPP, which bridges stomatal traits to ecosystem function. Finally, the new framework, which scales up stomatal traits from species to community level weighted by LAI, makes it possible to better parameterize the complex models pertaining to ecosystem carbon and water cycles in the future.

Methods
Site description. The NSTEC is a unique forest belt driven mainly by a thermal gradient, and includes almost all forest types found in the northern hemisphere 32 (Fig. 1); therefore, it offers an ideal condition to test the ecological and evolutionary responses of plants to environmental changes 32 . Nine natural forests along the NSTEC were selected to conduct field sampling, including Jianfengling (JF), Dinghu (DH), Jiulian (JL), Shennongjia (SN), Taiyue (TY), Dongling (DL), Changbai (CB), Liangshui (LS), and Huzhong (HZ) (Fig. 1 and Supplementary Table S1). These ecosystems span latitudes from 18.7 to 51.8 °N, with MAT ranging from − 3.67 to 23.15 °C and mean annual precipitation (MAP) ranging from 472.96 to 2265.80 mm. The soil type varies from tropical red soils with low organic matter to brown soils with high organic matter. Correspondingly, vegetation type varies among tropical monsoon rainforest, subtropical evergreen forest, temperate deciduous forest, and cold-temperate coniferous forest. To minimize anthropogenic disturbances, we set up the sampling plots within well-protected national nature reserves in each forest type, where the vegetation is relatively homogenous and is strongly representative of the given forest ecosystems.
Field sampling and measurement. The field survey was conducted during July and August of 2013. Three or four experimental plots (30 × 40 m) were set up in each forest ecosystem. The geographic information (latitude, longitude, and altitude), plant species identity, and community structure were determined for each plot. The number, height, diameter at breast height (DBH) of all trees with DBH ≥ 2 cm (basal stem diameter for shrubs), and aboveground live-biomass of all herbs were measured. Subsequently, 20 fully expanded sun leaves were collected from four individuals of each plant species using a long chain saw or by two professional climbers, and all of the leaf samples for a given plant species from one plot were mixed together as a replicate. Collected fresh leaves were stored in sealed plastic bags and scanned to determine leaf area within 4-8 h after collection. In total, 904 species-at-site combinations, consisting of 760 plant species from 427 genera and 154 families, were collected. Species were also grouped in different PFTs by the possible groups: growth form (herbs, shrubs, and trees), leaf type (coniferous vs. broadleaved trees), and leaf habit (evergreen vs. deciduous broadleaved trees). At each plot, soil samples (0-10 cm) were collected randomly from 30-50 points using a 6-cm-diameter auger, and then mixed thoroughly.
We randomly selected ten fresh leaves for each plant species and scanned them using a scanner (CanoScan LiDE 110, Japan), then measured the area by using Photoshop CS (Adobe Systems, San Jose, USA). Finally, the leaves were oven-dried at 70 °C for 48 h and then weighed to calculate SLA (m 2 kg −1 ) 33 .
Leaf SD L (stomata mm -2 ) and SL L (μ m) were measured from surface impressions of the mid-blade abaxial leaf surface made with clear nail polish (methods described in Wang et al. 34 ). All stomatal measurements were conducted with electronic image analysis equipment (COIC XSZ-HS3 and MIPS software, Optical Instrument Co., Ltd., Chongqing, China). In total, 4284 leaf cuticle images were obtained (see examples in Supplementary Fig. S1).
Calculation of community-level stomatal traits. Here, we defined the total stomata number per unit ground surface area within a community as community-level SD C (stomata mm -2 ) and total Scientific RepoRts | 5:14454 | DOi: 10.1038/srep14454 stomatal length per unit land area as community-level SL C (μ m mm -2 ). Community-level stomatal values were calculated as the sum of species-specific SD Li and SL Li weighted by their LAI (m 2 m -2 ) in each plot (30 × 40 m). To estimate LAI, we firstly calculated the foliar biomass of woody species based on the mean DBH (basal stem diameter for shrubs) and height through allometric equations, which were obtained from Chinese Ecosystem Research Net (CERN) database (http://159.226.111.42/pingtai/cernc/index.jsp), published studies, and our previous field measurements (unpublished data) (available in Supplementary  Table S8). The use of species-specific equations is preferred in the calculation of foliar biomass; however, this is virtually unrealistic in practice, especially for tropical rainforests. Therefore, we used the biomass equations from the same genera, similar PFT, or mixed-species equations of a forest to estimate species' foliar biomass when their allometric equations were unavailable. In total, 246 species-specific equations with R 2 -values ranging from 0.52 to 1.00 were used in this study. For herbs, the foliar biomass was represented by aboveground biomass. Then, species-specific foliar biomass multiplied by SLA was used to calculate the LAI i 35 . Finally, community-level SD C and SL C were calculated according to equations (1) and (2). The frameworks to scale up stomatal traits from species to community level were also illustrated in Supplementary Fig. S2.
Ecosystem NPP. The ecosystem NPP (g C m -2 yr -1 ) used in this study was derived from Luo 36 , who estimated the biomass and NPP of 668 stands from 7 major forest types in China. In Luo's published dataset, NPP was computed as the sum of the annual net increment of root, stem, branch and leaf components of trees, and those of shrubs and herbs. For trees, the annual net increments of stems, branches and roots were obtained by multiplication of the biomass of different components and their growth rate.
The annual net increment of leaves was estimated as the leaf biomass divided by leaf age of different trees. The average NPP of the shrub and herb layers was estimated by dividing their biomass by average stand age. A more detailed description to estimate forest NPP is presented in Ni et al. 35 and Luo 36 . From this dataset, we extracted NPP data based on latitude, longitude, and forest type. Because NPP data in Luo's dataset was based on dry matter, the dry matter was converted into units of C, using a conversion factor of 0.5 37 .
To test the accuracy of these data, we compared ecosystem NPP and LAI with the MODIS products in a 1 × 1 km grid (http://daac.ornl.gov/cgi-bin/MODIS/GLBVIZ_1_Glb/modis_subset_order_global_col5. pl). MODIS NPP data from 2000 to 2010 and LAI data from July to August of 2013 for nine sampling sites were selected. The observed closely relationships between measured NPP and LAI and MODIS products proved our data reliable to some extent ( Supplementary Fig. S3).
Climate data. The climatic variables in this study, including MAT and MAP, were extracted from the meteorological database produced by CERN 38 . The climate dataset at a 1 × 1 km spatial resolution was generated from the data of 740 climate stations of the China Meteorological Administration during 1961 and 2007, using the interpolation software of ANUSPLIN.
The monthly averaged insolation (kWh m -2 day -1 ) of nine sampling sites over a 22-year period (1983-2004) was obtained from the NASA Surface Meteorology and Solar Energy dataset (http://eosweb. larc.nasa.gov/sse/) to represent light intensity in this study.
Analysis of soil properties. Soil water content (SWC, %) was determined after fresh soil was dried at 105 °C for 24 h. Other fresh soil samples were air-dried and sieved with the roots removed by hand, then ground to pass through a 2-mm sieve. Soil total N concentration (mg g -1 ) was analyzed using an elemental analyzer (Vario MAX CN Elemental Analyzer, Germany). Soil total P concentration (mg kg -1 ) was measured by the ammonium molybdate method using a continuous-flow analyzer (AutoAnalyzer3 Continuous-Flow Analyzer; Bran Luebbe, Germany) after H 2 SO 4 -H 2 O 2 -HF digestion 39 . Data analysis. Stomatal data were log 10 -transformed prior to analyses to obtain the approximate normality and homogeneity of residuals. We analyzed stomatal traits at the species level (averaged within species at each site) and community level (weighted by species-specific LAI), respectively.
One-way ANOVA with Tukey's post-hoc test was used to compare differences in species-level SD L and SL L among different PFTs. Then, the bivariate analyses of SD L -SL L and SD C -SL C relationships were performed using ordinary least squares (OLS) linear regressions, and differences in the slopes and elevations of SD L -SL L relationships among PFTs were evaluated by SMA estimation with the R package smatr.
In SMA analyses, we tested the heterogeneity of regression slopes, and estimated the common slopes when non-heterogeneity of the slopes was demonstrated (P > 0.05). Then, differences in the elevation of regression slopes (i.e. y-axis intercept) were tested using ANOVA with Tukey's post-hoc tests.
To explore the latitudinal patterns of stomatal traits and their environmental drivers, we performed regression analyses between stomatal traits and latitude at the species and community levels, respectively. Then, we analyzed each of the stomatal traits using linear mixed-effect models with the residual maximum likelihood (REML) method in the R package lme4. In these analyses, we treated environmental variables and PFT as fixed effects, and the random-effect term for site was included to account for the non-independence of species occurring at the same site 23 . This term also allowed us to determine how variance in stomatal traits was partitioned within and across sites. To avoid problems of collinearity among the environmental variables, only variables with insignificant Pearson's correlation (P > 0.05, Supplementary Table S4) were included. The environmental variables that had significant effects (P < 0.05) on stomatal traits and interaction terms between PFT and the environmental variables were included in the final model (Supplementary Table S5). Models with lower Akaike's Information Criterion (AIC) values were chosen as the final best-fit models 40 .
Finally, linear regressions were performed to describe the relationships between ecosystem NPP and community-level SD C and SL C , and the relative contributions to the spatial variation of NPP were assessed by a general linear model.
All analyses were conducted using SPSS 13.0 statistical software (SPSS Inc., Chicago, IL, USA, 2004) and R software (version 2.15.2, R Development Core Team 2012).