Plant litter decomposition is a key process for ecosystem carbon and nutrient cycling. Growing evidence suggests that substantial variation in litter decomposition occurs at a fine scale, but the contributions of microorganism- and fauna-driven decomposition to this variation, and the relative control of biotic and abiotic drivers over this variation, remain virtually unexplored. To address this knowledge gap, we used a spatially explicit network of 113 evenly spread plots in a 9-ha subtropical coniferous forest to evaluate the variation and controls of microorganism- and fauna-driven decomposition at a scale where macroclimate, dominant vegetation, and litter quality are kept constant. Despite keeping dominant decomposition drivers constant, the variation in decomposition was larger than that commonly reported in regional studies and amounted to ca. a third of the variation previously reported at the global scale. Furthermore, while abiotic factors, including topographic conditions and soil fertility, could explain 31% of the variation in microorganism-driven decomposition, they could not explain any variation in the fauna-driven decomposition, suggesting contrasting sensitivities of microorganism- and fauna-driven decomposition activity to microenvironmental factors.

Synthesis: Our study highlights the need to consider the local-scale variation in litter decomposition rate using a spatially explicit approach, in order to better identify the factors driving the microbial and faunal pathways of organic matter turnover.

Study area
We conducted this experiment in a subtropical forest located in the Yingzuijie National Nature Reserve (26°46' N - 26°59' N, 109°49' E - 109°58' E, 270 - 938 m a.s.l), within Huitong County, Hunan province, China. All fieldwork was conducted with the authorization of the Hunan Huitong National Research Station of Forest Ecosystem. The site has a humid mid-subtropical monsoon climate characterised by a mean annual temperature of 16.5 °C and a mean annual precipitation of 1,200 mm. The soil is classified as Anthrosols and Alisols (IUSS Working Group WRB, 2022). The stand structure corresponds to a forest plantation with Chinese fir as the dominant tree species (30-year-old stems and 29 m2 ha-1 basal area) with an initial planting density of 2 m × 2.3 m after clear-cut and burning. During the stand development, native woody species regenerated, including Camellia oleifera, Itea omeiensis, and Loropetalum chinense. We focused our study on a 9-ha (300 m × 300 m) area, which is divided into 225 plots of 20 m × 20 m each. To evaluate the effect of local-scale variation in environmental conditions and litter decomposition, we selected 113 non-adjacent plots from the 225 plots for our study (Fig. 1).

Litter decomposition experiment
To compare the response of litter decomposition, driven by microorganisms and soil fauna, to environmental variation at local scale, we quantified the decomposition rates of Chinese fir leaf litter by using litterbags with fine (0.1 mm) and large (2.0 mm) mesh sizes in all 113 plots and 10 analytical replicates, leading to 2260 litterbags (113 plots × 2 mesh sizes × 10 analytical replicates). To obtain the needle litter, we collected freshly fallen litter from trees in the immediate surroundings of the study site. The collected litter was air-dried at room temperature (25 ℃), and only complete needles were used. Variation in chemical properties of the Chinese fir litter used in this study was minimal (see detailed initial chemical determination and properties described in Yin et al. (2023). We prepared 15 cm × 15 cm litterbags with nylon mesh to allow the access of microorganisms only, we used 0.1 mm mesh for both sides; to allow the access of microorganisms and micro- and mesofauna, we used a 2 mm mesh for the upper side to allow soil animal access, and a 0.1 mm mesh for the side facing the soil surface to prevent losses of litter fragments. We then filled the litterbags with 5 g (± 1 mg) of air-dried litter. To convert air-dry mass into oven-dry mass for initial litter, we weighed initial air-dry litter subsamples, dried them at 60°C for 48 h, and reweighed them to obtain a conversion factor. In September 2019, ten pairs of litter bags, each pair consisting of two litter bags with different mesh, were randomly placed in each 20 m × 20 m plot on top of the forest floor and incubated for 383 days in situ. To obtain a robust estimate of microorganism- and fauna-driven decomposition for each 20 × 20 m plot, the ten pairs were placed in ten subplots, which were randomly selected out of the 16 possible 5 × 5 m subplots. Upon sampling, the remaining litter was taken out of the litter bags and gently rinsed with tap water to remove soil particles. 
The samples were dried at 60 °C to constant mass and weighed. We used percentage mass loss (ML) as a measure of decomposition rate, which we calculated with the following equation:
ML=(M_0-M_t)/M_0×100    (1)
where Mt is the remaining mass, M0 is the initial mass of the litter after correction for moisture (i.e., equivalent mass at 60 °C). 
We calculated the ML for the 0.1- and 2-mm bags to represent the decomposition rate driven by microorganisms (MLM) and the decomposition rate driven by both fauna and microorganisms (MLF&M), respectively. The difference (MLF&M - MLM) between them was calculated as the fauna-driven decomposition (MLF) (Coq et al., 2010). To overcome the non-independence of within-plot replicates and to avoid pseudo-replication, we used the mean values of the ML from the 10 analytical replicates in each plot for statistical analyses.
To compare our range of variability in decomposition rates with that observed in other studies, we also calculated the litter decomposition rate constant (k) with the following equation: 
k = -ln⁡〖(M_t/M_0 )/ t〗      (2)
where Mt is the remaining mass, M0 is the initial mass of the litter after correction for moisture (i.e., equivalent mass at 60 °C), and t is the number of days. Similar to ML, we calculated the k for the 0.1- and 2-mm bags to represent the decomposition rate driven by microorganisms (kM) and the decomposition rate driven by both fauna and microorganisms (kF&M), respectively. The difference (kF&M – kM) between them was calculated as the fauna-driven decomposition (kF). 

Topographic condition measurements and calculations
To characterise the topographic conditions of each plot, we generated an elevation raster with a pixel size of 20 × 20 m corresponding to each plot. This was obtained by conducting kriging interpolation to the elevation data collected at the four corners of each plot (Burrough, 1986; McBratney & Webster, 1986). Slope, aspect, and topographic wetness index (TWI) were calculated in ArcGIS 10.8 (ESRI, 2018) based on the elevation raster data. To obtain a meaningful measure of aspect, we transformed the direction of downhill slope faces (aspect0) in compass degrees (0-360°) into a variable ranging from 0 to 1 (the higher the value, the more southward the aspect) with the following equation:
Aspect=sin(π×[aspect]_0/360)    (3)
TWI represents the spatial distribution of surface saturation and potential water accumulation in plots and has been widely used to estimate spatial variability of local soil hydrological and physical properties (Seibert et al., 2007; Xin et al., 2016). TWI was formulated as follows:
TWI = ln(a/tanβ)    (4)
where a is the upslope contributing area per unit contour length (or Specific Catchment Area, SCA) and tanβ is the local slope gradient for estimating a hydraulic gradient. The estimates of a and β were calculated following the method proposed by Qin et al. (2011), which has been shown to reflect local terrain and drainage conditions well. 

Soil sampling and analysis
To determine how well soil chemical and microbial characteristics predicted decomposition, we collected one soil core (0-10 cm in the topsoil layer) using an auger (2 cm in diameter) in May 2019 in each of the 16 possible 5 × 5 m subplots within each 20 × 20 m plot. This led to 16 cores, which we then pooled and homogenised to form a composite sample. We removed visible debris and stones from the soil samples, sieved these soil samples through a 2-mm mesh, and divided them into three subsamples. One subsample was stored at 4 °C to determine the contents of available nutrients, another was air-dried at room temperature to determine the total contents of nutrients, and the remaining one was freeze-dried and stored at -80 °C for microbial determination. 
To determine soil fertility, we followed standard methods described by Liu et al. (1996). Briefly, we measured soil organic carbon (SOC) and soil total nitrogen (TN) using an element analyser (Vario Macro Elementar, GmbH, Germany). Total phosphorus content (TP) was determined by inductively coupled plasma mass spectrometry (ICP-Q, Thermo Scientific, America) after digestion with a mixed solution of HNO3: HF (with a volume ratio of 5:1). Soil carbon-to-nitrogen (C:N) and carbon-to-phosphorus (C:P) ratios were calculated for each plot. Soil available nitrogen (AN) was extracted using a 2 mol L-1 KCl solution, and soil available phosphorus (AP) was extracted by a 0.03 mol L-1 NH4F-0.025 mol L-1 HCl solution, then determined by a continuous-flow analyser (AA3, SEAL Analytical, Australia). 
To determine soil microbial communities, soil microbial DNA was extracted from 0.3 g of sample using the MoBio PowerSoil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, United States). The bacterial 16S rRNA gene and fungal ITS2 region were amplified using the primers 515F/909R (Tamaki et al., 2011) and ITS4/gITS7F (Ihrmark et al., 2012), respectively. Sequencing libraries were generated using TruSeq® DNA PCR-Free Sample Preparation Kits and sequenced on an Illumina MiSeq platform with the 2 × 250 bp V2 Reagent Kits. We then used the Quantitative Insights Into Microbial Ecology (QIIME, V1.9.0) pipeline to process DNA sequence data. Bacterial and fungal sequences were independently clustered into operational taxonomic units (OTUs) at a 97% identity threshold using UPARSE software (Edgar, 2013). The results of OTU abundance tables were rarefied to an even number of sequences per sample (7108 and 24325 sequences per sample for bacteria and fungi, respectively), corresponding to the minimum number of sequences for a single soil sample (details in Zhai et al., 2023). All the sequencing values were deposited in the National Center for Biotechnology Information (NCBI database) with the BioProject accession no. PRJNA752698 (https://www.ncbi.nlm.nih.gov/sra/PRJNA752698). 

Vegetation survey and plant attribute calculation
To evaluate the effect of vegetation attributes on decomposition, we estimated, for each plot the basal area (BA) using diameter at breast height (DBH) of all woody plants as a proxy of vegetation biomass, species richness (SR) as a proxy of diversity, and the coefficient of variation of DBH (CVD) as a proxy of forest structural diversity (Ali, 2019). Since Chinese fir (CF) and regenerated woody plants (RWP) species may affect decomposition differently due to their different biomass, diversity, and position in the forest (Rago et al., 2021), we calculated the basal area and CVD for Chinese fir and RWP separately, and considered the RWP species richness, leading to five indicators (BACF, CVDCF, BARWP, CVDRWP, and SRRWP). 
To do so, we measured the diameter at breast height (DBH) ≥ 1 cm of all woody plants in each plot following a standard field protocol (Condit, 1998), between September 2019 and February 2020. We calculated the basal area using the following equation:
BA=∑_0^m∑_0^n[(π×[(DBH]_n]_ /2)]^2 )/(A)  ]    (4)
where BA is the relative basal area (m2/ha), m is the species number of RWPs (or 1 when calculating the basal area of Chinese fir) in the specific plot, n is the number of trees of the species m in each plot, 〖DBH〗_n is the DBH (m) of the tree n, and A is the area of each plot (0.04 ha). 

Data analysis
We conducted a linear mixed effects model using the ‘lmer’ function from the ‘lme4’ (Bates et al., 2015) package to compare the contribution of different groups of decomposers (i.e., fauna plus microorganisms, microorganisms, and fauna) on decomposition rate (ML and k), with ‘plot’ as a random effect. The ‘anova’ and ‘ranova’ functions from the ‘lmerTest’ package (Kuznetsova et al., 2017) were used to test the significance of fixed and random effects. A Tukey HSD test was used to determine significant differences between MLF&M, MLM, and MLF, and between kF&M, kM, and kF across all plots. In addition, we used Pearson’s correlation analysis to evaluate relationships between MLF&M, MLM, and MLF, and between kF&M, kM, and kF. We used the ‘PCA’ function from the ‘FactoMineR’ package (Lê et al., 2008) to visualise how different groups of environmental predictors were related (including 4 topographic conditions indicators, 4 soil fertility indicators, and 5 plant attributes). For each principal component analysis (PCA), component variables were centred and standardised before ordination. In addition, we used a PCA to visualise soil bacterial and fungal community composition based on the OTUs composition matrix in each plot, using Aitchison distance after CLR (centred logratio transform)-matrix transformation with the “vegan” and “compositions” packages (Boogaart et al., 2024; Oksanen et al., 2022). 
To determine the relative importance of topographic conditions, soil fertility, plant attributes, fungal community, and bacterial community composition on MLF&M, MLM, and MLF, separately, we used individual and multiple regression models to determine how much variation of each response variable could be explained by these explanatory variables. Before multiple regression model analyses, each PC1 of topographic condition (PC1TC), soil fertility (PC1SF), plant attribute (PC1PA), fungal community (PC1FC), and bacterial community (PC1BC) were standardised by the ‘scale’ function to eliminate unit dimension effects. Full regression models were executed using the 'lm' function with the R syntax: ML ~ PC1TC + PC1SF + PC1PA + PC1FC + PC1BC. Model selection was performed with the 'dredge' function of the 'MuMIn' package (Bartoń, 2022), which ranks all candidate models (i.e., all combinations of explanatory variables included in the full model) based on the lowest corrected Akaike information criterion (AICc). We then estimated the relative importance (Σwi) of each predictor by summing AIC weights in models that included the explanatory variable, following Burnham and Collins (2015). Last, we reported the R2 of the most parsimonious models. We ran these analyses similarly on kF&M, kM, and kF. All data analyses and visualization were performed using R, v.4.1.2 (R Core Team, 2021).