Mixed Stands of Larix principis-rupprechtii and Betula platyphylla Provide Higher Ecosystem Multifunctionality Than Corresponding Pure Forests

1. Introduction

Ecosystem multifunctionality (EMF) refers to the ability of ecosystems to simultaneously provide a range of functions and services [1]. In forest ecosystems, this concept is referred to as forest multifunctionality [2], which encompasses regulating, provisioning, supporting, and cultural functions [3]. Timber production, a key objective of forest management [4], has led most research to focus on quantifying stand productivity (SP) [5]. The intricate interactions between various trophic levels, as well as between above- and below-ground components, form the foundation for the realization of various ecosystem functions [6]. Plant functional diversity significantly influences these interactions among plants and other trophic levels [7], serving as a key driver across multi-trophic levels. Soil organic carbon (SOC) accumulation and nutrient recycling [8] are essential for supporting and regulating forest ecosystems. Compared to other ecosystems, there are fewer studies on EMF, especially in complex, dynamic systems like forests [9]. Historically, research has tended to focus on specific aspects of forest ecosystem functions. Therefore, a comprehensive evaluation of EMF, alongside an analysis of its influencing factors, is essential. Differences exist in the functions provided by different stand types [10]. For example, Tian et al. [8] revealed that mixed forests offer higher ecological functions and EMF compared to pure forests. Monoculture plantations are susceptible to both biotic and abiotic factors [11], whereas mixed forests exhibit higher biodiversity and SP [12], enabling them to provide more ecosystem functions and services than pure forests [13].

Biotic (biodiversity) and abiotic (topography, climate, soil, etc.) factors influence EMF [14]. Biodiversity includes aboveground components, such as species, structural, functional, and phylogenetic diversity [15], as well as belowground components like soil microbial diversity. It primarily affects EMF through two mechanisms: the niche complementarity effect and the mass ratio effect [16]. The former demonstrates that increased biodiversity leads to greater functional diversity among tree species, improving resource use and space efficiency, thereby enhancing ecosystem functions [17]. In contrast, the mass ratio effect suggests that the trait of dominant species dictate the performance of forest ecosystem [18]. For example, thicker trees have been shown to positively influence EMF [19]. Li et al. [20] identified the niche complementarity effect as a potential mechanism for EMF. However, Yue et al. [21] found that in their study of biodiversity-EMF relationships in secondary forests in Northeast China, variations in multifunctionality were primarily influenced by the mass ratio effect.

Climate significantly influences tree growth, with mean annual temperature (MAT, ℃) and mean annual precipitation (MAP, mm) in the growing season being key variables [22]. These factors impact ecosystem functioning both directly and indirectly. Directly, they accelerate the activity of consumers and decomposers; indirectly, they alter community composition [23]. As confirmed by Wang et al. [24], higher MAP was negatively correlated with EMF, while higher MAT directly limited the EMF of Pinus yunnanensis secondary forests. Ecosystems with higher EMF levels are generally more resilient to environmental changes [25], whereas pure forests with a single tree species tend to be more susceptible to environmental change [26]. Changes in soil stoichiometric characteristics can influence the rates of various ecosystem functions. Changes in pH can cause soil acidification and inhibit the decomposition of organic carbon, affecting plant growth and ecosystem functions [23]. Different types of litter can improve soil quality, typically benefiting ecosystems [27]. Additionally, altitude and slope can influence forest diversity and function by altering soil physicochemical properties and microclimate conditions [28]. These findings suggest that EMF levels vary across stand types [1]. Therefore, studying changes in ecological function and EMF across forest types is crucial for optimizing EMF and ensuring sustainable forest management practices [29].

This study investigated three forest types: Betula platyphylla pure forest (BP), B. platyphylla- Larix principis-rupprechtii mixed forest (BL), and L. principis-rupprechtii pure forest (LP). Considering the influence of both biotic and abiotic factors, a theoretical framework outlining the relationships between these factors was constructed (Figure 1). The primary objectives were to: (1) quantify and analyze the relationships between forest productivity, functional diversity, nutrient cycling, soil fertility maintenance, and variations in EMF across these forest types; (2) elucidate the roles of biotic and abiotic factors in EMF at different thresholds; and (3) identify the main mechanism driving EMF. The following three hypotheses were tested: H1: Mixed forests exhibit higher productivity than pure forest; H2: Biodiversity and environmental factors jointly determine forest EMF; and H3: The primary mechanism driving EMF is niche complementarity.

2. Results

2.1. Comparisons of Single Function and EMF Across Forest Types

Compared to pure forests, mixed forests exhibited significant advantages in SP, FDis, FRic, and RaoQ (p < 0.05; Figure 2a). BP was significantly higher than LP in SF, NC, and FEve (p < 0.05; Figure 2a). The EMF of BL was the highest among the three forest types (Figure 2a). As the thresholds increased, the number of functions declined. BL outperformed pure forests across all thresholds (Figure 2b), especially at higher thresholds. This indicated that mixed forests offered greater benefits in maintaining EMF.

2.2. Effects of Biotic and Abiotic Factors on Single Function

NC and SF were markedly influenced by both biotic and abiotic factors (p < 0.05; Figure 3). Sdba exhibited a positive correlation with FEve, NC, and SF (p < 0.05; Figure 3a), but a negative correlation with other individual functions. TD was negatively correlated with all individual functions. MAP and soil nutrients had a significant positive influence on NC and SF (p < 0.05; Figure 3b), while MAT negatively affected all functions except SF.

2.3. Effects of Biotic and Abiotic Factors on EMF

The linear regression model showed that, except for TN and LPC, most indicators had significant positive effects on EMF at the 25% threshold (p < 0.05; Figure 4a). At the 50% threshold, TP, CWM_SLA, LNC, and SP exhibited a strong positive relationship with EMF (p < 0.05; Figure 4b). At the 75% threshold, an increase in the nutrient content of leaves significantly enhanced EMF (p < 0.05; Figure 4c). However, At the 90% threshold, no function exhibited a notable effect on EMF (Figure 4d).

The structural equation model results illustrated that, at different thresholds (from 25% to 90 %), the predictors accounted for 97%, 93%, 31%, and 29% of the variation in EMF, respectively (Figure 5). SoilPC1 exhibited a notable direct positive correlation with EMF at the 25%, 50%, and 75% thresholds (p < 0.05; Figure 5), and it indirectly influenced EMF through Sdba (Figure 6), which exerted a more substantial overall effect on EMF. Climate had a significant positive correlation with EMF at lower threshold levels (p < 0.05; Figure 5). Nevertheless, pH had no notable effect on EMF at all thresholds (Figure 5). TD and RaoQ significantly enhanced EMF at lower thresholds (25-50%) (p < 0.05; Figure 5a-b), while Sdba was negatively correlated with EMF at the 25% threshold (p < 0.05; Figure 5a).

Figure 5. Path diagram of the structural equation model between EMF and driving factors. The figure displays the effects of biotic and abiotic factors on forest EMF. Panels (a) - (d) represent EMF at 25%, 50%, 75%, and 90% thresholds (EMF_25, EMF_50, EMF_75, EMF_90), respectively. pH is soil pH, SoilPC1 is soil nutrient content, Climate is a composite variable of MAT and MAP. RaoQ is Rao’s quadratic entropy index, TD is tree diversity and Sdba is the standard deviation of basal area. In the figure, red line represents a significant path (p < 0.05), black line represents non-significant path, solid line indicates positive impact, and dotted line represents negative impact. Model fit statistics: Chi-Squared = 3.814, P = 0.282, AICC = -111.471. *, p < 0.05, * *, p < 0.01, * * *, p < 0.001.

Figure 5. Path diagram of the structural equation model between EMF and driving factors. The figure displays the effects of biotic and abiotic factors on forest EMF. Panels (a) - (d) represent EMF at 25%, 50%, 75%, and 90% thresholds (EMF_25, EMF_50, EMF_75, EMF_90), respectively. pH is soil pH, SoilPC1 is soil nutrient content, Climate is a composite variable of MAT and MAP. RaoQ is Rao’s quadratic entropy index, TD is tree diversity and Sdba is the standard deviation of basal area. In the figure, red line represents a significant path (p < 0.05), black line represents non-significant path, solid line indicates positive impact, and dotted line represents negative impact. Model fit statistics: Chi-Squared = 3.814, P = 0.282, AICC = -111.471. *, p < 0.05, * *, p < 0.01, * * *, p < 0.001.

Figure 6. The direct and indirect influence of abiotic and biotic factors on EMF at different thresholds. Panels (a) - (d) represent EMF at 25%, 50%, 75%, and 90% thresholds (EMF_25, EMF_50, EMF_75, EMF_90), respectively. All abbreviations are as explained in Figure 5. Green bars indicate indirect effects, blue bars represent direct effects, and yellow bars represent total effects.

Figure 6. The direct and indirect influence of abiotic and biotic factors on EMF at different thresholds. Panels (a) - (d) represent EMF at 25%, 50%, 75%, and 90% thresholds (EMF_25, EMF_50, EMF_75, EMF_90), respectively. All abbreviations are as explained in Figure 5. Green bars indicate indirect effects, blue bars represent direct effects, and yellow bars represent total effects.

3. Discussion

3.1. Comparisons of Single Function and EMF Across Forest Types

In this study, mixed forest showed a significantly higher SP than pure forest (p < 0.05; Figure 2a), supporting our hypothesis H1. This phenomenon may be attributed to the contrasting shade tolerance and canopy structure between B. platyphylla and L. principis-rupprechtii. The niche complementarity between these species reduces interspecific competition compared to intraspecific competition, thereby enhancing resource utilization efficiency and stand productivity [30]. Similarly, Feng et al. [26] found that mixed coniferous and broad-leaved forests exhibit higher stand productivity due to stronger interspecific complementary effects. Therefore, mixed larch and birch forests may be a more effective way to improve stand productivity. FD, a key component of biodiversity, represents the distribution of plant functional traits within multidimensional niche spaces [1] and has received significant attention in EMF studies. Species mixture effects tend to be more pronounced in conifer-broadleaf mixtures due to differences in leaf morphology, as well as variations in SLA, LDMC, LNC, and LPC [31]. The FDis, FRic, and RaoQ values of BL were significantly higher than those of pure forests, a result consistent with Li et al. [32], further emphasizing the importance of FD in mixed forest ecosystems. In our study, BP and BL exhibited significantly higher soil fertility maintenance functions compared to LP (P < 0.05; Figure 2a). This may be due to the wider root distribution in broad-leaved and mixed forests compared to pure coniferous forests, as well as the increased microbial activity in the soil, which enhances litter decomposition [33]. Higher SOC levels are also associated with greater SF [34]. Soil and leaf nutrients jointly regulate forest nutrient cycling capacity. In this study, BP exhibited significantly higher NC than LP (p < 0.05; Figure 2a). This can be attributed to higher nitrogen and phosphorus levels, as well as the larger SLA and higher nutrient content of leaves, which improve photosynthesis capacity [31]. As a broad-leaved tree species, birch leaves store more organic matter. Nevertheless, no significant effect of stand composition on NC was observed, which was different from the conclusion of Li et al. [35]. This may be due to the predominant presence of larch in the selected mixed forests. The EMF of BL was the highest among the three forest types (Figure 2a) and outperformed pure forest at all thresholds (Figure 2b), confirming hypothesis H1. Our results align with most previous studies comparing artificial pure forests and mixed forests. Mixed forests can offer a greater variety of ecosystem functions [29]. Higher tree species and stand structure diversity are essential for maintaining EMF [1]. Given that tree species compositions in plantations are often simplified, selecting suitable tree species for mixing is crucial to improve the EMF of forests.

3.2. Effects of Biological Factors on EMF

Under different thresholds, various aspects of diversity can serve distinct functions. Stand structure attributes (Sdba) were positively correlated with EMF_50, EMF_75, and EMF_90, while a notable negative correlation was found with EMF_25 (p < 0.05; Figure 3). Additionally, Sdba significantly influenced multiple single functions, such as NC and FC (p < 0.05; Figure 4a), highlighting the important role of stand structure in maintaining EMF. These results can be attributed to the fact that trees of different sizes have different resource requirements, and stand structure attributes allow for the occupation of different growth spaces, thereby reducing direct competition and forming a heterogeneous canopy structure [6,36]. This spatial differentiation facilitates more efficient resource use, including light and space, thus enhancing key forest ecosystem functions [20].

Biodiversity often plays a positive role in enhancing forest ecosystem functions [37]. At 25%, 50%, and 75% thresholds, a positive association was found between TD and EMF (Figure 5), although the strength of these correlations differed across various thresholds. The multiple functions of forest are mainly driven by tree species diversity. As tree species diversity decreases, its effect on EMF diminishes, and the mass ratio effect weakens [1]. At the 90% threshold, a negative correlation was observed between TD and EMF, likely due to the fact that with the increase of threshold, only some indicators have higher performance, leading to lower resource availability and intensified competition, which ultimately reduces EMF [36].

Functional diversity indicators can reflect more differences in traits and niches [38]. TD largely influences ecosystem function by altering the composition of functional traits among trees [39]. In this study, RaoQ was used as the measure of FD. The results showed a positively correlated between RaoQ and EMF at all thresholds (Figure 5), with mixed forest exhibiting higher EMF compared to pure forest (Figure 2). This result aligns with findings by Li et al. [20] and supports the niche complementary effect [40]. The leaf morphology of B. platyphylla and L. principis-rupprechtii differed significantly, and mixing these species results in more efficient light utilization, positively impacting ecosystem function and stability [41]. Linear regression analysis revealed that specific tree leaf traits were significantly associated with EMF (p < 0.05; Figure 4) at low and medium thresholds (<75%). As the thresholds increased, EMF exhibited a decreasing trend. According to the mass ratio hypothesis, the characteristics of dominant species significantly shape ecosystem functions [6]. Leaf traits, which reflect a tree’s capacity to capture light energy and nutrients, promote resource complementary and contribute to ecosystem function [42]. Highly diverse ecosystems allow a greater variety of functional traits to complement each other, enhancing species’ efficiency in resource utilization.

At the 25% threshold, structural attributes negatively impacted EMF (p < 0.05; Figure 5). However, both TD and FD exhibited strong positive influences on EMF_25 and EMF_50 (p < 0.05), with no significant impact on higher EMF. At lower EMF levels, these two effects exerted a synergistic influence on EMF [8], which is contrary to Hypothesis H3. As thresholds increased, the number of functions decreased, and biodiversity did not significantly influence EMF, weakening the niche complementarity effect. At the 75% threshold, the mass ratio effect emerged as the primary driving mechanism for EMF. Therefore, in order to achieve multi-functional forest management, it is essential to enhance biodiversity by selecting suitable tree species for mixing and increasing stand structure complexity. This approach would help increase EMF.

3.3. Effects of Abiotic Factors on EMF

The relationship between environment and EMF is both dynamic and complex [43], particularly under the influence of global climate change [8]. Our study indicated a significant direct association between climate and EMF at low thresholds (p < 0.05; Figure 5). Adequate water and light can promote plant growth and improve soil microbial activity, thereby enhancing the EMF of the stand. Notably, previous studies have revealed that the effect of biodiversity on EMF is likely to be influenced, or even altered, by environmental factors [23]. Specifically, climate mediated this effect through TD and stand structure attributes, altering the relationship between biodiversity and EMF at the 25% and 90% thresholds (Figure 5-6).

Soil plays a crucial role in tree growth, serving as a main driver of EMF [23]. At the 25% thresholds, TP, AP, and AN exhibited significant positive impacts on EMF (p < 0.05; Figure 4a). At the 50% thresholds, TP and AP maintained similar effects (p < 0.05; Figure 4b). However, as the threshold increased, the influence of soil nutrients on EMF diminished. At the 75% and 90% thresholds, this influence became negligible (Figure 4c-d). Additionally, the SEM results indicated that soil nutrients positively impacted EMF at all thresholds, with a significant increase observed alongside FD (p < 0.05; Figure 5). These findings show that richer soil nutrients promote tree growth, thereby enhancing overall forest ecosystem function (Manuel Villa et al., 2020). At different thresholds, a notable negative correlation was found between SoilPC1 and Sdba (p < 0.05; Figure 5), indicating that soil nutrients may indirectly affect forest ecosystem function by affecting stand structure [44].

4. Materials and Methods

4.1. Study Area and Sample Plots

The research was conducted in the Mulan Forest Farm (41°35′~42°40’ N,116°32’~118°14’ E) and the Saihanba Mechanical Forest Farm (42°22′~42°31’ N,116°53’~117°31’ E), both located in Weichang County, Hebei Province, China [45]. The Mulan Forest Farm is characterized by hilly and mountainous terrain, with elevations ranging from 710 to 1975 m. This region experiences significant diurnal temperature fluctuation and is subject to a semi-arid to semi-humid continental monsoon climate [46]. The MAT is -1.4℃, and the MAP ranges from 400 to 500 mm, predominantly occurring from June to August. Both forest farms have the same major tree species, including L. principis-rupprechtii, B. platyphylla, Acer pictum, and Quercus mongolica. The Saihanba Mechanical Forest Farm has elevations between 1010 and 1940 m. The MAT is -1.3℃, and the MAP is 479 mm.

According to the ninth national forest inventory, the study area is predominantly covered by L. principis-rupprechtii, which occupies 38182.08 hm² (51.1% of the total forest area), followed by B. platyphylla, covering 16039.33 hm² (21.4%). The unit volume of L. principis-rupprechtii is the highest, at 172.5m³/hm², followed by B. platyphylla at 124.2 m³/hm². From July to September 2023, a total of 36 plots were investigated using a stratified random sampling approach (Figure 7). Each plot had similar site conditions. A total of 10 BP, 14 BL, and 12 LP plots were established, each covering an area of 30 m × 30 m. BP is natural secondary forest, LP is plantation forest, and BL is a mixed forest formed by replanting larch in a degraded birch forest. The main management practices applied were near-nature management, which aimed to optimize stand structure.

For each plot, trees with a diameter at breast height (DBH) exceeding 5 cm were measured. The recorded survey parameters involved trees species, DBH, tree height, crown width, and relative coordinates (X, Y). Table 1 shows the details of the sample plots.

4.2. Data Collection

4.2.1. Tree Core Sampling

For tree core sampling, trees within each plot were categorized into 2 cm diameter classes. The number of trees cores to be collected from each plot was determined based on the diameter distribution [46]. Two cores were carefully extracted from each tree at breast height (1.3 m), one in the eastern and one in the northern direction, using growth cones and preserving the bark as intact as possible. In total, 840 core samples were collected. These samples were then dried, fixed, and polished. Cross-dating was performed [46], and the widths of the annual rings were measured with an accuracy of 0.001mm using the WinDENDRO image analysis system [47]. The annual ring width series was analyzed using COFECHA 3.0 software [48].

4.2.2. Functional Traits

Leaves from each tree species were collected using pruning shears, and five indicators related to functional diversity were measured: leaf area (LA, cm²), specific leaf area (SLA, cm²/g), leaf dry matter content (LDMC, g/g), leaf nitrogen content (LNC, g/kg), and leaf phosphorus content (LPC, g/kg) [49]. From each sample tree, five mature, healthy, and undamaged leaves were selected [50]. The collected fresh leaves were promptly transported to the laboratory for analysis. Measurements of leaf traits were conducted following standardized protocols [51].

4.2.3. Soil Physicochemical Properties

Five soil collection points were set up at equal intervals along the diagonal of each sample plot. The surface litter was removed prior to excavating the soil profile. At a depth of 0-20 cm, soil samples were collected using a ring knife and aluminum box [52]. After thorough mixing, the dried samples were passed through a 20-mesh sieve for further analysis, including measurements of pH, SOC, total nitrogen (TN), total phosphorus (TP), total potassium (TK), available phosphorus (AP), available nitrogen (AN), and available potassium (AK). All soil traits were quantified following established methods [8,26].

4.3. Quantification of Forest Single Functions

4.3.1. Stand Productivity

Stand productivity (SP) characterizes a forest’s capacity to use energy and materials from the external environment, which is closely related to various forest functions [53]. Therefore, SP was selected as one of the main ecosystem functions for analysis in this study. Annual ring width growth was converted into basal area increments (BAI, cm²/year) using the following equation to quantify SP [54]:

$$BA{I}_i={\displaystyle \frac{\left(\pi {R_t}^2-\pi {R_{t-T}}^2\right)}{T}}$$

(1)

BAI_i represents the average annual BAI of the ith sampled tree over a period of T years(cm²/year), T is the length of study period, R_t and R_t-T denote the tree’s radius at t years and t-T years, respectively. To calculate the BAI, weighting is applied based on the trees’ diameter class, utilizing the ratio of DBHc to DBHss. DBHss refers to the average of the squared diameters within a specific diameter class, while DBHc is defined as the squared diameter of the central tree.

4.3.2. Functional Diversity

Functional traits are important indicators of how plants respond to environmental changes, as they are closely related to resource utilization and affect other ecosystem functions [55]. Functional diversity (FD), derived from these traits, plays a crucial role in ecosystem functioning. Communities with higher FD tend to have higher functional performance [56]. As a classical measure of FD, we analyze the following indicators:

Functional evenness (FEve)：This index assesses the evenness of species distribution within functional space (Eq. 2-4) [57]. A higher FEve indicates more uniform resource utilization across species.

$$FEve={\displaystyle \frac{{\displaystyle \sum _{i=1}^{S-1}\min \left(PE{W}_i\frac{1}{S-1}\right)-\frac{1}{S-1}}}{1-\frac{1}{S-1}}}$$

(2)

$$PE{W}_i={\displaystyle \frac{E{W}_i}{{\displaystyle \sum _{i=1}^{S-1}E{W}_i}}}$$

(3)

$$EWi={\displaystyle \frac{dist\left(i,j\right)}{w_i+w_j}}$$

(4)

EW_i is the evenness weight, dist(i, j) is the Euclidean distance between species i and j, w_i is the relative abundance of species i, and PEw_i refers to the the partial weighted uniformity of the branch.

Functional richness (FRic)：FRic reflects the niche that the community occupies and the extent of resources utilization [58]. A higher FRic indicates a greater resource utilization rate in the community. FRic is calculated using equation Eq. 5:

$$FRic={\displaystyle \frac{SFic}{Rc}}$$

(5)

In Eq. 5, SFic is the niche space occupied by species in community i, and Rc denotes the absolute range of trait c.

Functional dispersion (FDis): FDis calculates the greatest dispersion of functional abundance distribution within the community (Eq. 6-7). A higher FDis suggests greater trait disparity among species in the community [59].

$$FDis={\displaystyle \frac{{\displaystyle \sum _{i=1}^S(w_i}\times dis{t}_i)}{{\displaystyle \sum _{i=1}^S{w}_i}}}$$

(6)

$$c_k={\displaystyle \frac{{\displaystyle \sum _{i=1}^S\left(w_i\times x_{i,k}\right)}}{{\displaystyle \sum _{i=1}^S{w}_i}}}$$

(7)

S is the total number of species in the plot, w_i is the relative abundance of species i, dist_i is the distance from species i to the center of mass c_k, x_i,k is the kth trait value of species i.

Rao’s quadratic entropy (RaoQ)：RaoQ measures FD by a similarity matrix based on species functional attributes [60].

$$RaoQ={\displaystyle \sum _{i=1}^{S-1}{\displaystyle \sum _{i=i+1}^S{d_i}_j{p}_i{p}_j}}$$

(8)

S is species richness, d_ij is the overlap of the probability density function of two trait values, p_ip_j is calculated by multiplying the proportions of species i (p_i) and j (p_j).

4.3.3. Maintenance of Soil Fertility

SOC is an essential element of the carbon pools in terrestrial ecosystems. It is mainly derived from the decomposition of litter [61], and the secretion and activity of plant roots [62]. SOC content can reflect soil fertility and quality. The weighted average of SOC at a depth of 0-20 cm was used as an indicator to quantify the fertility maintenance function (SF) of forest soil [28].

4.3.4. Nutrient Cycling

The capacity for nutrient cycling (NC) in forest ecosystems, defined as the movement, transformation, and utilization of nutrients, revolves around the cycling processes between plants, litter, and soil [63]. Plant leaves capture carbon through photosynthesis, absorb nutrients such as nitrogen and phosphorus from the soil via their roots, and return these nutrients to the soil through litter decomposition [64]. Consequently, the contents of LNC, LPC, AN, AP, TN, and TP were utilized as quantitative indicators of nutrient cycling in this study [55]. These indicators represent vital nutrients for plant growth and are involved in the nutrient cycling process both above and below ground in forest ecosystem [65].

$$NC={\displaystyle \frac{{\displaystyle \sum _{i=1}^n{N}_{i,k}}}{n}}$$

(9)

Where n is the number of plots selected for different forest types, N_i,k is the kth nutrient index content of the ith plot.

4.4. Quantifying Ecosystem Multifunctionality

The single function approach, turnover approach, averaging approach, and threshold approach are widely utilized for quantifying EMF [1]. The single function approach is limited by its lack of qualitative information, whereas the turnover approach fails to offer a direct assessment of EMF [66]. The averaging approach tends to ignore the correlation between each index [67]. Therefore, this study adopts the multi-thresholds approach for calculating EMF. Following preliminary analysis, eight representative indicators were identified: SP, CWM_SLA, CWM_LDMC, SOC, TN, AN, TP, AP, LNC, and LPC. These indicators are essential for driving forest function, productivity, and soil nutrient cycling [38].

The multi-thresholds approach is expressed as the percentage of the maximum observed function value, representing the number of functions that exceed a set threshold [36]. In order to avoid the bias caused by over-weight assignment of some categories [68], all selected indicators were first standardized [36]. Subsequently, cluster analysis was performed to determine the optimal number of clusters, which was found to be 3 categories (Supplementary file, Figure S1). Based on this, EMF values were calculated for the 25%, 50%, 75%, and 90% thresholds using the weighting method, namely EMF_25, EMF_50, EMF_75, and EMF_90.

$$EMF={\displaystyle \sum _{i=1}^n\left(r_i\left(f_i\right)>t_i\right)}$$

(10)

Where n represents the number of ecosystem function indicators, f_i is the measured value of the ith indicator in the plot, r_i converts fi into a positive value, t_i is the thresholds.

4.5. Selection of Driving Factors of EMF

4.5.1. Biotic Factors

The standard deviation of basal area (Sdba) was used to represent stand structure properties [36]. Tree diversity (TD) was characterized by the overall weighted coefficient of variation (CV_D) in DBH, maximum tree height (MH), and SLA [8]. The equations listed below were used for the calculations:

$$C{V}_D＝{\displaystyle \frac{S{D}_D}{M{\mathrm{ean}}_D}}\times 100\%$$

(11)

$$CWM＝{\displaystyle \sum _{i=1}^{n}B{A}_i\times Trai{t}_i}$$

(12)

where SD_D and Mean_D represent the standard deviation and average value of DBH, respectively, and CWM refers to the community trait. Several metrics for species i: BA_i is the relative basal area, Trait_i is the trait value, and n is the number of species.

4.5.2. Abiotic Factors

The following factors were analyzed as representatives of the environment: soil nutrients (TN, TP, TK, AK, AN, and AP), pH, MAT, and MAP. The soil nutrients were first analyzed using principal component analysis (PCA) to obtain a set of nearly orthogonal variables [36], with the first principal component (SoilPC1) serving as the representative indicator (Supplementary file, Table S1). The contribution of climatic factors to tree growth is most significant during the growing season, which extends from June to August [69]. Consequently, we focused on the effects of MAT and MAP from June to August on EMF. MAT and MAP were estimated using Climate AP v3.10 [70], based on longitude, latitude, and elevation data.

4.6. Statistical Analysis

Before statistical analysis, all indicators were normalized to ensure data normality [36]. One-way ANOVA was conducted to evaluate differences in single functions and EMF among the three forest types [71]. The ‘FD’ package [72] was used to compute CWM, FEve, FDis, FRic, and RaoQ. The interactions between single functions and their driving factors were tested by Pearson’s correlation coefficient. The ‘vif’ function in the ‘car’ package was used to compute the variance inflation factor (VIF) [73] to address multicollinearity among variables. All VIF values were below 5, indicating that multicollinearity was not present. Finally, the ‘piecewiseSEM’ package was employed to investigate the factors of EMF and their interactions under thresholds of 25%, 50%, 75%, and 90% [74]. Fisher’s C, AIC, and p-value (p>0.05) suggested that the model provided a good fit [1]. All statistical analyses were performed using R 4.4.0 [75].

5. Conclusion

The study explored the influence mechanisms of stand type, biodiversity, and environmental factors on EMF under different thresholds. By quantifying the differences in single functions and EMF across the three stand types, we found that most individual functions in the mixed forest exceeded those in the pure forest. For example, SP was notably higher in the mixed forest compared to the corresponding pure forest. FD and SF in BL were considerably greater than in LP, although no significant difference was observed between BL and BP. Furthermore, the NC in BP surpassed that of LP, with no significant variation between LP and BL. At different thresholds, the EMF of BL was notably greater than that of BP and LP. Forest EMF was determined by both biotic and abiotic factors. At low thresholds (< 50%), the mass ratio and niche complementarity hypothesis jointly influenced EMF. As the threshold increased, the effects of these two effects on EMF weakened. At the 75% threshold, the mass ratio effect became the main driver of EMF. Environmental factors could indirectly affect EMF by altering functional traits, structural attributes, and tree species composition. The mixed species of L. principis-rupprechtii and B. platyphylla, with their differing functional traits, enhance the EMF’s capacity to resist climate change. Therefore, assessing the various functions and driving factors of larch forests in North China is crucial for biodiversity conservation and improving stand stability. In the context of multi-functional management of L. principis-rupprechtii forest, selecting birch species for mixing, increasing functional diversity, protecting the natural forest structure, and implementing near-natural management practices are key strategies to enhance EMF.