Trait variation and functional diversity maintenance of understory herbaceous species coexisting along an elevational gradient in Yulong Mountain, Southwest China.

Characterizing trait variation across different ecological scales in plant communities has been viewed as a way to gain insights into the mechanisms driving species coexistence. However, little is known about how changes in intraspecific and interspecific traits across sites influence species richness and community assembly, especially in understory herbaceous communities. Here we partitioned the variance of four functional traits (maximum height, leaf thickness, leaf area and specific leaf area) across four nested biological scales: individual, species, plot, and elevation to quantify the scale-dependent distributions of understory herbaceous trait variance. We also integrated the comparison of the trait variance ratios to null models to investigate the effects of different ecological processes on community assembly and functional diversity along a 1200-m elevational gradient in Yulong Mountain. We found interspecific trait variation was the main trait variation component for leaf traits, although intraspecific trait variation ranged from 10% to 28% of total variation. In particular, maximum height exhibited high plasticity, and intraspecific variation accounted for 44% of the total variation. Despite the fact that species composition varied across elevation and species richness decreased dramatically along the elevational gradient, there was little variance at our largest (elevation) scale in leaf traits and functional diversity remained constant along the elevational gradient, indicating that traits responded to smaller scale influences. External filtering was only observed at high elevations. However, strong internal filtering was detected along the entire elevational gradient in understory herbaceous communities, possibly due to competition. Our results provide evidence that species coexistence in understory herbaceous communities might be structured by differential niche-assembled processes. This approach -- integrating different biological scales of trait variation -- may provide a better understanding of the mechanisms involved in the structure of communities.

Introduction

Functional traits are increasingly being employed to measure and understand species responses to abiotic change, biotic interactions and the link between diversity and ecosystem processes (Diaz et al., 2007, Kraft et al., 2015, Carmona et al., 2016). Characterizing trait variation and linking community trait patterns to causal mechanisms are major research focuses in functional ecology (Messier et al., 2010, Sutherland et al., 2013, Vilà-Cabrera et al., 2015). Many trait-based studies estimate trait values for species while often not adequately accounting for intraspecific trait variation (Carmona et al., 2016). Adequately quantifying the amount of variation in traits is necessary to understand community assembly (Swenson et al., 2007, Kraft et al., 2015), functional trait diversity (Díaz and Cabido, 2001) and species coexistence (Wright et al., 2004, McGill et al., 2006, Ackerly and Cornwell, 2007).

While environmental gradients are often employed in ecological research, we have little understanding of how traits vary across different biological scales, such as within an individual, within a species, among species, and among communities, and further how these relate to species persistence and diversity (Messier et al., 2010, Violle et al., 2012, Carmona et al., 2016). For example, Messier et al. (2010) reported that two commonly studied leaf traits (i.e. leaf mass area and leaf dry matter content) showed no variation at a plot level when analyzed across six biological scales, highlighting the role of environmental filtering in determining plant community structure.

Integrating inter- and intraspecific trait variation has provided important insights into the processes of community assembly (Jung et al., 2010, Lepš et al., 2011, Siefert, 2012, Siefert et al., 2015, Luo et al., 2016). ‘Internal’ and ‘external’ filters have been distinguished from ‘traditional’ ecological filter for community assembly. Internal filters represent all local processes, including micro-environmental heterogeneity or density-dependent processes that regulate species coexistence within the community. External filters can select specific species from regional species pools through environmental filtering (Violle et al., 2012, Taudiere and Violle, 2015). Trait variance ratios (T-statistics metrics, including TIP/IC, TIC/IR, TPC/PR) provide a useful method to decompose the trait variance across biological scales and can be used to draw inferences about hypotheses of community assembly related to mechanisms that select for key traits and promote trait similarity vs. those that select for dissimilar species (Taudiere and Violle, 2015). Specifically, the ratio TIP/IC is the variance within a single species relative to the total variation in the community over all species, which was designed to estimate of the effects of internal filters on trait overlap among coexisting species (Violle et al., 2012, Bello et al., 2013). Meanwhile, the ratios TIC/IR and TPC/PR can be used to assess the importance of external filters on communities at individual and species levels, respectively (Violle et al., 2012).

The TIP/IC metric represents a useful estimator of the trait overlap between species (Violle et al., 2012). Studying the relationship between TIP/IC and species richness might provide useful insights into the mechanisms of community assembly (Violle et al., 2012). For example, a negative richness-TIP/IC relationship would indicate niche-based coexistence via partitioned niche space. However, if trait overlap is not related to species richness, then we should expect that either communities are assembled according to the neutral theory, or that the measured traits are not reflecting assembly mechanisms. Finally, a positive richness-TIP/IC relationship may be more complicated to explain. This pattern can reflect two distinct hypotheses. First, the ‘individual variation’ theory explicitly identifies intraspecific variation as the main driver of local diversity (Clark et al., 2007, Clark, 2010). In contrast, if equalizing mechanisms, such as herbivory or disease, minimize average fitness differences between species, then coexistence may only require small niche differences between species, allowing for greater within-species variability (Chesson, 2000, Adler et al., 2007, Mayfield and Levine, 2010, Bagousse-Pinguet et al., 2014). Thus, it is also important to consider changes in TIP/IC to disentangle possible interactions between biotic and abiotic forces driving community assembly (Violle et al., 2012, Bagousse-Pinguet et al., 2014).

Elevational gradients provide the necessary variation to test ecological hypotheses (Körner, 2007). For our purposes, the elevational gradient supplies substantial environmental changes over a relatively short distance, and allows us to observe how community assembly mechanisms shape functional trait patterns. In this study, we partitioned the variance of understory herbaceous functional traits across four biological scales (individual, species, plot, elevation). We also explored whether the variation in functional traits was correlated with species richness. Furthermore, we determined whether species and functional diversity of coexisting herbaceous species in understory communities changed consistently along a 1200-m elevational gradient in Yulong Mountain, located in a global biodiversity hotspot with dramatic variation in topography, habitat and climate (Luo et al., 2016). Specifically, we addressed: 1) how functional traits of understory herbaceous species varied across different biological scales; 2) how plant species and functional diversity vary along an elevational gradient; and 3) how trait variation at different ecological scales relates to community assembly and species diversity. This study provides an understanding of the distribution of diversity and community assembly of understory plant communities along an elevational gradient.

Material and method

Study area

This study was conducted on the southern edge of Yulong Mountain, located near Lijiang in northwest Yunnan province, China. Yulong Mountain belongs to the Hengduan Mountains region, which is known as the central global biodiversity hotspot of “Mountains of Southwest China” (Myers et al., 2000). This region, with a dry season (November to May) and rainy season (June to October), is dominated by the southwest monsoon and Indian Ocean monsoon. The mean annual temperature is 12.8 °C and the annual precipitation is 935 mm (Feng et al., 2006).

The vegetation in this region is highly heterogeneous along the elevational gradient. Evergreen temperate conifers, sclerophyllous broad-leaved and coniferous/sclerophyllous broad-leaved mixed forests dominate the forests in this region. The dominant understory herbaceous species are Roscoea cautleoides, Viola delavayi, Rubus fockeanus, Ranunculus yunnanensis, Polygonum paleaceum, Ainsliaea foliosa and Triplostegia glandulifera.

Field survey

We placed fifteen 20-m × 50-m forest plots along elevational gradients (Luo et al., 2016). Each plot was divided into ten subplots (10 m × 10 m). Then we randomly set ten 1-m2 quadrats nested within each subplot for the inventory of herbaceous vegetation (Fig. 1). Herbaceous species composition and abundance within the quadrats were recorded during the growing season in 2013. Species richness and composition in quadrats within plots were combined. In total, 175 understory herbaceous species were collected and the species richness of each plot is shown in Table 1.

Fig. 1

Map of the study region in Lijiang, Yunnan, China (a). Plot locations (white triangle) are shown along the elevational gradients in Yulong Mountain (b). Each plot was subdivided into ten subplots, with ten 1-m2 quadrats where the understory herb species were monitored (gray squares) (c).

Table 1

The plot information and number of samples for understory herbaceous species in Yulong Mountain.

Plot code	Latitude	Longitude	Elevation (m)	Species richness	Number of samples
					Leaf traits	Heightmax
Plot01	26.9998	100.1980	2670	73	864	172
Plot02	26.9981	100.1989	2650	66	538	205
Plot03	26.9991	100.1984	2665	74	844	169
Plot04	26.9995	100.1911	2960	51	686	83
Plot05	26.9997	100.1909	2950	48	702	92
Plot06	26.9992	100.1906	2965	48	706	102
Plot07	27.0045	100.1819	3250	43	652	50
Plot08	27.0006	100.1808	3260	25	342	54
Plot09	27.0003	100.1780	3280	29	392	72
Plot10	27.0097	100.1757	3524	31	370	78
Plot11	27.0101	100.1756	3550	31	362	88
Plot12	27.0107	100.1751	3540	28	374	79
Plot13	27.0183	100.1761	3840	23	360	33
Plot14	27.0186	100.1775	3850	22	266	54
Plot15	27.0179	100.1770	3830	21	238	45
Total				175	7696	1376

Trait data

Four functional traits were selected in this study: maximum height (Heightmax), leaf thickness, leaf area, and specific leaf area (SLA). These four traits are believed to influence different plant functional strategies. Height is related to light availability and competition (Wright et al., 2005); leaf thickness plays a key role in determining leaf physical strength and longevity (Wright et al., 2004, Pérez-Harguindeguy et al., 2013); leaf area relates to energy and water balance and tolerance to environmental stress (Cornelissen et al., 2003, Pérez-Harguindeguy et al., 2013). SLA captures species strategies for acquiring, using and conserving resources, including light, nutrients and water (Wright et al., 2004). For height, we measured the maximum height for each species that occurred in each quadrat. For leaf trait measurements, 14 to 20 leaves were used from seven to ten non-damaged and mature individuals from each species in each plot and these measures took petioles into account. In total, 1376 individuals were measured for height and 7696 leaves were sampled in this study (Table 1). Leaf thickness (mm) was measured at the center of the lamina by avoiding major leaf veins by using electronic digital calipers. Leaves were scanned with a CanonScan LiDE 210 (Canon Inc., Tokyo, Japan) and the leaf area was calculated using Image J (Abràmoff et al., 2004). SLA was calculated as leaf area divided by leaf dry mass (after the leaf was dried to a constant weight at 70 °C).

Data analysis

To check the completeness of our sampling, species richness accumulation curves (method Chao{Chao, 1992 #1717}) (Chao and Lee, 1992) for sampled herbaceous species was applied along the number of quadrats. Furthermore, to visualize the composition difference among different elevations we used Nonmetric Multidimensional Scaling (NMDS) based on Bray–Curtis distance.

To explore the variance components of functional traits across four nested scales (individual, species, plot and elevation), we fitted a linear mixed model using a restricted maximum likelihood (REML) method and a variance component analysis. Confidence intervals of 95% for the percentage of variance explained at each nested scale were calculated by bootstrapping (999 runs with sample size randomly with replacement for specific trait). Sample size was 1376 for Heightmax, 7696 for leaf thickness, 7696 for leaf area and 7696 for SLA. Since large trait values have a greater influence on the arithmetic mean, making them more prone to sampling error (Bland and Altman, 1996), all functional traits data were normalized by log10 transformations.

Trait distribution of species in different elevations were estimated with a kernel density (Mouillot et al., 2005). The overlap as the integral of the intersection area of the kernel distribution of the traits was calculated, and the average pairwise overlap of traits at different elevations was also compared. To gain more insight into how trait distributions varied along the elevational gradient, the distribution of each trait was described through the mean, range, variance, kurtosis, and standard deviation of the distribution of neighbour distances (sdNND). To explore the influence of spatial autocorrelation between the plots, we calculated the Moran's I value of each trait distribution. The correlation strengths between trait distribution and elevation were estimated using spatial autoregressive (SAR) models correcting for spatial non-independence among plots (Kissling and Carl, 2007).

Furthermore, to compare the differences in mean species richness and functional diversity between the five elevations, ANOVA and Tukey's HSD post hoc pairwise comparisons were performed. Functional diversity including functional richness (FRic), functional evenness (FEve) and functional dispersion (FDis) (Laliberté and Legendre, 2010) were calculated. FRic represents functional space occupied by number of species, FEve corresponds to how regularly distribution in functional space (Mason et al., 2005), and functional dispersion is the mean distance in multi-dimensional functional trait space of individual species to the centroid of all species within a community (Laliberté and Legendre, 2010). Functional dispersion and functional evenness are both unbiased by species richness (Laliberté and Legendre, 2010, Mandle and Ticktin, 2015).

We computed three observed trait variance ratios (TIP/IC, TIC/IR and TPC/PR) for each plot and each elevation (Violle et al., 2012). The three ratios were calculated as follows:

TIP/IC = σ2IP/σ2IC, the ratio of within-population variance to total within-community variance, reflects internal filtering affecting individuals and niche packing among the species of the community.

TIC/IR = σ2IC/σ2IR, the ratio of community-wide variance to total variance in the regional pool, assessed external filtering strength when accounting for individuals.

TPC/PR = σ2PC/σ2PR, the ratio of community-wide variance to total variance in the regional pool, assessed external filtering strength when accounting for species (no intraspecific variation).

Where σ2IP is the variation of trait values among individuals within population, σ2IC is the variation of trait values among individuals within community, σ2IR is the variation of trait values among individuals within the regional pool, σ2PC is the variation of population mean trait values within community, and σ2PR is the variation of population mean trait values within the regional pool.

To compare observed patterns with null models, three null models were created as following: 1) TIP/IC: randomization of individual trait values within the community; 2) TIC/IR: randomized without replacement of individual trait values belonging to all plots; 3) TPC/PR: assigned a plot-level value to each individual and drawn without replacement of plot level trait values belonging to all plots. More details about null models can be viewed in supporting information (Table S1). All these null models were run 999 times. The standardized effect sizes (SES) were calculated as: SES = (Iobs−Inull)/sd null, where Iobs is the observed value, Inull is the mean of the null model, and sdnull is the standard deviation of these null values. To interpret the relationship between TIP/IC and species richness, linear regressions were used. All analyses were conducted in R version 3.2.4 (R Development Core Team, 2015) and the following packages vegan (Oksanen et al., 2011), varComp (Qu, 2015), nlme (Pinheiro et al., 2015), spdep (Bivand, 2012), cati (Taudiere and Violle, 2015), and FD (Laliberté and Shipley, 2011) were used for data analysis in this study.

Results

The shape of the rarefaction curve describing the accumulated species richness was saturated across the increasing number of quadrats (Fig. S1). The variation partition in leaf traits across four ecological scales revealed that very little variation was explained at the elevation and plot levels (Table 2), despite the fact that species composition was different across the elevation gradient (Fig. 2). Interspecific variation accounted for 82%, 83%, and 68% of total variance of leaf thickness, leaf area and SLA, respectively; while intraspecific variation represented 17%, 10%, 28% of the total variance of leaf thickness, leaf area and SLA, respectively (Table 2). The partitioning revealed that the largest proportion (44%) of variance in height was explained at the intraspecific level, followed by 38% of variance at the interspecific level and 18% at the elevation level (Table 2).

Table 2

Variance partitioning of the nested linear models for Heightmax, leaf thickness, leaf area and SLA across nested ecological scales. Individuals nested within species, species nested within plot, and plot nested within elevation. All data were log10 transformed prior to analyses, sample size was 1376 for Heightmax, 7696 for leaf thickness, 7696 for leaf area and 7696 for SLA.

Scale	% variance of functional traits [95% CI]
	Heightmax	Leaf thickness	Leaf area	SLA
Elevation	18 [2, 34]	0 [0, 0]	7 [6,8]	0 [0,0]
Plot	0 [0, 0]	1 [0, 2]	0 [0, 0]	4 [4,5]
Species	38 [34, 42]	82 [80, 84]	83 [79–87]	68 [52,84]
Individuals	44 [33,55]	17 [11, 23]	10 [2–16]	28 [20,36]

Fig. 2

Nonmetric Multidimensional Scaling (NMDS) for the dissimilarity using Bray–Curtis distance in understory herb composition among different elevational gradients. Circle dot indicate study plot, cross bars indicate the species composition (Stress value is 0.093).

Four traits showed similar distribution patterns at different elevations (Fig. 3). Pairwise trait overlaps between the elevations tended to have relatively high overlap (Table 3). Furthermore, for most traits there was not a significant relationship between trait distribution and elevation (Fig. S2). Species richness decreased dramatically at high elevations (Fig. 4a), however, the measures of functional diversity (functional evenness and functional dispersion) were consistent across the elevation gradient (Fig. 4b).

Fig. 3

Distributions of traits (height, leaf thickness, leaf area and SLA) constructed from kernel density for different elevations.

Table 3

Pairwise overlaps in trait values compared among elevations based on kernel density.

Traits	Height	Leaf thickness	Leaf area	SLA
Pairwise elevations
2650–2950	0.76	0.78	0.69	0.79
2650–3250	0.83	0.88	0.74	0.89
2650–3550	0.82	0.94	0.68	0.96
2650–3850	0.93	0.96	0.66	0.97
2950–3250	0.65	0.9	0.77	0.91
2950–3550	0.64	0.93	0.67	0.94
2950–3850	0.79	0.97	0.68	0.97
3250–3550	0.66	0.84	0.61	0.81
3250–3850	0.70	0.92	0.64	0.86
3550–3850	0.82	0.84	0.57	0.83

Fig. 4

The changes of herbaceous species richness (A) and functional diversity (B) along elevational gradients. Error bars denote SE, different letters represent significant differences from LSD (least significant difference) comparisons (P < 0.05).

The TIP/IC metric, which compares within-plot intraspecific to interspecific variation, was significantly lower than null expectations for all four traits (Fig. 5A and D). The individual (TIC/IR) and species (TPC/PR) level trait variance ratios that compare plot patterns to the regional species pool were significantly lower than null expectations at high elevations (Fig. 5B and C). Conversely, there were no significant patterns for the individual (TIC/IR) and species (TPC/PR) level trait variance ratios at low elevations (Fig. 5B and C), meaning that they were not different from null expectations based on randomization.

Fig. 5

Standardized effect size (SES) of TIP/IC (A), TIC/IR (B) and TPC/PR (C) for the four functional traits along elevational gradients and three trait variance ratios across all plots (D). The boxes indicate the confidence interval of the null model for each trait variance ratio. Each colored dot represents the SES value of one community when it is deviated from null model. The crossed circles and segments represent the mean and the standard deviation of the SES values for a given trait variance ratios and a given trait, respectively. For a given trait variance ratio, the mean of the SES (crossed circle) is significantly different from the null distribution if not embedded within the colored.

The observed TIP/IC values for Heightmax were positively related to species richness (R = 0.28, P = 0.04) (Fig. 6), indicating that greater height overlap occurred in high species richness plots. However, no significant relationships were found for leaf thickness, leaf area or SLA (Fig. 6).

Fig. 6

The relationships between TIP/IC and species richness for four traits including height, leaf thickness, leaf area and SLA. Solid lines indicate when significant.

Discussion

Our results reveal that despite extremely high species turnover across an elevation gradient, functional diversity remained surprisingly consistent. Furthermore, the distribution of trait values within plots relative to the regional species pool was not significantly different than null expectations at low elevations. These two findings reveal that community functional diversity is robust to environmental gradients and comprised of redundant species at low elevations. However, despite these larger scale patterns, intraspecific trait variation within communities was an important component of overall trait variation. Trait overlap between species within community deviated significantly from null models, and the overlap between species increased with greater species richness for height. Meanwhile, the individual and species level trait variance ratios were significantly lower than null expectations at high elevations, which indicated that external filtering greatly influenced individuals and species at high elevations.

Trait variation at different ecological scales

Trait-based approaches have been increasingly adopted to disentangle the different processes in community ecology (Kraft et al., 2008, Messier et al., 2016). The relative importance of interspecific and intraspecific variability might depend on the species and traits considered (Albert et al., 2010), thus different patterns of trait variation can be observed at different ecological scales or different regions (Roche et al., 2004, Albert et al., 2010, Hulshof and Swenson, 2010, Messier et al., 2010, Auger and Shipley, 2013, Kang et al., 2014). For instance, the relative importance of interspecific SLA variation was 63% in Costa Rican tropical dry forest (Hulshof and Swenson, 2010), however, the relative importance was only 21% in Panamanian tropical forest (Messier et al., 2010). Until now no study has investigated understory herbaceous trait variation at different ecological scales. Our results showed that there was almost no variation in foliar traits at large scales. This lack of variation likely reflects the importance of filtering processes acting on understory plants where the most important limiting factor should be light availability. While the intraspecific variability of leaf traits was lower than the interspecific variability, it was not negligible. The relative importance of interspecific variation was 82%, 83%, 68% for leaf thickness, leaf area and SLA, respectively, while the relative of importance of intraspecific variation was 17%, 10%, 28% for the same traits. Thus, the contribution of interspecific variation to total variation for the three leaf traits is much higher than that of intraspecific for the herbaceous species.

In this study, maximum height also exhibited more plasticity than leaf morphological traits, with intraspecific variation accounting for 44% of total variation for height. Further, intraspecific variation was an important component at the largest scale, accounting for 18% of total variance across the elevation gradient. We previously observed similar patterns for canopy trees in this region (Luo et al., 2016), revealing that these patterns are remarkably consistent across different life history strategies. Further, these results reinforce the fact that height, which has been correlated with the vertical light niche (Falster and Westoby, 2003), is the important factor in this system. Thus, intraspecific variation in height helps ensure constant fitness in light-limited environments.

Trait variation and species diversity along an elevation gradient

The three ratios of trait variances at different ecological scales were proposed to test the trait variance in community assembly (Violle et al., 2012). Higher TIP/IC suggests an increase in trait overlap between species (Violle et al., 2012). Overall, our results shed light on the positive relationship between TIP/IC ratio and species richness for height; as such, species diversity was associated with greater overlap in height between species. A positive relationship between trait overlap and species richness implies that light competition might be a potential driver in structuring grasslands (Grime, 2006, Bagousse-Pinguet et al., 2014). On the one hand, increasing of intraspecific variance of height may decrease the opportunity for individuals from taller species to overcome smaller species in the process of asymmetric light competition (Semchenko et al., 2012). On the other hand, previous studies have shown that individuals from different species tended to have similar heights in order to decrease the gap in competitive ability (Mayfield and Levine, 2010, Spasojevic and Suding, 2012). In this study, the competition for light availability could also be an important determinant of species coexistence in understory, especially in high species richness plots at low elevation. Meanwhile, species richness was not related to functional diversity, as ‘equalizing fitness’ processes may shape plant height in structuring understory herbaceous communities as well.

No significant linear relationship was observed between TIP/IC and species richness for leaf thickness, leaf area or SLA. However, different trends emerged for leaf thickness and SLA. The SLA TIP/IC ratio increased with increasing species richness, while leaf thickness decreased with increasing species richness. This indicates that different functional traits were potentially associated with different processes in communities.

Functional maintenance and species coexistence

Using trait variance ratios of the understory herbaceous species, the presence of internal and environmental filtering on these communities were able to be determined (Violle et al., 2012, Taudiere and Violle, 2015). The magnitude of TIP/IC indicated the strength of internal filters (Violle et al., 2012, Neyret et al., 2016), and in our case the observed TIP/IC was significantly lower than randomized values for all four traits. Thus, we conclude that species exhibit more packing in trait space than null expectations, and this indicates that biotic interactions (such as competition) or micro-environmental heterogeneity is determining local community composition. The effect of internal filtering on communities at low elevations was stronger than at high elevations. This possibly resulted from the stronger competition in more species-rich plots.

The individual (TIC/IR) and species (TPC/PR) level trait variance ratios were slightly lower than expected from null models at high elevational plots, therefore species in these plots might be limited by environmental filters. It has been reported that environmental filtering can shape the functional composition of grasslands (Siefert, 2012, Moraes et al., 2016). However, in this study, functional diversity, especially functional evenness and functional dispersion, was not related to species richness, indicating common limiting environmental factors like shading from trees that were not correlated with elevation. The assemblages of understory herbaceous communities might be structured by differential niche-based processes. For instance, strong internal filtering was detected in understory herbaceous communities along entire elevational gradients, and external filtering was also observed at high elevations. It is essential that future work combine information about abiotic (e.g. climatic, edaphic and topographic) variables to explore how environmental factors shape community assembly.

Maintaining functional diversity and structure in herbaceous assemblages along elevational gradients may be critical for ecosystem services. A major focus of future studies should be the identification of how species loss alters functional diversity and structure under future climatic change or human disturbance at high elevations. This would be valuable information for conservation strategies for a subalpine forest system in a global biodiversity hotspot.