2.1 Study area, sites and forests
Our study area is located in the northeast of China, covering the Liaoning and Jilin Provinces (Fig. 1), which is classified as a temperate continental climate. The mean annual air temperate is 2.8°C, and the coldest and warmest monthly mean are − 13.7°C in January and 19.7°C in July (Hao et al. 2007). The mean annual participation is 700 mm which mostly occurs during the growing seasons (from June to September). Broad-leaved Korean pine (P. koraiensis) mixed forest, well-known for their high species diversity, productivity and complex stand structure, was the dominant vegetation type in the study region. The predominant soil type of the studied area is Eutric cambisols according to the FAO soil classification system (Hao et al. 2007). Some areas had been subjected to variable intensities of human disturbances, but there were no severe human disturbances since 1998 when the forest protection policy was strictly implemented (Dai et al. 2004). Thus, forests recovering from disturbances include stands with different successional stages in the study area (Chen et al. 2014).
To quantify the effects of abiotic factors, disturbance intensity, plant trait diversity and composition on soil water storage capacity, we established 11 permanent forest sites during 2012 ~ 2013 (Table 1), and then re-investigated those forest plots at an interval of five years using the standard field approaches (Yuan et al. 2018). Within each site, all trees with a diameter at breast height (DBH) > 1 cm were marked, measured, recognized to species level, and positioned in contiguous 20×20 subplots (Hao et al. 2007). In this way, our data covered 22,766 stems in total belonging to 81 species, 46 genera, and 26 families across 260 subplots (Yuan et al. 2018). In addition, the topography of each subplot was measured by assessing the elevation of four corners of each subplot using Electronic Distance Measuring Device (Hao et al. 2007), and then the mean elevation, convexity and slope were calculated for each subplot (Harms et al. 2001). Elevation of the studied plots ranged from 640.4 to 1023.1 m.
Table 1. The descriptive summary of the studied sites and plots.
Site
|
Site size (ha) (dimension, m)
|
No. of subplots
|
Latitude longitude
|
Elevation (m)a
|
Species richness
|
DBH (cm)a
|
|
Plots of low-disturbance level
|
|
|
|
L1
|
1 (100*100)
|
25
|
127.98N; 42.35E
|
877.5
(875.1; 879.7)
|
25
|
8.68
(1.0, 115.2)
|
|
L2
|
1 (100*100)
|
25
|
126.23N; 43.18 E
|
727.5
(724.2; 731.1)
|
35
|
8.98
(1.0, 152)
|
|
L3
|
1 (100*100)
|
25
|
127.88N; 42.25E
|
998.8
(995.6; 1002.1)
|
38
|
7.56
(1.0, 110.0)
|
L4
|
1 (100*100)
|
25
|
127.91N; 42.21E
|
1107.2
(1105.6; 1108.3)
|
21
|
8.64
(1.0, 96)
|
Plots of medium-disturbance level
|
|
|
|
M1
|
1 (100*100)
|
25
|
127.86N; 42.48E
|
1010.9
(1016.3; 1023.1)
|
34
|
7.83
(1.0, 96.5)
|
M2
|
1 (100*100)
|
25
|
127.94N; 44.01E
|
721.6
(698.4;743.7)
|
41
|
8.98
(1.0, 75.0)
|
M3
|
0.8 (80*100)
|
20
|
124.79N; 40.91E
|
834.1
(817.0; 851.0)
|
40
|
7.22
(1.0, 60.9)
|
Plots of high-disturbance level
|
|
|
|
H1
|
1 (100*100)
|
25
|
126.48N; 42.36E
|
758.6
(749.6; 764.8)
|
47
|
3.4
(1.0, 75.0)
|
H2
|
1 (100*100)
|
25
|
128.17N; 42.19E
|
652.9
(640.4; 666.2)
|
45
|
6.35
(1.0, 68.5)
|
H3
|
1 (100*100)
|
25
|
130.16N; 43.39E
|
717.7
(705.6; 726.2)
|
34
|
6.3
(1.0, 70.0)
|
H4
|
0.6 (60*100)
|
15
|
124.90N; 41.33E
|
892.2
(873.1; 909.4)
|
43
|
7.15
(1.0, 53.8)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
a Mean value and range (min, max) of each 20*20m subplots. DBH, diameter at breast height.
2.3 Quantification of plant functional trait diversity and composition indices
For quantifying the multiple facets of plant diversity, six functional traits were measured which were closely associated with forest growth, recruitment and death (Yuan et al. 2019), including maximum height (MH), leaf area (LA), specific leaf area (SLA), leaf carbon content (LCC), leaf nitrogen content (LNC), and leaf phosphorus content (LPC). The detailed measurement approaches for these six functional traits are described in Yuan et al. (2016). Using the data of six functional traits, we measured three multidimensional functional trait diversity indices, namely, functional richness (FRic), functional evenness (FEve), and functional dispersion (FDis). FRic represents the multivariate traits space that species in the community occupied, FEve indicates the regularity of community traits distribution, and FDis is the mean distance of the species trait to the centroid of all species weighted by basal area (Laliberté and Legendre 2010; Villeger et al. 2008). Besides, we also calculated the community-weighted mean of each trait (CWMMH, CWMLNC, CWMLCC, CWMLPC, CWMLA, and CWMSLA) within each plot, weighted by the species' relative basal area, to represent the plant functional trait composition (Ali et al. 2017; Yuan et al. 2019). The CWM and FD indices were calculated by FD package in R 3.5.3 (Laliberté and Legendre 2010).
2.4 Quantification of soil water storage capacity and soil organic carbon
In 2018, we randomly selected three soil sampling sites in the midpoints between the central point and four corners in each 20 m×20 m subplot. In each sampling site, two soil corers using stainless cylinders of 100 cm3 in volume were selected for the bulk soil density and capillary water storage contents measurement after removing large debris. The corers containing large roots were abandoned for the precise analyses of data. Subsequently, five soil cores (3.8 cm in diameter, 10 cm deep) at each sampling point were collected, pooled, and transferred to the laboratory with plastic zipper bags for further chemical analyses. Each soil sample was further divided into two parts, i.e., one for soil organic carbon analysis using the dichromate oxidation method (Lu 1999), and another one for soil moisture measurement after 12 hr dried at 105°C.
Soil water storage capacity was measured through soil porosity, and this could be divided into capillary porosity, non-capillary porosity and total porosity corresponding to soil capillary storage (CW), soil non-capillary storage (NCW), and soil saturated water storage (TSW) (Chen et al. 2019; Xia et al. 2017). The CW represents the soil water content in the soil water retention curve when the pressure is -33 kPa whereas the NCW reflects the water contents difference between 0 ~ -33 kPa (Ahuja et al. 1993; Liu et al. 2009). Moreover, these two parts also indicate soil water retention and infiltration capacity (Mo et al. 2011; Rabot et al. 2018). TSW, CW, and NCW were calculated using Eq. (1), Eq. (2) and Eq. (3) (Chen et al. 2019; Wu et al. 2016; Zhang et al. 2010):
Where TP is the total soil porosity (%) measured as TP=(1-BD/ds)×100, BD is the soil bulk density (g cm3), ds is the soil particle density (2.65 g cm-3); CP is the soil capillary porosity (%) measured as CP = (WC/V)×100, Wc is the soil capillary water contents (g cm3) (Liu et al. 2009), V is the volume of the soil core (cm3); NCP is the soil non-capillary porosity (%) measured as NCP = (TP – CP)×100; h is the height of soil top layer (0.2 m), TSW is the total soil water storage content (t hm-2), NCW is the non-capillary storage content (t hm-2), and CW is the capillary storage content (t hm-2).
Total soil porosity (TP) was measured based on the measured soil bulk density while assuming that soil particle density is 2.65g cm-3, WC was additional water weight after placing the stainless cylinders soil core in a tray with a 5-mm level of water until filter paper at the top of each core became moist (Liu et al. 2009), whereas non-capillary was the difference between TP and CP (Eq. (3)) (Wu et al. 2016).
2.5 Statistical analyses
We assessed the effects of disturbance on soil water storage capacity (i.e. TSW, NCW and CW) using a two-way ANOVA with Tukey’s test as a post hoc analysis to assess the significant differences among disturbance levels. As such, we also tested the differences for associated variables, which may influence soil water storage capacity (i.e. above-ground and below-ground variables), among disturbance intensity levels.
Before testing the conceptual model in Fig. S1, we assessed the spatial autocorrelation in the response variables (TSW, NCW and CW) among subplots using the generalized least-square models (GLS) with and without spherical autocorrelation. Our GLS analysis indicated that the models without spherical autocorrelation structures usually showed the lower Akaike Information Criterion (AIC) values compared to spherical autocorrelation models (Table S2, S3 and S4), suggesting that there was no strong spatial autocorrelation among subplots.
Based on Pearson’s correlation (Table S5), the best combination of variables (Tables S6 and S7), and the conceptual model (Fig. S1), we identified elevation, CWM of traits, functional trait diversity, and soil organic carbon as the best factors influencing soil water storage capacity (i.e. CW, NCW, and TSW). To test the conceptual model in Fig. S1, we used structural equation modeling (SEM) because it allows us to test the multiple theories, direct and indirect effects. To critically evaluate the best-fitted SEM (Hoyle 2012), we used the highest Bentler’s Comparative Fit Index (CFI > 0.90), the root mean square error of approximation (RMSEA ≤ 0.05), chi-square (X2) test (P-value > 0.05), standardized root mean square residual (SRMR ≤ 0.05) and lowest AIC value. To compare the mass ratio hypothesis and niche complementarity effects on soil water storage capacity, we selected CWMSLA and FEve as the representative of plant functional trait composition and diversity, respectively, in SEMs. The direct, indirect, and total effects of predictors on response variables were calculated. The relative importance of each predictor was calculated as the percent of the given predictor standardized coefficient to the sum of standardized coefficients of all predictors (Yuan et al. 2019).
To complement the results from SEMs, we also conducted partial regressions and simple bivariate regressions (Grace et al. 2016). The SEM analysis was conducted using the “lavaan” package (Rosseel 2012). All analyses were conducted in R. 4.0.2 (Team RDC, 2019). All predictors in our research were standardized to a mean of 0 and standard deviation of 1, and the response factors (i.e. CW, NCW and TSW) were natural-log transformed before further analyses. Summary of variables used in the analysis is provided in Table S1.