Linking Vegetation, Soil Carbon Stocks, and Earthworms in Upland Coniferous–Broadleaf Forests

: Linking vegetation, soil biota, and soil carbon stocks in forests has a high predictive value. The speciﬁc aim of this study was to identify the relationships between vegetation, earthworms, and soil carbon stocks in nine types of forests dominating autonomous landscape positions in a coniferous–broadleaf forest zone of the European part of Russia. Mountain forests were selected in the Northwest Caucasus, while plain forests were selected in Bryansk Polesie and on the Moskva-Oka plain. One-way analysis of variance (ANOVA) and v -tests were used to assess the impact of different factors on soil C stocks. To assess the contribution of vegetation, litter quality, and earthworms to variation of carbon stocks in organic (FH-layer) and mineral layer (0–50 cm), the method of hierarchical partitioning was performed. The highest C stocks in the organic horizons were associated with the low-quality litter, i.e., with a low base saturation, high acidity, and wide C/N ratio. The highest soil C stocks in the mineral layers were found in mixed forests with the highest richness of plant species, producing litterfall of different quality. The C stock in the organic horizon was negatively related to the biomass of worms that process the litter, while the carbon stock in the mineral layers was positively related to the biomass of worms whose life activity is related to the mineral layers. These ﬁndings demonstrated the substantial inﬂuence of plants producing a litter of different quality, and of earthworms, belonging to different functional groups, on soil C stocks in coniferous–broadleaf forests.


Introduction
Forests play a huge role in climate regulation, including through their ability to absorb greenhouse gases and store carbon, both in phytomass and in soils. Soils contain 40% or more of the site carbon stock [1]. Soil carbon stocks depend on abiotic factors of soil formation, such as climate and soil-forming rocks. It has been shown that soil carbon stocks are positively correlated to average annual temperature, annual precipitation, and, therefore, to net primary productivity [2]. It has also been demonstrated that the levels of carbon accumulation in soils depend on particle size distribution in the soil-forming rocks and their chemical composition, as well as on forest management regimes and forest disturbances, first of all, fires [3]. However, soil carbon stock also depends on biotic factors, such as vegetation and soil biota. It has been shown that in forests of different types, variation in soil carbon stocks was related to the species composition of tree plants [1] and ground vegetation [4]. However, the question arises, what is the combined effect of different woody and herb plant species on carbon stocks? Although it was not possible to identify significant differences between carbon accumulation in soils under broadleaf trees in monocultures in common garden experiments [5], it was shown that a greater variety of tree species leads to an increase in organic carbon stocks in the soil [6]. Having looked into the influence of the richness of boreal forest tree species on their productivity and carbon dynamics, Shanin et. al [7] showed that mixed forests are more productive than monodominant ones.
Vegetation composition determines the amount, quality, and rate of decomposition of plant litterfall, its horizontal soil distribution, and the distribution of carbon compounds within the soil profile [8]. Plants that produce low-quality, slow-decomposing litterfall with a high content of secondary metabolites (lignin, tannins, etc.) and low content of nutrients (nitrogen, calcium, magnesium, etc.), contribute to the formation of a thick forest soil organic horizon since such litterfall is slowly processed by soil biota [9][10][11]. This litter is mainly decomposed by microorganisms, saprotrophic micro-and mesofauna, living in the litter, which could lead to increased carbon dioxide emissions from the soil surface and losses of soil carbon [12]. Plants that form high-quality litterfall with a high content of nutrients and low content of secondary metabolites contribute to a more active decomposition of plant residues by soil biota [13]. Such litterfall is actively processed by the soil biota, including macrofauna, especially earthworms, functioning in different soil layers. Their activity contributes to the fixation of carbon in soil due to processes of humus formation and the transfer of soil organic matter to deeper horizons as a result of digging [12]. Through their intestinal enzymes, earthworms also activate a pool of soil microorganisms [14][15][16]. Various combinations of plant litter quality and composition of soil macrosaprophages can affect the balance between carbon sequestration in soil, its emission, and removal with soil waters. An experiment with litterfall of 16 common European tree species showed that the litter mass loss was significantly higher in the presence of earthworms than in that of isopods. The effect of functional diversity of the litter material was highest in the presence of both macrosaprofages' groups, supporting the concept that litter diversity effects are most pronounced in the presence of different detritivore species [17]. Earthworms' trophic activity leads to a decrease in litter amount, but due to humus formation and vertical transfer of organic matter by earthworms, carbon accumulates in soil mineral horizons [18]. The findings of other, earlier studies demonstrated the opposite effect of earthworms on carbon stock in soil or the absence of such an effect [19,20]. It is important to note that the studies did not always take into account the belonging of worms to different functional groups, which could be the reason for different conclusions. However, the belonging of earthworms to individual groups functioning in different soil horizons and having different effects on the soil's organic matter is of great importance [21][22][23]. To give an unambiguous explanation of how earthworms affect carbon accumulation in soil, it is necessary to take into account both their combined influence and the influence of individual functional groups.
The aim of this study is to investigate the relationships between vegetation, soil carbon stocks, and earthworms in forests dominating autonomous landscape positions in a coniferous-broadleaf forest zone of the European part of Russia.
The main scientific questions are: (1) What characteristics of vegetation and earthworms are the most informative predictors of variation in soil carbon stocks in coniferous-broadleaf forests?
(2) Does higher plant and earthworm functional diversity contribute to higher soil carbon stocks in coniferous-broadleaf forests?

Study Areas
The plots for the study were established in coniferous-broadleaf forests at autonomous positions of landscapes in 2016-2019. Mountain forests were selected in the Northwest Caucasus and plain forests on the Moskva-Oka plain and in Bryansk Polesie ( Figure 1). In each region, forests belonging to the types dominating autonomous landscapes were selected. In each type of forest, 3 plots with a size of 0.25 ha were established, in total 27 sample plots. Northwest Caucasus (NC). The study was carried out in coniferous-broadleaf forests in the upper reaches of the Pshekha River and the Belaya River (43 • 59 -44 • 03 N, 39 • 42 -40 • 08 E). The altitude is 650-700 m above sea level. The average annual temperature is 10.2-11.1 • C, the average precipitation is 1080-1145 mm [24].
The soil cover is dominated by Dystric Cambisols [27], with varying degrees of gleying on the eluvium of clayey shales. The signs of gleying are more pronounced in forests dominated by beech, which has a relatively low transpiration rate [28] [29]. The average annual temperature is 4.6 • C, the annual precipitation is 675 mm [24]. A significant area is occupied by the secondary forests, developing after logging.
Forests are dominated by spruce (Picea abies (L.) H. Karst) and lime (Tilia cordata Mill.) [30,31]. At the stage of selecting the objects of study, forests with a predominance of lime (3 plots) and a predominance of birch (3 plots) were identified. Further studies did not reveal any differences in the characteristics of vegetation and soils between these types of forests, and therefore, when statistical processing of the data, these forests were combined into one type-lime-birch (MO1, 6 plots). Spruce-broadleaf (MO2, 3 plots) forests were also selected as the objects of the study. According to the Braun-Blanquet approach, the MO2 forest communities belong to the association Rhodobryo rosei-Piceetum abietis Korotkov 1986; and the MO1 forest communities belong to the Mercurialo perennis-Quercetum roboris Bulokhov [32,33]. The soil cover is dominated by Albic Retisols [27].
The modern forest cover of BP is represented by secondary communities, and among them, pine forests (BP1, 3 plots) and pine-broadleaf (BP2, 3 plots) forests are widely distributed in the autonomous positions of landscape in the study area [34][35][36] The detailed physico-geographical, geobotanical characteristics of all objects of study, as well as a description of macrosaprophages and the physicochemical characteristics of soils, have been also presented earlier [23,[36][37][38].

Methods
Vegetation study. Geobotanical description of the forest communities was made on 20 × 20 m square plots within each permanent sample plot. In each region, 11 to 14 geobotanical descriptions were made for each type of forest, 100 descriptions were completed in total. A complete floristic list was compiled for all vegetation plots, taking into account the layered structure of the forest vegetation and the abundance of species. The coverage-abundance score was assessed using the Braun-Blanquet scale. Latin names of vascular plants and mosses were given according to the World Flora Online [39].
Plant species of all layers were divided into two functional groups taking into account litterfall quality: fast and slowly decomposed. As for the herb layer, the first functional group (woody plants) includes all woody species, namely trees, shrubs, and dwarf shrubs, and the second functional group (herb plants) includes all other species.
To estimate stem volume in each sample plot, the diameter and height of each tree higher than 1.5 m were measured. Stem volume was calculated using volume tables [40,41]. Then, carbon in stem biomass was calculated using the equation in accordance with the officially accepted national guidelines [42]: where CP i is carbon in stem biomass of the j-th tree species, t C ha −1 ; V i is stem volume, m 3 ha −1 ; KP i is conversion factors for converting stem volume to C, t C m −3 Macrofauna study. To count soil macrosaprophages, soil samples with an area of 25 × 25 cm and a depth of 30 cm were excavated and hand-sorted [43]. Fieldwork was carried out in the spring and summer seasons of 2016-2019. In each region, 30 to 48 soil samples were collected to count the earthworms in each forest type. A total of 438 soil samples were analyzed.
Earthworms were placed in 95% ethanol. Millipedes, insect larvae, and mollusks were placed in 70% ethanol. The species identification of earthworms was carried out according to the field guide of T.S. Vsevolodova-Perel [44]. The functional groups of earthworms are given by classification according to Vsevolodova-Perel [44]: epigeic, epi-endogeic, endogeic, and anecic.
Insects, mollusks, crustaceans, and millipedes were identified up to supraspecific taxa (using field guides by [45][46][47]). The number of macrosaprophages was counted and their biomass was measured: invertebrates with a full intestine placed in ethanol were weighed. To assess the functional diversity of macrosaprophages, the Simpson's diversity index (D) was calculated: where N is total biomass of functional groups, g; and Ni is biomass of functional groups of the i-th groups, g. Soil study. In each of 27 sample plots in the nodes of the regular network with a 10 m pitch, individual samples were taken from the soil organic and mineral horizons of the soil up to a depth of 50 cm using a soil drill in 2016. In addition, the soil samples were also taken in each type of forest in 2019. The samples of the organic horizon were taken using a 0.25 × 0.25 m frame, and they were divided into two layers: L-layer (litter layer) and FH-layer. The individual samples were combined into 3 samples in accordance with the soil horizons (A, AE/E, B) in each plot. From each horizon, 157 samples were analyzed, resulting in a total of 785 samples. In addition, on each sample plot, samples were taken in one soil profile 100 cm deep to characterize the soil-forming rock.
In the laboratory, soil samples were dried and sieved through a 2-mm sieve. The carbon (C) and nitrogen (N) content was determined in all samples using a CHN analyzer (EA 1110 (CHNS-O)). The actual acidity (pH) in the water extract from the L-layer samples was determined using a potentiometer. Exchangeable acidity in the L-layer samples was determined by extraction with 1N KCl followed by titration (pH = 7.0) [48]. Bio-available (ammonium acetate extraction, pH 4.65, [49]) metals (Ca, K, Mg, Na) in the L-layer samples were determined by atomic absorption spectrometry (AAnalyst 800 spectrometer). Particlesize distribution in soil-forming rock was determined according to national methodical guidelines [50].
The bulk density and thickness of mineral horizons were measured. In the laboratory, all samples were dried at 105 • C and weighed. The C stock in the organic horizon (FH-layer) was calculated by multiplying the sample weight by the C content. The C stock (Cst) in soil mineral profiles was calculated by multiplying bulk density, carbon content, and mineral layer thickness [42]: where BD is bulk density, g/cm 3 ; C is carbon content, g/kg; and H is layer thickness, cm.
Since the studied soils were characterized by a sharp vertical differentiation of soil profile into genetic horizons, the calculations of carbon stock for fixed mineral layers (0-30 and 30-50 cm layers, and the sum of them) were based on its estimation in the horizons, according to officially accepted national guidelines [42].
Statistical analysis. Statistical processing was performed using the R statistical software [51].
One-way analysis of variance (ANOVA) and v-tests [52] were used to assess the impact of different factors on soil C stocks. The tests and the associated descriptive statistics were calculated using the function catdes of the FactoMineR package [53].
Altogether 187 species on 98 plots were used to determine the ordination pattern of vegetation. Non-metric multidimensional scaling (NMDS) was performed on logtransformed species abundances (percent of coverage) using the metaMDS function of the vegan package of the R statistical environment.
To assess the contribution of vegetation, litter quality, and earthworms to variation of carbon stocks in organic (FH-layer) and mineral layer (0-50 cm), the method of hierarchical partitioning was performed using the hier.part package [54]. The C stocks in the organic (FH-layer) and mineral layer (0-50 cm) were regarded as the dependent variables. The following were regarded as predictors: (1) stem wood stock, the projective cover of woody plants in the herb layer, the projective cover of herb plants in the herb layer, the projective cover of shrub layer and tree layer; (2) species richness (SR) of woody plants in the herb layer, SR of herb plants in the herb layer, SR in the shrub layer, SR in the tree layer, and total SR of shrub and tree layers; (3) N content, C/N ratio, pH values, base saturation (BS), Ca, K, Mg content in the L-layer; and (4) the contribution of biomass of worms associated with the organic horizon (epigeic) and with the mineral layers (endogeic) to variation of the soil C stock was assessed. The contribution of biomass of worms belonging to epi-endogeic and anecic groups, trophically associated with the organic horizon and topically associated with the mineral soil horizons, in a variation of the C stocks in both the organic and mineral soil horizons, was also taken into account.

Texture of Soil-forming Rocks
Results of the one-way analysis of variance demonstrated that more than 80% of the variation in the contents of clay fraction (<2 µm) in the samples from soil-forming rocks is explained by the belonging of plots to regions ( Table 1). The lowest content of clay fraction was found in BP forests (1.6%), and the highest in NC forests (23.1%). The content of clay fraction ranged from 15.9% to 31.3% in all types of NC forests (Table S1), from 17.8% to 27.8% in MO forests (Table S2), and from 0.5 to 2.6% in BP forests (Table S3). Table 1. One-way ANOVA of the region and forest type factors, affecting soil forming rock, vegetation, litter quality, and soil macrosaprophages.

Vegetation
Significant differences were found between NC forests, on the one hand, and MO and BP forests, on the other hand, in terms of floristic composition (Figure 2), stem wood stock, the projective cover of vegetation layers, and SR (Table S1). More than 60% of the variation in the projective cover of the herb layer was attributed to a factor such as region, while for other variables, the proportion of the variance attributed to this factor was less than 31% (Table 1). The highest projective cover of the herb layer dominated by Carex pilosa Scop. and Oxalis acetosella L. was observed in MO forests. The lowest projective cover was found in NC forests, where the herb layer is poorly developed (5%-20%) because of shading.
The highest SR (27) was found in MO forests, whereas in NC and BP forests it was significantly lower: 18 and 22, respectively. The SR of woody plants in the herb layer of NC and BP forests was on average 9, which is higher than the total SR of woody plants in shrub and tree layers for these regions, which was 4 and 7 species on average, respectively. In MO forests, the SR of woody plants in the herb layer and in the shrub and tree layers together was on average 6 and 9, respectively (Table S5).
The largest stem wood stock was recorded in NC forests (125-724 t/ha), whereas in MO and BP forests, it ranged from 103 to 269 t/ha. Northwest Caucasus. NC1 and NC2 forests with herbs have a high level of similarity and differ significantly from the old-growth intact NC3 forests with sparse ground vegetation in their floristic composition ( Figure 2) and in ecological and phytocoenotic variables (Table S2). Stem wood stock in NC1 and NC2 forests ranged from 125 to 215 t/ha, and in NC3 forests reached 724 t/ha.
Average values of the total SR and the projective cover of the herb layer in the younger NC1 and NC2 forests were higher than in NC3 forests (Table S2) due to the differences in light availability under the canopy, which was lower in the old-growth forests. Lower canopy closure and height of the tree layer in NC1 and NC2 forests contribute to more favorable conditions for the development of ground vegetation. The highest SR was observed in NC1 and NC2 forests and varied from 13 to 33 species, the lowest being in NC3 forests, ranging from 3 to 14 species.
The SR of woody plants in the herb layer was significantly higher than in the shrub and tree layers together in all forest types. The most significant difference was observed in NC1 forests (11.5 and 4 species on average, respectively) and in NC2 forests (12 and 6 species on average, respectively). In NC3 forests, SR of woody plants in the herb layer was higher than in the shrub and tree layers together (on average 4.1 and 1.4 species, respectively), with almost all woody plants in the herb layer being trees (3.9). Because of shading in NC3 forests, the ground vegetation, and the undergrowth were very poorly developed, the average projective cover of the herb layer was 5%-10%. Of all the trees in the herb layer, only Abies nordmanniana and Fagus orientalis reached the shrub and tree layers, while the other tree species created a mass effect [55] and had no chances for development outside the herb layer.
Moskva-Oka plain. MO1 forests and MO2 forests differed significantly in their floristic composition ( Figure 2) and ecological and phytocenotic variables (Table S3). Stem wood stock in MO1 forests ranged from 121 to 166 t/ha, whereas in MO2 forests, it was from 171 to 184 t/ha. Compared to MO1 forests, the herb layer of MO2 forests had lower cover (89.7% and 73.1%, respectively), but higher species diversity (22.1 and 35.6, respectively) ( Table S3). Differences in a projective cover of the herb layer were due to light conditions. The higher species diversity of the herb cover of MO2 forests was due to the contribution of both boreal and nemoral species, whereas in MO1 forests, nemoral species absolutely dominated, and boreal species were practically absent. The MO forests were characterized by a low contribution of woody species in the herb layer.
The SR of woody plants in the herb layer was almost always lower than in the shrub and tree layers together. The clear differences between the number of tree species in tree and herb layers were observed in both MO1 forests (on average 4 and 8 species, respectively) and in MO2 forests (7 against 11 species, respectively).
Bryansk Polesie. The BP1 forests differed significantly from the BP2 forests and BP3 forests in the species composition of vascular plants ( Figure 2) and in ecological and phytocenotic variables (Table S4). Stem wood stock in young BP1 forests ranged from 103 to 116 t/ha, whereas in older BP2 forests and BP3 forests it was from 180 to 268 t/ha. The SR and projective cover of the herb layer in BP3 forests were higher than in BP2 forests and BP1 forests (Table S4). The highest SR was observed in BP3 forests and varied from 32 to 39 species, the lowest in BP1 forests, ranging from 10 to 19 species, and BP2 forests occupied an intermediate position with SR amounting to 20 to 26 species. This is due to the complex mosaic structure of BP3 forests with treefall gaps and the high contribution of broadleaf trees and shrubs.
The SR of woody plants in the herb layer was almost always higher than in the shrub and tree layers together. The most significant difference was observed in BP3 forests (11 and 8 species, respectively) and in BP2 forests (10 and 7 species, respectively). In BP1 forests, the difference between the SR of woody plants in the herb layer and in the shrub and tree layers together was also found (on average 7 and 5 species, respectively).

Litter Quality
The region factor contributed significantly to the variance of the litter quality variables, such as BS, N, and pH in the L-layer (Table 1). A comparison between the objects of the three regions (Table S1) showed that the lowest litter quality, with a relatively low BS (86%) and pH (5.5), was typical for BP forests, while NC forests were characterized by a high BS (93%), which is evidence of the high litter quality. The high pH values in the litter were found in the MO plots (5.9).
Forest types, except for MO forests, contributed significantly into variance of BS, N, pH, and the C/N ratio in the L-layer (Table 1).
Northwest Caucasus. When comparing litter quality, it was found (Table S2) that BS varied slightly, from 94% in NC1 forests to 92% in the NC3 forests. The lowest litter C/N ratio (23) was found in NC1 forests. The pH values and N contents in the L-layer of all forests were comparable (Table S2). The contribution of forest types to the variance of the litter quality was around 10%.
Moskva-Oka plain. The litter quality variables, such as BS, N, pH, and C/N, in MO1 forests and MO2 forests were comparable (Table S3).
Bryansk Polesie. When comparing the litter quality of different types of forest (Table S4), it was found that the BP3 forests' litter had the highest BS values (91%), pH values (6.1), and N content (2%), whereas BS, pH, and N contents in the L-layer of BP1 forest's litter were significantly lower (79%, 4.7, and 1.4%, respectively). The highest C/N ratio in the litter was found in BP1 forests (34) and the lowest in BP3 forests (21). The BP2 forests occupied an intermediate position, and in terms of litter quality, forests of this type were closer to broadleaf forests (Table S4). The contribution of forest types to the variance of the litter quality in BP forests was more than 50% (Table 1).

Soil Macrosaprophages
The composition and functional diversity of the most numerous macrosaprophage taxa between regions varied less significantly than their biomass. The Simpson's diversity index (D) was 1.5, 1.6, and 1.2 in NC, MO, and BP forests, respectively, but differences in the biomass of macrosaprophages between the forests of different regions ranged from 2.5 to 15 times (Table S1). In forest types from three regions, the earthworms of the Lumbricidae family contributed significantly to total biomass: 80% in NC forests, 83% in BP forests, and 97% in MO forests. Earthworms were not found in one type of BP forest only, namely, in BP1 forests (pine forests). Statistically significant differences between the regions in the biomass of epi-endogeic, endogeic, and anecic earthworms were found (Table 1). In NC forests, fauna of earthworms was characterized by a predominance of endemic species and by species with a limited range, whereas the fauna of MO and BP forests was mainly represented by cosmopolitan species. In NC forests, a high contribution to biomass was also made by millipedes of the Julidae family-up to 10%. In all three regions, small larvae of mosquitoes and flies (order Diptera), two-legged millipedes (class Diplopoda), woodlice (order Isopoda), and gastropods (class Mollusca) were found in the litter.
No differences in the composition between other groups of macrosaprophages were revealed. In three types of forest, other macrosaprophages inhabit the litter: suborder Oniscidea, fam. Julidae, order Polydesmida, fam. Limacidae, fam. Ectobiidae (imago), fam. Tipulidae (larvae), whose average biomass was 2.3 g/m 2 . A high contribution (3%-10%) to the total biomass of all macrosaprophages, other than earthworms, was made by the white vegetable-eating millipede Pachyiulus krivolutskyi (fam. Julidae)-the dominant representative of the order Julida in Russia, endemic in the Western Caucasus.
Moskva-Oka plain. In MO forests, the total biomass of macrosaprophages averaged 29.0 g/m 2 . The average biomass of earthworms was 28.4 g/m 2 . Eight species of earthworms belonging to four functional groups were identified; they were all cosmopolitan. As compared to the other regions, MO forests demonstrated significant differences in the biomass of epi-endogeic, endogeic, and anecic groups of earthworms (Table 1,  Table S3). The highest contribution to biomass was made by endogeic and epi-endogeic worms (Table S3).
The MO1 forests provided habitat for seven species of earthworms: epigeic D. octaedra, L. castaneus, epi-endogeic L. rubellus, endogeic A. caliginosa, A. rosea, anecic L. terrestris, and A. longa. Although four groups of Lumbricidae were represented, the biomass of endogeic species significantly exceeded the biomass of worms of other groups (Table S3). In addition to earthworms, other macrosaprophages live in the litter: larvae of Diptera, woodlice (order Isopoda), slugs of fam. Limacidae, and two-legged millipedes (fam. Julidae), and their total biomass was 0.3 g/m 2 .
Bryansk Polesie. In BP forests, the total macrosaprophage biomass averaged 2.2 g/m 2 . The fauna of earthworms was represented by cosmopolitan species, with the exception of the east Asian species E. nordenskioldi, for which Bryansk Oblast marks the western border of its habitat. Among the functional groups of worms, epigeic earthworms (Table S1) and epi-endogeic earthworms dominated in biomass. No earthworms were found in BP1 forests. Moreover, no earthworms were found in these forests even in their potentially favorable habitats, i.e., depressions, pits, and dead wood in late stages of decomposition. In the L-layer, other macrosaprophages were found: larvae of order Diptera, imago of fam. Ectobiidae, larvae, and imago fam. Geotrupidae, with a total biomass of 0.36 g/m 2 . Three species of Lumbricidae were found in the BP2 forests: epigeic D. octaedra, D. r. tenuis, and epi-endogeic E. nordenskioldi. Earthworm biomass was low (2.5 g/m 2 ) and dominated by epigeic earthworms (Table S4). Three species of Lumbricidae were found in BP3 forests: epigeic D. octaedra, D. r. tenuis, and epi-endogeic L. rubellus. Earthworm biomass was also low (2.2 g/m 2 ) and dominated by epi-endogeic earthworms (Table S4). The population of other macrosaprophages in BP2 and BP3 forests did not differ from that of BP1 forests neither in the composition of the main taxa, nor in the biomass (Table S4).

Soil Carbon Stock
The results of the one-way analysis of variance demonstrated a significant contribution of such factors as region and forest type to the variation of C stocks in the organic and mineral layers of the soil ( Table 1). The region factor provided 37% of the variation in C stocks in the FH-layer and more than 40% of the variation in C stocks of mineral horizons ( Table 2).
A comparison between the three regions showed that in the FH-layer, the lowest C stock was in MO, while BP forests were characterized by the highest C stock: 0.9 and 5.2 t ha −1 , respectively ( Figure 3; Table 2). However, NC forests were characterized by the highest C stock in the 30-cm and 30-50 cm mineral layers (83 and 26 t ha −1 , respectively) ( Figure 4; Table 2). The highest total soil C stock, calculated taking into account the FH-layer and the 50-cm mineral layer, was found in NC. The share of the FH-layer in the total C stock decreased from 19% in BP to 3% in NC and 2% in MO.  Northwest Caucasus. The highest C stock in the organic horizons was found in NC2 forests (2.5 t ha −1 ), and the lowest in NC1 forests (1.7 t ha −1 ) (Figure 3; Table S2). The NC3 forests accumulated the least amount of C in the 0-30 and 30-50 cm layers (58 and 16 t ha −1 , respectively), while NC2 forests were characterized by the highest C stock in these layers (98 and 35 t ha −1 , respectively) and the NC1 forests occupied an intermediate position (91 and 28 t ha −1 in the 0-30 and 30-50 cm layers, respectively) ( Figure 4; Table S2). Forest types provided 7 % of the variation in FH-layer C stocks and more than 40% of the variation in C stocks of mineral horizons ( Table 2).
Moskva-Oka plain. In MO, differences in the soil C stock between the forest types dominated by different tree species were also evident ( Table S3). As compared to MO2 forests, MO1 forests were characterized by lower C stock in the FH-layer (0.1 t ha −1 ) ( Figure 3) and mineral layers 0-30 and 30-50 cm (60 and 14 t ha −1 , respectively) ( Figure 4). Forest types provided 26% of the variation in C stocks of the FH-layer and about 30% of the variation in C stocks of mineral horizons (Table 2). Asterisks indicate the significance of the differences between the mean in category and overall mean: *-p ≤ 0.1, **-p ≤ 0.01, ***-p ≤ 0.001. Bryansk Polesie. When considering BP forests, it turned out that the highest C stock in the FH-layer was in BP1 forests where it was more than two times higher than in BP3 forests: 8.6 against 3.4 t ha −1 , respectively ( Figure 3; Table S4). The BP3 forests were characterized by the low C stocks in the mineral layer 0-30 (33 t ha −1 ) (Figure 4; Table S4). The BP2 forests were characterized by the highest C stock in the mineral layers 0-30 (51 t ha −1 ). The C stock in the 30-50 cm layer was comparable in all three forest types and ranged from 7 to 9 t ha −1 . Forest types in BP provided 48% of the variation in C stocks of the FH layer and more than 20% of the variation in C stocks of mineral horizons (Table 1).

Contribution of Vegetation and Earthworms to Soil Carbon Stock Variation
Of all the variables included in the hierarchical analysis, ones with the highest contribution to the variation of C stocks in the organic and the 0-50 cm mineral layer were selected: (1) stem wood stock; (2) the SR of woody plants in herb layer and total species richness in shrub and tree layers; (3) the C/N ratio and BS in L-layer; and (4) the biomass of functional groups of earthworms associated with the organic or mineral layers.
The hierarchical analysis showed that about 13% of the variation of C stocks in the FH-layer at all study sites was made by biomass of earthworms associated with the organic horizon, and about 18% and 23% of the variance was due to the C/N ratio and BS in the litter (the L-layer), respectively. The individual independent contribution of these variables was 21%, 30%, and 38%, respectively. The correlation of C stocks in the FH-layer with the earthworm biomass and BS was negative and with the C/N ratio was positive. The biomass of functional groups of earthworms associated with the organic horizon had an I/J ratio of less than 1 (Figure 5), which means that this indicator was closely related to other predictors (in particular, BS). The richer the litterfall and, correspondingly, the litter is in nutrients, and the higher the biomass of earthworms, the faster the organic horizon decomposes, which leads to a decrease in its stock and, accordingly, carbon stocks in it. The individual independent contribution of predictors that have an insignificant and independent effect on C stock, such as SR of woody plants in the herb layer, SR of plants in the shrub and tree layers, and stem wood stock, to the variation of carbon stocks was no more than 5% ( Figure 5). The biomass of earthworms trophically and topically associated with the mineral soil horizons (endogeic, epi-endogeic, and anecic groups) contributed significantly to the variation of C stocks in the 50 cm mineral layer (20% of variance), and the contribution of BS was 21%. The contribution of SR of woody plants in the herb layer was also high (16%). The individual independent contribution of these variables was 36%, 38%, and 9%, respectively. The SR of woody plants in the herb layer had the highest independence, as indicated by the highest I/J ratio (Figure 4). Unlike in the organic horizon, the correlation with the earthworm biomass and BS becomes positive. With the increase of BS in the litter layer, the SR of woody plants in the herb layer, and the earthworm biomass in the soils, the C stocks in the mineral horizons of the soils increased. The individual independent contribution of predictors that have an insignificant independent effect, such as total SR in shrub and tree layers and tree biomass, to the variation of carbon stocks was no more than 6%.

Climatic Factors and Soil Carbon Stock
The leading role of climate conditions in the soil C stocks was emphasized at the regional level [3]. Soil carbon stocks are positively correlated to average annual temperature, annual precipitation, and, therefore, to net primary productivity [2].
The influence of climatic factors on the soil carbon stock could be mediated by vegetation, its species composition and productivity. Warmer 5 • C and wetter 500 mm climatic conditions of NC forests compared to BP and MO forests explain the pronounced differences in the species composition of regional floras, as well as ecological and phytocenotic variables of forest communities. These differences between NC, BP, and MO forests are due to the different biology of the dominant trees and the duration of the growing season, which leads to differences in primary productivity. The more productive communities produce higher litterfall amounts. According to the existing data, the annual litterfall in MO forests ranged from 2.5 to 4.4 t ha −1 [56], in BP forests from 3.1 to 4.4 t ha −1 [57], and in NC forests from 3 to 12.2 t ha −1 [58]. Stem wood stocks increased from 156 t ha −1 in MO forests to 300 t ha −1 in the NC forests. The highest stocks of stem wood, as well as the highest soil carbon stocks, were characteristic of NC forests (Figure 3) growing in more favorable climatic conditions. Thus, the influence of climate on C stocks in the mineral horizons of soils was confirmed. Carbon stocks in the organic horizon did not follow similar patterns ( Figure 2) and were related to forest types.

Texture of Soil-Forming Rocks and Soil Carbon Stock
Compared to the NC and MO forests, developing on soil-forming rocks with high content of fine particles, the C stocks in the soil mineral layers in BP forests formed on sandy rocks was significantly lower. The impact of fine particles on carbon stocks could be explained by the formation of organo-mineral complexes with varying degrees of stability [59], which is associated with the qualitative composition of clay fractions and their absorption capacity [60,61]. However, it is worth noting that compared to NC and MO regions, the content of clay in the soil-forming rock of BP region was an order of magnitude lower, but the stocks of soil carbon in forests were lower only by a factor of two. This means that other factors, first of all, biotic ones, are of great importance in the soil C stock variation.

Vegetation and Soil Carbon Stock
Forest plants provide the soil with organic matter, primarily through above-and underground litterfall. To assess the role of vegetation in the accumulation of carbon in the soil, it is important to take into account not only the quantity but also the quality of the litter.
According to hierarchical analysis, the contribution of litter quality to the C stocks in the FH-layer significantly exceeded the contribution of stem wood stock (Figure 4), which indirectly reflects the amount of tree litterfall.
Plants belonging to different species can produce litterfall of different quality. Plant species can be combined according to the quality of their litterfall into functional types [62]. The highest C stocks in the 50-cm mineral layers were found in forests with higher diversity of plant species with mixed litterfall, such as NC2 forest type and in MO2 forest type ( Figure 3). The explanation for that could be the fast decomposition of the high-quality litterfall of nemoral herbs and deciduous trees causing the intense C fluxes into the mineral layers. However, at the same time, there was a relatively high C stock in the FH-layer in these types of forests, which could be explained by a contribution of fir trees to the shrub layer in NC2 forest type and of spruce trees to the tree layer in MO2 forest type ( Figure 2). It is well-known that litterfall of these conifers decomposes slowly which is an explanation for the accumulation of the organic horizon mass. In the BP region, as well as in MO and NC, the highest stocks in the 50-cm layer were also found in the forests with mixed litterfall (BP2 forest type), where, along with pine needles, the litterfall of deciduous trees and nemoral herbs contributed significantly to total litterfall amount.
Compared to NC2 forest type, the lower C stock in the soil mineral layers of NC3 forest type dominated by beech and fir could be explained by almost undeveloped herb layer (3% of cover) and single deciduous trees with high-quality, fast decomposing litterfall.
Along with the litterfall quality, an informative predictor of C stock in the FH-layer and the 50-cm mineral layer was the SR of woody plants in the herb layer. This parameter can be considered an integral indicator that takes into account the diversity of woody species of all vegetation layers, since the renewal of all species occurs directly in the herb layer. This indicator also demonstrates the full potential of woody plant diversity in the upper layers, that is, it reflects not only the modern composition of the upper layers, but also woody species potential and, therefore, the past history, when other species of woody plants could dominate in the upper layers and contribute to C stocks in the mineral layers of soil. A positive relationship between the soil C content and tree diversity in the herb layer was shown earlier [63]. According to hierarchical partitioning analysis, unlike in the herb layer, total SR in the tree and shrub layers makes a less significant contribution to the variation of soil carbon stocks.
In NC and BP forests, the SR of woody plants in the herb layer is almost always higher than in the shrub and tree layers together. This is due to favorable conditions for seed germination and the introduction of the entire range of the species which are typical of these areas, as well as due to a richer composition of shrubs. The mass effect phenomena [55] indicate a high heterogeneity of conditions in NC and BP forests.
In MO forests, the SR of woody plants in the herb layer is almost always less than the total SR in the shrub and tree layers. In these forests, Carex pilosa has a high projective cover, which could prevent the seeds of woody plants from reaching the soil and their germination. Moreover, the existing data [64] indicate that Carex pilosa suppresses competitors, mainly due to its powerful and well-developed root system.
Thus, the species diversity of plants influences the soil C stock: the highest C stock in the mineral layers of soil was found in mixed forests with a high diversity of woody and herb plants with different quality of the litterfall.

Litter Quality, Earthworms, and Soil Carbon Stock
Litter quality is related to the composition of vegetation and is an informative predictor of the soil C stock. According to hierarchical analysis, BS and the C/N ratio in the litter, being the indicators of its quality, are the informative predictors of the C stock in the FH-layer and mineral 50-cm layer.
The activity of soil saprophages processing litterfall is of great importance for soil carbon accumulation ( Figure 4). Earthworms, millipedes (Diplopoda), isopods (Oniscidea), molluscs, and larvae of Diptera [9,44,[65][66][67] are the main macrosaprophages that ensure the decomposition of litterfall in coniferous-broadleaf forests of the European part of Russia and regulate soil C stock. Earthworms usually significantly predominate in biomass among other groups of macrosaprophages [66,[68][69][70]. The informative indicator of the activity of soil saprophages is the biomass of different functional groups of worms.
High C stocks in FH-layer of BP forests could be explained by the low biomass of earthworms whose activity is limited in sandy soils [71].
It is noteworthy that the biomass of the functional groups of worms associated with the organic horizon (epigeic, epi-endogeic, and anecic) is closely related to litter quality. Low litter quality with a low BS, high acidity, and C/N ratio is a limiting factor for earthworms in BP1 forests. Pine needle litterfall is N-poor (often below 0.4%) [13], boreal dwarf shrub litterfall is rich in polyphenolic compounds [72], and green mosses have a low content of nutrients [73]. In MO and NC forests, as well as in BP3 forest type, plants with high-quality litterfall prevail. In these forests, litterfall is relatively rich in N: birch 0.7%, beech 0.9%, aspen 1.0%, hornbeam 1.1%, oak 1.2%, maple 1.3%, and lime 1.5% [74]. The richer the litterfall is in nutrients, the faster the litter is decomposed by soil biota, which leads to a decrease in the organic horizon mass and, accordingly, the C stocks in it.
The functional composition of earthworms strongly depends on the composition of the litterfall. Mixed coniferous-deciduous litterfall formed by a high variety of plant species is favorable both in trophic and topical terms: the hard-to-decompose fractions advantage the increasing organic horizon mass, which serves as a habitat for earthworms, and the easy-to-decompose fractions provide a trophic resource of high quality [20,75]. In such forests, as a rule, there is a rich functional composition of earthworms associated with the organic horizons and mineral layers of the soil. In mountain NC2 forest type and in plain MO2 forest type, the high C stocks in the FH horizon and in the 50-cm mineral layer could be explained by the effects of mixed litterfall with both fast and slowly decomposed components (Figure 2). A poor quality litterfall of spruce and fir trees explains the accumulation of the organic horizon mass, which is a habitat for saprophages that are functionally related to this horizon [11]. At the same time, the significant contribution of herbs, along with the high-quality litterfall of deciduous trees, ensures the active processing of the litterfall by soil saprophages. These forests also demonstrated the highest diversity of functional groups of earthworms. The functional groups of earthworms of the organic horizon feeding the organic matter, thus reducing the C stock in the FH layer. While epi-endogeic worms contribute to humus formation, endogeic earthworms provide organic matter migration into mineral horizons and active bioturbation of mineral soil layer [11,76,77]. In MO1 forest type and NC3 forest type with a low functional diversity of earthworms and low plant species richness, the C stocks in the mineral layers were low (Figure 3). The litter of homogeneous character does not advantage the diversity of earthworms in these forests [23].
In sandy soils, there are less favorable conditions for the activity of soil saprophages as compared to loamy soils [78]. In BP forests, the sandy soil is not conducive to the spread of worms, and saprophages are also represented by arthropods, which together with worms process the litterfall here. The highest C stocks in FH-layer are typical of BP1 forests, which is associated with low-quality litterfall and, consequently, low biomass of macrosaprophages. In BP2 and BP3 types, the litterfall quality is higher due to the contribution of nemoral herbs and deciduous trees, which results in increasing the activity of macrosaprophages causing a decrease in the FH mass and C stock in it.
The C stocks in the mineral horizons of soil of BP forests are mainly regulated by microbiota processing the litterfall due to the lack of endogeic earthworms, but abiotic factors, such as water regime, could also be important for soil C stock in these forests. In BP3 type of forests, the relatively low carbon stocks in the 0-50 cm layer may be associated with the stagnation of groundwater and the processes of gleying, which could contribute to the soil profile washability and C leaching outside the forest ecosystem. Widespread changes in the hydrological regime are currently seen in large areas of the European part of Russia [79].
Thus, there are evident links between the litter quality, earthworms, and soil C stocks in coniferous-broadleaf forests, but abiotic factors, such as particle size distribution in the soil-forming rocks and water regime, also influence soil C stocks.

Conclusions
This study confirmed the influence of climate and soil-forming rocks on the soil C stock and demonstrated the links between vegetation, earthworms, and soil C stocks.
The soil C stocks are regulated by the litter quality and activity of earthworms. The important predictors of the soil C stock variation are the litter quality indicators, such as base saturation and the C/N ratio. The highest C stocks in the organic horizons were associated with the low-quality litter, i.e., with a low base saturation, high acidity, and high C/N ratio. The highest soil C stocks in the mineral layers were found in the mixed forests, such as pine-broadleaf, spruce-broadleaf, and hornbeam-fir-beech ones, with the highest richness of plant species producing litterfall of different quality. Low-quality slowly decomposed litterfall of conifers and beech could contribute to the accumulation of the organic horizon mass, while fast decomposition of high-quality litterfall of deciduous tree species, such as lime, birch, and hornbeam, as well as of herbaceous plants, by soil biota could advantage the active down profile migration of carbon compounds into the mineral layers.
The carbon stock in the organic horizon was closely and negatively related to the biomass of the functional groups of worms that process the litter, while the carbon stock in the mineral layers was closely and positively related to the biomass of functional groups of worms whose life activity is related to the mineral part of the soil.
The content of clay fractions in soil, which is the stabilizing agent of organic matter, as well as the litter quality, regulates the activity of earthworms. The biomass and functional diversity of earthworms in forests on sandy soils were significantly lower compared to those on loamy soils. The highest biomass and functional diversity of earthworms in the mixed forests with the highest richness of plant species, producing litterfall of different quality, was explained by favorable topical and trophic conditions for worms.
Mixed litterfall, consisting of both fast and slowly decomposed components, produced by plants belonging to different woody and herb species, is a prerequisite for high activity of earthworms of different functional groups, and for a high soil C stock in coniferous-broadleaf forests. Thus, higher plant and earthworm functional diversity contribute to higher soil carbon stocks in coniferous-broadleaf forests.