The study site is situated within Tiantong National Forest Park (29^°48′N, 121^°47′E) in Zhejiang Province, eastern China (Figure 1A). This region experiences a typical subtropical monsoon climate, characterized by warm, moist summers and dry, cold winters (Wang et al., 2007). The mean annual temperature is 16.2 ^°C; the warmest and coldest months are July and January, with a mean temperature of 28.1 ^°C and 4.2 ^°C, respectively. This region supports EBLF, primarily dominated by species belonging to the Fagaceae and Theaceae families (Ren et al., 2021).

This study was conducted in a 4.84-ha undisturbed core area (220 m × 220 m) of a 20-ha permanent forest plot established following the ForestGEO protocol (Condit, 1998; Yang et al., 2016). All free-standing trees with a diameter at breast height ≥ 1 cm were tagged, mapped, measured, and identified to species every five years. For this study, a total of 121 quadrats (10 m × 10 m) were uniformly selected to survey and collect tree and understory data (Figure 1B). The woody plant community is notably dominated by species such as Eurya loquaiana and Litsea elongata; detailed information on Importance Values (IV) and leaf habits is provided in Supplementary Table 1.

A total of 45 fern species were recorded in the summers of 2024 and 2025 (Supplementary Table 2). The cover of herbaceous plants in each quadrat ranges from 4% to 95%, and the cover of ferns ranges from 3% to 83%. Ferns were the dominant component of the herbaceous layer, accounting for an average of 69.6% of the total cover in each plot. The distributions of relative fern cover and species richness across the study plots are illustrated in Supplementary Figure 1. The dominance ranking of community, listing the top 20 species ranked by Importance Value (IV), is provided in Supplementary Figure 2.

Six functional traits were measured: fern height (FH), leaf area (LA), specific leaf area (SLA), leaf dry matter content (LDMC), chlorophyll content (SPAD), and actual quantum yield of photosystem II (Φ_PSII). FH serves as an indicator of plant performance (Violle et al., 2007), and variations in height among individuals may reflect asymmetric competition for light (Blondeel et al., 2020; Herben and Austin, 2016). Smaller LA is generally associated with reduced water loss and enhanced tolerance to nutritional stress (Niinemets et al., 2007). SLA, defined as the ratio of total leaf area to total leaf dry mass, reflects a trade-off between mass-based photosynthetic capacity and leaf life span (Wright et al., 2004). LDMC, calculated as the ratio of leaf dry mass to fresh mass, is positively correlated with leaf lifespan and the plant’s ability for nutrient retention (Wright et al., 2004). SPAD plays a vital role in absorbing, transferring, and converting light energy during photosynthesis. Its content and composition are influenced by the light environment (Boardman, 1977; Croft et al., 2017). Φ_PSII quantifies the efficiency of light energy use by the photosystems II (Maxwell and Johnson, 2000).

Traits were measured on mature and healthy fronds from five randomly selected, spatially distinct individuals (or ramets) per fern species in each plot. All trait measurements followed the standard protocols of Pérez-Harguindeguy et al. (2013); for species with dimorphic fronds, all measurements were taken exclusively on the sterile fronds. In total, 2273 samples from 45 fern species were collected to measure their leaf functional traits. FH and SPAD were measured in the field. FH was measured as the vertical distance from the plant base to the apex of the youngest fully expanded leaf, without stretching the axis. SPAD was recorded on fully expanded leaves using a SPAD-502 Plus chlorophyll meter (Konica-Minolta, Japan), with three measurements per leaf averaged (Marenco et al., 2009). In the laboratory, fresh mass and leaf area were determined within 12 hours using an analytical balance and a Li-Cor Portable Area Meter Li-3000 (Li-Cor Biosciences, Nebraska, USA). Φ_PSII was assessed using a highly portable modulated chlorophyll fluorometer (MINI-PAM-II) after 30 minutes of dark adaptation. Leaves were then dried at 75°C for 72 h to obtain leaf dry mass. SLA and LDMC were calculated following standard protocols.

Given that plant functional traits commonly covary due to underlying physiological and evolutionary constraints, a Principal Component Analysis (PCA) was performed on the trait dataset to extract the primary axes of functional variation and reduce dimensionality. The first two principal components (PCs) cumulatively explained 65.7% of the total variation (Figure 2). The first axis (PC1, explaining 40.0% of the variance) was defined as ‘Resource-use Strategies’, representing a gradient from acquisitive to conservative strategies, and was negatively correlated with SLA and positively correlated with LDMC. The second axis (PC2, explaining 25.7% of the variance) represented ‘Plant Growth’, driven mainly by positive loadings for FH and LA.

Abiotic factors were quantified from four topographic and nine soil properties (Table 1). Topographic variables (elevation, slope, aspect and convexity) were calculated for each 10 m × 10 m quadrat (Yang et al., 2011). Soil samples were collected and analyzed for pH, water content, total nitrogen (TN), soil organic carbon (SOC), total phosphorus (TP), potassium (K), calcium (Ca), magnesium (Mg), and manganese (Mn) following John et al. (2007). Point measurements of soil properties were spatially interpolated to all quadrats using standard block kriging (Goovaerts, 1997; John et al., 2007).

To reduce dimensionality and multicollinearity, independent Principal Component Analyses (PCA) were performed for the topographic (4 variables) and soil (9 variables) datasets. The first two principal components (PCs) from the topography and soil PCAs were retained for subsequent analyses; these collectively accounted for 73.8% and 83.3% of the total variance in their respective datasets (Figure 3).

The first topographic PC axis (PC_Topo-ele) was primarily positively correlated with elevation, while the second axis (PC_Topo-asp) was associated with aspect (Figure 3A). The first soil PC axis (PC_Soil-pH) correlated positively with pH and negatively with TN, SOC, and soil moisture. The second soil PC axis (PC_Soil-P) was primarily driven by TP and Ca content (Figure 3B).

The overstory biotic environment was characterized using four variables: tree stand density, basal area, and the first two axes of a non-metric multidimensional scaling (NMDS) ordination of tree species composition (NMDS1 and NMDS2; Table 1). Stand density was calculated as the number of individual trees (DBH ≥ 1 cm) per quadrat. Total basal area was the sum of the cross-sectional stem areas at breast height for all trees in a quadrat (Condit, 1998; Forrester and Bauhus, 2016). An NMDS ordination was performed based on a matrix of species Importance Values (IV) using the “metaMDS” and “MDSRotate” functions in vegan package in R (Fortin and Dale, 2005; Oksanen et al., 2024). The first two NMDS axes (NMDS1 and NMDS2), representing the primary gradients in community composition, were extracted for analysis (Supplementary Figure 3). NMDS1 captured a gradient from simple, low-stature communities to complex, multi-layered canopies in moist habitats. NMDS2 primarily reflected a shift from evergreen to deciduous species, indicating changes in the light environment beneath the forest canopy.

We calculated three types of community-weighted mean (CWM) for each trait and the two PCA axes in each plot, weighted by the relative cover of fern species (Garnier et al., 2004; Wieczynski et al., 2019; Garbowski et al., 2024).

CWM_fixed was calculated based on a fixed mean trait value per species averaged across all plots, reflecting only species turnover effects (Equation 1).

where pi is the weight of the i-th species (cover), S is the number of species, and x_i is the fixed mean trait value (or PC score) of the i-th species averaged across all plots where the species is found.

CWM_specific was calculated using plot-specific mean trait values, incorporating both ITV and species turnover (Equation 2).

where x_im is the specific mean trait value (or PC score) of the i-th species, which is valid just for a given plot sampled.

CWM_intravar, representing the community-level trait variation attributable solely to ITV, was derived as the difference: CWM_intravar = CWM_specific − CWM_fixed (Siefert et al., 2015). In this equation, the sign of CWM_intravar (positive or negative) indicates the direction of the community’s trait shift away from the fixed species mean, while its magnitude reflects the strength of the ITV effect (Table 2).

To partition the sources of community-level trait variation, the relative contribution of ITV and species turnover was quantified following the sum of squares decomposition method (Lepš et al., 2011). This method decomposes the total sum of squares of the CWM_specific values (SS_specific) which represents the total community-level trait variation, into three components: the sum of squares explained by species turnover (SS_fixed), the sum of squares explained purely by intraspecific trait variation (SS_{intraspecific}) and the sum of squares explained by the covariation between them (SS_cov), so that SS_specific = SS_fixed + SS_{intraspecific} + SS_cov. We calculated the percentage contribution of each component to the total explained variation and visualized the results using stacked bar plots. We performed this decomposition for each of the six functional traits.

We used Generalized Least Squares (GLS) models to evaluate the effect of eight environmental predictors (PC_Topo-ele, PC_Topo-asp, PC_Soil-pH, PC_Soil-P, stand density, total basal area, NMDS1, and NMDS2) on community functional composition while accounting for spatial autocorrelation. Separate models were fitted for each of the two functional axes (PC1 and PC2) and the six traits, using CWM_specific, CWM_fixed, and CWM_intravar as response variables. All eight predictors were standardized to interpret their relative importance on a comparable scale (Gross et al., 2017). Potential multicollinearity was assessed by calculating variance inflation factors (VIF), and all VIF values were found to be below 5 (Supplementary Table 3). The overall explanatory power of each model was assessed using Nagelkerke’s pseudo-R-squared.

To further quantify the independent and shared effects of these environmental factor groups (topographic, soil, and biotic), we performed a variance partitioning analysis based on the Nagelkerke’s pseudo-R² values derived from the GLS models (Pinheiro and Bates, 2000; Pinheiro et al., 2024). All statistical analyses were conducted in R (version 4.4.1, R Core Team, 2024). Analyses relied on several packages, including nlme, tidyverse, FactoMineR, vegan, ggplot2 and eulerr. Data exploration, model evaluation, and graphical visualization closely followed the comprehensive framework proposed by Govaert et al. (2024), which provides a robust standard for analyzing trait-environment relationships. We adapted their R code workflow to specifically address the ecological context of understory fern communities.

The total variance in CWM traits was partitioned into components of species turnover and ITV (Figure 4). The contribution of species turnover ranged from 23% to 59%, while that of ITV ranged from 34% to 65%. Species turnover made a greater contribution to the total variance of LA, LDMC, and SPAD, whereas ITV accounted for a larger proportion of variance in FH, SLA, and Φ_PSII. Notably, we observed a negative covariation between species turnover and ITV for SLA and SPAD, resulting in the sum of the independent contributions of species turnover and ITV exceeding 100% of the total community-level trait variance.

Regarding environmental explanatory power, the GLS models revealed distinct drivers for functional community composition. A consistent pattern emerged across both the integrated functional axes (Figure 5) and the six individual traits (Supplementary Figure 4): models based on species turnover alone (CWM_fixed) exhibited consistently higher explanatory power than models including ITV (CWM_specific). For the integrated functional axes, the environmental models explained 32% of the variation for ‘Resource-use Strategies’ (increasing to 44% when considering turnover alone), and 14% for ‘Plant Growth’ (vs. 67% for turnover). Models for CWM_intravar (ITV alone) showed limited environmental explanatory power across integrated functional axes (16% and 6%). For the individual traits, the reduction in explanatory power when including ITV was most evident for SPAD (pseudo-R² decreased from 0.32 to 0.11) (Supplementary Figure 4). Models for CWM_intravar (ITV alone) showed limited environmental explanatory power across all traits (pseudo-R² < 0.11).

Variance partitioning revealed that topographic, soil, and biotic factors influenced the CWM of plant functional traits through both their unique effects and shared effects (Figure 5). For CWM_fixed models, the primary environmental drivers differed markedly between the two axes. The unique effect of biotic factors was the strongest explanatory component for ‘Resource-use Strategies’ (17.7%), while the unique effect of soil factors was dominant for ‘Plant Growth’ (23.0%), followed closely by the unique effect of biotic factors (20.0%). For CWM_specific models, the relative importance of these drivers shifted. The unique effect of soil factors became the largest contributor for ‘Resource-use Strategies’ (14.0%), whereas for ‘Plant Growth’, the unique effect of biotic factors explained the largest proportion of variance (5.7%).

Regarding individual traits, the relative importance of environmental drivers varied considerably (Supplementary Figure 4). For CWM_fixed models, a clear divergence was observed between trait types: biotic factors were the primary drivers for leaf physiological traits, explaining the largest proportion of variation for SPAD (21.7%) and SLA (18.5%). In contrast, soil factors played a leading role for morphological traits, particularly for LA (13.4%) and FH (11.1%). For CWM_specific models, the dominant drivers shifted for certain traits. Notably, the primary driver for SLA shifted from biotic to soil factors (7.9%), while the driver for FH shifted from soil to biotic factors (9.4%). Finally, for LDMC and Φ_PSII, the effects of all three environmental factor groups were negligible (R²≈0) in both model types.

Regarding abiotic drivers, soil factors were the primary predictors for both integrated functional axes (Figure 6) and individual traits (Supplementary Figure 5). PC_Soil-P (positively associated with soil phosphorus and calcium) was negatively associated with the ‘Resource-use Strategies’ axis across both the CWM_fixed and CWM_specific models, while PC_Soil-pH was positively associated with the ‘Plant Growth’ axis in the CWM_fixed model. Reflecting these broad patterns at the individual trait level, PC_Soil-P exhibited the most widespread effects: increasing values were associated with increased SLA, SPAD, and Φ_PSII, and decreased FH, LA, and LDMC. These relationships were generally consistent between the CWM_fixed and CWM_specific models, although the positive association with SLA was significant only in the CWM_specific model, and some associations weakened or became non-significant (e.g., SPAD) when ITV was included. Regarding topographic factors, PC_Topo-ele showed a significant positive correlation with ‘Resource-use Strategies’ axis across both models. In contrast, PC_Topo-asp primarily influenced individual morphological traits, showing associations with increases in FH and LA and a decrease in SPAD in the CWM_fixed model.

Biotic factors also significantly influenced the functional axes and fern traits (Figure 6; Supplementary Figure 5). Stand density showed a consistent negative correlation with the ‘Resource-use Strategies’ axis in the CWM_specific model, whereas total basal area was positively associated with this axis in the CWM_fixed model. NMDS1 exerted a consistently negative influence on the ‘Resource-use Strategies’ axis across both the CWM_fixed and CWM_specific models. At the individual trait level, stand density exhibited consistent negative correlations with FH and LA across both model types. SPAD and Φ_PSII were significantly correlated with stand density only in the CWM_fixed model, while LDMC showed a significant response only in the CWM_specific model. Total basal area influenced traits exclusively in the CWM_fixed model (positive with LA, negative with Φ_PSII). NMDS1 was strongly correlated with four traits in the CWM_fixed model (negatively with LA and LDMC; positively with SLA and Φ_PSII). In the CWM_specific model, only the negative correlation with LDMC remained significant. NMDS2 only showed consistent positive relationships with photosynthetic traits (SPAD and Φ_PSII) in CWM_fixed models.

We found that while ITV accounted for a significant proportion (34–65%) of total fern community trait variance (Figure 4) as some previous studies (Siefert et al., 2015; Wang et al., 2022), its inclusion in community-weighted means (CWM_specific) consistently weakened the explanatory power of trait-environment models compared to models based on species mean traits alone (CWM_fixed) (Figure 5; Supplementary Figures 4, 6). This attenuation of model fit is directly attributable to the limited environmental predictability of the ITV component itself, as evidenced by the low pseudo-R² values for CWM_intravar models across both integrated functional axes and individual traits (Figure 5; Supplementary Figure 4). This suggests that the high magnitude of ITV does not track environmental gradients in the same direction or strength as species turnover. Instead, it likely arises from species-specific (idiosyncratic) behaviors (Westerband et al., 2021; Albert et al., 2010). Our species-level analysis supports this: only a minority of species showed significant trait responses to any given gradient, and the few responding species were influenced by different environmental factors (Supplementary Table 4; Supplementary Figure 7). For example, the SLA of Diplopterygium glaucum responded significantly only to topography, while that of Woodwardia japonica responded to elevation and stand density (Supplementary Figure 5). Similarly, for single trait like LDMC, Diplopterygium glaucum showed no significant relationship with any measured environmental factors, whereas Woodwardia japonica responded to soil P (Supplementary Figure 7). This combination of prevalent non-responsiveness and divergent responses among the few responsive species results in a community-level ITV signal that is largely decoupled from the dominant environmental gradients, explaining its weak aggregate relationship.

This pattern suggests that community-level trait responses are primarily driven by shifts in species composition (turnover) rather than coordinated phenotypic adjustments. One potential interpretation for this dominance of turnover may lie in the hypothesis of evolutionary conservatism in ferns, where the rigid physiological constraints inherent in this ancient lineage may limit individual-level plasticity. This aligns perfectly with findings in other fern communities, where strong species turnover along elevation gradients is interpreted as direct evidence for this environmental determinism (Kessler et al., 2007; Merino et al., 2023). Previous research has further found that species turnover occurs significantly between closely related species, which suggests that the niche of a single species is relatively fixed, making it difficult to cross environmental gradients through self-adjustment and thus necessitating their replacement by pre-adapted related species (Lehnert, 2011; Lehnert et al., 2013; Lehnert, 2014). Alternatively, the high but idiosyncratic ITV we observed may simply reflect a scale mismatch, where individuals respond to fine-scale microenvironmental heterogeneity that was not captured at our 10×10 m plot scale (Albert et al., 2010). Finally, the observed pattern may also be interpreted as a stabilizing mechanism. The fact that the combined contributions of turnover and ITV frequently exceed 100% (Figure 4) reveals a negative covariation between these two components. This suggests that ITV acts as a functional buffer rather than a mere source of noise; as environmental shifts drive fern species replacement, idiosyncratic phenotypic adjustments within species may move in the opposite direction, thereby compensatory stabilizing the community’s overall functional structure (Lepš et al., 2011; Luo et al., 2016; Ferrara et al., 2024). Thus, the apparent independence of ITV from environmental gradients may represent an emergent property of a stabilizing mechanism that maintains functional continuity across varying environments.

Soil gradients, particularly the phosphorus-richness gradient (PC_Soil-P), emerged as a major driver of community functional traits. Along this gradient, we observed a clear, community-wide shift along the Leaf Economics Spectrum (Wright et al., 2004), from a resource-conservative to a resource-acquisitive strategy. This transition was captured by the ‘Resource-use Strategies’ axis (PC1) and manifested as a coordinated adjustment of traits, characterized by a significant increase in SLA and a concurrent decrease in LDMC (Figures 2, 6). The fact that individual trait responses (SLA and LDMC) closely mirrored the shifts in the integrated PC1 axis further reflects the strong biological coupling within the fern leaf economics syndrome (Supplementary Figure 5). This suggests that understory ferns are not merely adjusting traits independently to track soil phosphorus, but are constrained by fundamental trade-offs that force a coordinated functional transition. These coordinated trait shifts strongly suggest that phosphorus acts as a primary limiting nutrient constraining the functional strategies of this fern community. This finding is supported by independent research, which identified this subtropical forest as “highly phosphorus-limited” based on soil respiration experiments (Liu et al., 2019). Our results demonstrate that this understory fern community responds to soil nutrient gradients with a sensitivity similar to that observed in trees (Jones et al., 2013), challenging the generalized view that understory herbs are primarily light-limited (Wang et al., 2022). It highlights the critical role of soil resource limitation even in light-poor understory environments.

Biotic factors derived from the tree community exerted multiple influences on understory fern traits. Tree density showed the most consistent effects across models, indicating that microclimatic variation—primarily understory light availability—acts as a stable biotic filter shaping fern functional strategies (Sercu et al., 2017; Murakami et al., 2022; Jin et al., 2024). Specifically, ferns in denser stands exhibited smaller and thinner fronds, alongside higher SPAD values (Figure 6; Supplementary Figure 5). This conservative pattern is logical in a high-density understory, which creates a deep and persistent low-light environment. Our results show the fern community shifting towards a persistence-oriented strategy: it effectively “forgoes” costly vertical competition (smaller FH, LA) to minimize structural and respiratory costs, while simultaneously maximizing capture efficiency in the deep shade by increasing chlorophyll content (higher SPAD) (Anderson et al., 1995; Johnson et al., 2000; Rundel et al., 2020).

Total basal area, which reflects forest structural size and canopy biomass, affected several traits only in the CWM_fixed model. This suggests that as the forest matures, it imposes powerful, long-term structural filters on the understory, such as the accumulation of deep, recalcitrant litter layers and intense, persistent root competition from large, established trees (Çömez et al., 2021; Müller et al., 2021; Ou et al., 2015). These conditions filter for species pre-adapted to these substrate and microhabitat conditions, while excluding species that rely on different forest-floor environments. This explains why the relationship is visible only in the CWM_fixed model.

Tree species composition also served as a significant biotic filter. NMDS1, which captured a gradient from simple to complex, multi-layered canopies (Supplementary Figure 3), was associated with shifts in fern stature, leaf structure, and photosynthetic traits. NMDS2, reflecting a shift from evergreen to deciduous species, especially pioneer deciduous species like Alangium kurzii and Clerodendrum cyrtophyllum, showed consistent positive relationships with photosynthetic traits (SPAD and Φ_PSII). The disappearance of most NMDS–trait relationships in the CWMspecific model further reinforces that these biotic influences operate predominantly through species turnover rather than ITV.