Patterns of relative magnitudes of soil energy channels and their relationships with environmental factors in different ecosystems in Romania

The percentage compositions of soil herbivorous, bacterivorous and fungivorous nematodes in forests, grasslands and scrubs in Romania was analysed. Percentages of nematode abundance, biomass and metabolic footprint methods were used to evaluate the patterns and relative size of herbivory, bacterial- and fungal-mediated channels in organic and mineral soil horizons. Patterns and magnitudes of herbivore, bacterivore and fungivore energy pathways differed for a given ecosystem type and soil depth according to the method used. The relevance of herbivore energy channel increased with soil depth due to higher contribution of root-feeders. Ectoparasites, sedentary parasites and epidermal cell and root hair feeders were the most important contributors to the total biomass and metabolic footprints of herbivores. Metabolic footprint method revealed the general dominance of bacterial-based energy channel in all five types of ecosystems. The influence of altitude and climatic factors on percentages of abundance, biomass and metabolic footprints of herbivores, bacterivores and fungivores decreased with soil depth, whereas the influence of humus content, cation-exchange capacity and base saturation increased. Vegetation, altitude, climate and soil physico-chemical characteristics are important factors that influenced the abundance, biomass and metabolic footprints of herbivores, bacterivores and fungivores.


Results
The percentages of abundances of bacterivores, fungivores and herbivores in organic and mineral soil horizons. The order of the percentage of nematode abundances was herbivore > bacterivore > fungivore in organic horizon of coniferous forests and scrubs and bacterivore > herbivore > fungivore in deciduous and mixed forests and grasslands (Fig. 1A). The contribution of bacterivores and herbivores was significantly higher than that of fungivores in organic horizon in both grasslands and scrubs. In addition, herbivores were significantly more abundant in coniferous forests than in deciduous forests, mixed forests and grasslands. The percentage of fungivores was significantly higher in deciduous forests than in grasslands and scrubs. The percentage of bacterivores was significantly higher in deciduous and mixed forests, grasslands and scrubs than in coniferous forests.
The order of the percentage of nematode abundances was herbivore > bacterivore > fungivore in mineral horizon in all five ecosystems ( Fig. 1B and C). Herbivores had significant contribution to nematode communities in coniferous forests and grasslands, whereas contribution of fungivores was significantly higher in deciduous and mixed forests as compared to grasslands (Fig. 1B).
The percentages of biomasses of bacterivores, fungivores and herbivores in organic and mineral soil horizons. The order of the percentage of nematode biomasses was bacterivore > fungivore > herbivore in organic horizon in coniferous forests, mixed forests, grasslands and scrubs and fungivore > bacterivore > herbivore in deciduous forests (Fig. 1D). The biomass of bacterivores and fungivores exceeded that of herbivores in coniferous forests. The percentages of the three trophic groups were significantly different from each other in mixed forests and grasslands.
The order of the percentage of nematode biomasses was herbivore > bacterivore > fungivore in the 0-5 cm layer of mineral horizon (Fig. 1E) in coniferous forests, and bacterivore > herbivore > fungivore in deciduous and mixed forests, grasslands and scrubs. The biomass of fungivores was higher in deciduous forests and scrubs, whereas that of bacterivores was greater in deciduous forests, mixed forests and grasslands.
The order of the percentages of nematode biomasses was herbivore > bacterivore > fungivore in the 5-10 cm layer of mineral horizon (Fig. 1F) in coniferous forests and grasslands and bacterivore > herbivore > fungivore in deciduous and mixed forests and scrubs. The biomass of fungivores was higher in coniferous and mixed forests, whereas that of bacterivores was greater in deciduous and mixed forests and scrubs.
The percentages of metabolic footprints of bacterivores, fungivores and herbivores in organic and mineral soil horizons. The order of the percentage of nematode metabolic footprints was bacterivore > fungivore > herbivore in organic horizon of coniferous forests, mixed forests, grasslands and Scientific RepoRts | 5:17606 | DOI: 10.1038/srep17606 scrubs and fungivore > bacterivore > herbivore in deciduous forests (Fig. 1G.). Bacterivorous metabolic footprints were greater than herbivorous and fungivorous metabolic footprints in coniferous forests and scrubs. Bacterivorous and fungivorous metabolic footprints were significantly greater than that of herbivorous in deciduous forests.
The order of the percentage of nematode metabolic footprints was herbivore > bacterivore > fungivore in the 0-5 cm layer of mineral horizon in coniferous forests, and bacterivore > herbivore > fungivore in deciduous and mixed forests, grasslands and scrubs (Fig. 1H). Bacterivorous metabolic footprint was significantly greater than herbivorous and fungivorous metabolic footprints in deciduous forests and scrubs. Bacterivorous and herbivorous metabolic footprints were greater than fungivorous metabolic footprint in grasslands.
The order of the percentage of nematode metabolic footprints was herbivore > bacterivore > fungivore in the 5-10 cm layer of mineral horizon in coniferous forests and grasslands, and bacterivore > herbivore > fungivore in deciduous and mixed forests and scrubs (Fig. 1I). Herbivorous metabolic footprint was significantly higher than that of bacterivorous and fungivorous in coniferous forests. Bacterivorous and herbivorous metabolic footprints were greater than that of fungivorous in grasslands and scrubs. , biomass (D-F), and metabolic footprint (G-I) of the bacterivorous, fungivorous and herbivorous nematodes in organic horizon and 0-5 cm and 5-10 cm layers of mineral horizon of coniferous forest, deciduous forest, mixed forest, grassland, and scrub ecosystems in Romania. Bars represent means ± standard error. In each ecosystem, different uppercase letters indicate significant (P < 0.05) differences among the three nematode trophic groups. For each nematode trophic group, different lowercase letters indicate significant (P < 0.05) differences among ecosystem types. The total number of samples in each soil horizon and ecosystem type are given in Table 2 The percentages of abundance, biomass and metabolic footprints of different groups of plant feeders in organic and mineral soil horizons. The percentages of abundance, biomass and metabolic footprints of plant feeding nematode groups in organic horizon, 0-5 cm and 5-10 cm layers in the five ecosystems were shown in Fig. 2.
The percentage of abundance of plant feeders ( Fig. 2A-C) showed that epidermal cell and root feeders was by far the dominant group in organic soil horizon, being particularly favoured in forest and scrub soils, as compared to grassland soils. Semiendoparasites were more abundant in coniferous forests, scrubs and grasslands. Higher proportions of ectoparasites were found both in grassland and scrub soils. The abundance of these three groups increased, in general, with soil depth tested.
The percentages of biomass of plant feeders ( Fig. 2D-F) showed higher proportions of ectoparasites in coniferous forests, scrubs and grasslands, followed by sedentary endoparasites, more numerous in coniferous forests, mixed forests and scrubs. The abundance of epidermal cell and root feeders increased, especially in mixed forests and grasslands. Scrub soils were characterized by the highest abundance of ectoparasites in the 5-10 cm layer of mineral soil horizon. Again, a general trend in increasing the percentages of biomass of ectoparasites, epidermal cell and root feeders and sedentary endoparasites with soil depth could be noted.
The percentages of metabolic footprints of plant feeders ( Fig. 2G-I) revealed an almost similar pattern with the one showed for the percentages of biomass.
Relationships between environmental variables and the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in organic and mineral soil horizons. Canonical correspondence analysis showed that the first two axes of the environmental variables explained 11.7%, 11.9%, and 15.6% of the variances of the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in organic horizon and in the 0-5 cm and 5-10 cm layers of mineral horizon, respectively (Fig. 3). Monte Carlo permutation test results showed that different environmental variables were positively or negatively correlated with the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in organic and mineral soil horizons.
Altitude, coniferous forests, annual average precipitation and mixed forests significantly correlated with the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in organic horizon (P < 0.05) ( Table 1 and Fig. 3A). Particularly, the altitude and annual average precipitation correlated negatively with fungivores (Fig. 3A). Mixed forests correlated positively with bacterivores and coniferous forests correlated negatively with the bacterivores (Fig. 3A). Bi-plots of canonical correspondence analysis (CCA) of the relationship between environmental variables (ecosystem types, soil physico-chemical properties, climatic variables, and altitude) and the percentages of abundance, biomass, and metabolic footprint of the bacterivorous, fungivorous and herbivorous nematodes in organic horizon (A) and 0-5 cm (B) and 5-10 cm (C) layers of mineral horizon. Ordination diagrams presenting species scores (triangles) and environmental factor scores (vectors). Solid vectors indicate ecosystem types, empty vectors indicate soil physico-chemical properties, climatic variables and altitude. Ba, bacterivore; Fu, fungivore; He, herbivore; Prec, annual average precipitation; Tem, annual average temperature; TN, total nitrogen; K, available potassium content; P, available phosphorus content; HC, humus content; CEC, cation-exchange capacity; SH, total hydrolytic acidity; SB, total exchangeable base; and V, base saturation. Environmental variables of organic horizon do not include soil physico-chemical properties. Monte Carlo permutation test results for the environmental variables were shown in Table 1.
Coniferous forests and grasslands significantly (P < 0.05) correlated with the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in the 0-5 cm layer of mineral horizon (Table 1 and Fig. 3B). Particularly, coniferous forests correlated positively with herbivores and negatively with the bacterivores (Fig. 3B). Grasslands correlated negatively with fungivores ( Fig. 3B).
Coniferous forests, altitude, humus content, available potassium content, soil pH, cation-exchange capacity, and mixed forests significantly (P < 0.05) correlated with the percentages of abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in the 5-10 cm layer of mineral horizon (Table 1 and Fig. 3C). Particularly, coniferous forests correlated negatively with bacterivores and soil pH correlated positively with bacterivores (Fig. 3C). Cation-exchange capacity correlated positively with bacterivores and negatively with herbivores (Fig. 3C). Mixed forests correlated negatively with herbivores (Fig. 3C). Humus content and available potassium content correlated positively with fungivores whereas altitude correlated negatively (Fig. 3C).

Discussion
Nematode abundance, biomass and metabolic footprint methods generated different results of relative size of herbivore, fungivore and bacterivore energy channels for a given ecosystem type and soil depth.
This outcome is different from the one reported by Zhao and Neher 25 ,who described similar patterns of soil energy channel in grasslands, croplands and forests when nematode abundance and biomass methods were used.
There was a clear distinction between the relative magnitudes of herbivore, bacterivore and fungivore energy channels when comparing the three methods. When nematode abundance method was used, the relative magnitudes of herbivore energy channel appeared much greater, but much lower for bacterivore and fungivore energy channels as compared to the biomass and metabolic footprint methods. As a result, nematode abundance method overestimated the herbivory mediated channel and underestimated the bacterial-and fungal-mediated pathways. Although epidermal cell and root hair feeders were the most abundant plant feeders ( Fig. 2A,B), their biomasses ( Fig. 2D-F) and metabolic footprints (Fig. 2G-I) were very low.  Fig. 3; -indicates that further addition of the variable via forward selection did not improve model fitting and the Monte Carlo permutation test did not perform for the respective variable; /indicates that soil physico-chemical properties were only assessed for mineral horizon.
Interestingly, biomass and metabolic footprint methods generated similar patterns among the five ecosystem types for a given energy channel and soil depth. This result is attributed to the fact that biomass is closely related to the carbon utilization and the magnitudes of soil energy channels. Therefore, biomass and metabolic footprint methods may be more appropriate for evaluating soil energy channels as compared to nematode abundance method because they are related to matter and energy flows in the ecosystems. On the contrary, nematode abundance method does not take into account body size and carbon utilization by nematodes with different functional characteristics 31 and using this method in evaluating soil energy channels may be biased due to large variation in nematode body sizes or biomasses.
The factors that influence the soil energy channels are complex (Supplementary information). Vegetation identity had been shown to have strong impacts on soil nematode fauna 23,31,32 and was thought to have major role in controlling energy channels in soil systems 24 . Consistently, ecosystem type-characterized by vegetation composition-affected the energy channels in the present study (Supplementary information). Particularly ecosystem type could directly and/or indirectly affect the energy channels via differentiation in soil resource (Supplementary information, Fig. S1B). In forests, basal resources such as leaf litter and root-derived products are considered to be important in structuring soil food webs 24 .
Different quality and quantity of basal resources related to tree composition and understory vegetation in forest ecosystems (coniferous, deciduous and mixed) resulted in different patterns and magnitudes of herbivore, fungal and bacterivore energy channels, in line with the first hypothesis. However, remarkably, there are no striking differences regarding the patterns and relative size of the three energy channels between forest types and between forests, grasslands and scrubs. This result suggests that grouping sites according to plant associations, which are strongly related to soil characteristics and parent rock, may be more appropriate when evaluating soil energy channels of different ecosystems.
The importance of the herbivore energy channel increased in order organic layer < mineral 0-5 cm layer < mineral 5-10 cm layer, the shift being more evident when comparing organic soil horizon to mineral 0-5 cm layer and less obvious among the two layers (0-5 cm and 5-10 cm) of mineral soil horizon. The increasing importance of herbivore energy channel in mineral soil horizon was in relation to the higher contribution of root-feeding nematodes, more abundant in deeper soil layers. The importance and relevance of herbivore energy channel in grassland and forest ecosystems has previously been addressed by Zhao and Neher 25 , suggesting that it is responsible for a considerable part of the below-ground energy flow. The results shown here are therefore in line with the previous findings and stress the increasing importance of herbivorous-mediated channel in relation to soil depth.
Distribution of different groups of plant feeders associated with soil horizons is thought to reflect the distribution of plant roots at various soil depths 33,34 . Vertical migration of certain plant feeders to more stable and uniform conditions in deeper soil layers is believed to be a strategy to avoid fluctuations in soil moisture and temperature in upper soil horizons despite a reduced availability of feeding sites 35 .
Evaluated by the metabolic footprint method, the five ecosystems were generally dominated by the bacterial-based energy channel, suggesting bottom-up regulation (Wang, pers. comm.), except the coniferous forests in the 5-10 cm layer of mineral horizon, where herbivorous-based energy channel prevailed. In addition, the relative sizes of the fungal-based energy channel in organic horizon of deciduous forests, the herbivorous-based energy channels in the 0-5 cm layer of mineral horizon of coniferous forests and grasslands, and the herbivorous-based energy channels in the 5-10 cm layer of mineral horizon of grasslands and scrubs were similar to the relative sizes of their corresponding bacterial-based energy channels of an ecosystem. In general, the size of the fungal-based energy channel in organic horizon was relatively greater than that of the herbivorous-based energy channel, except the coniferous forests and scrubs. This pattern was reversed in the 0-5 cm and 5-10 cm layers of mineral horizon, where herbivorous-based energy channel increased in importance. The general dominance of bacterial-based energy channel in all five types of ecosystems surveyed is congruent with previous findings 25,27,[36][37][38][39][40][41] , indicating resource-rich substrates.
Unravelling patterns of interactions between soil energy channels and environmental factors is both challenging and captivating because offer insights on how soil food webs are driven and function. As noted above, various factors influenced soil energy channels, although resource quantity and quality may be the main controlling factors for a given ecosystem type 23,24,31,32 at small scale. However, at regional scale different environmental factors interacted, resulting more complex effects on energy pathways (Supplementary information). Specifically, not only soil resource but also climate, ecosystem type and soil environment showed significant influences on soil energy channels in forests, scrubs and grasslands. In organic horizon, the proportion of fungivores declined with altitude and precipitation, in line with previous findings 42 , probably as result of resource limitation due to climatic influence. The negative influence of altitude on fungivores maintained as well in the 5-10 cm layer of mineral soil horizon. Bacterivores were favoured in mixed forests, but responded negatively to acidic substrate of coniferous forests, more favourable to plant-feeding nematodes in the 0-5 cm layer of mineral soil horizon. The negative correlation pattern of bacterivores with coniferous forests was found in both organic and mineral soil horizons, suggesting that soil pH resulting from high phenol content of coniferous litter [43][44][45] is a strong regulator of bacterial feeding nematode populations in this type of ecosystem. Soil pH has previously been found to be an important limiting factor for survival of bacterial feeding nematodes 46 . Although higher amounts of organic matter were reported in coniferous stands as compared to deciduous and mixed ones, the Scientific RepoRts | 5:17606 | DOI: 10.1038/srep17606 concentration of microorganisms was much higher in deciduous litter 47 . This point towards a direct link between the densities of bacteria and bacterivores in soils of coniferous forests.
As expected, the influence of altitude and climatic factors on percentages of abundance, biomass and metabolic footprints of herbivores, fungivores and bacterivores decreased with soil depth, whereas the influence of certain soil physico-chemical parameters (e.g. humus content, cation-exchange capacity, base saturation) increased.
The results bring new information on how different biotic and abiotic factors correlated with percentages of abundance, biomass and metabolic footprints of herbivores, fungivores and bacterivores at various soil depths. They reveal that vegetation, altitude, climate and soil physico-chemical characteristics are important factors that influence the abundance, biomass and metabolic footprints of herbivores, fungivores and bacterivores in forests, scrubs and grasslands, in line with the second hypothesis.

Concluding remarks
The analysis of trophic diversity of nematode communities in five ecosystem types (coniferous forests, deciduous forests, mixed forests, grasslands and scrubs) in Romania revealed different patterns and magnitudes of bacterivore, fungivore and herbivore energy channels in organic and mineral soil horizons.
Different relative sizes of herbivore, fungivore and bacterivore energy pathways for a given ecosystem type and soil depth were assessed when percentages of nematode abundance, biomass and metabolic footprints methods were used. Biomass and metabolic footprint methods generated similar patterns among the five ecosystem types for a given energy channel and soil depth. These two methods are considered to be more appropriate for assessing soil energy channels as they are more connected to C and energy flows in soil systems.
The role and relevance of herbivore energy channel increased in all ecosystem types with soil depth due to higher contribution of root-feeders in mineral soil horizon as compared to organic soil horizon. Ectoparasites, sedentary parasites and epidermal cell and root hair feeders were the most important contributors to the total biomass and metabolic footprints of herbivores.
Metabolic footprint and biomass methods were more appropriate in evaluating the relative magnitudes of soil energy channels and both of them revealed the general dominance of bacterial-based energy channel in all five types of ecosystems.
The influence of altitude and climatic factors on percentages of abundance, biomass and metabolic footprints of herbivores, fungivores and bacterivores decreased with soil depth, whereas the influence of humus content, cation-exchange capacity and base saturation increased.
The results showed that vegetation, altitude, climate and soil physico-chemical characteristics are important factors that influence the abundance, biomass and metabolic footprints of bacterivores, fungivores and herbivores in forests, scrubs and grasslands.
Resource quantity and quality may be the main driving factors of soil energy channels in ecosystems at small scale, but multiple interactions between various environmental factors at regional scale result in more complex effects on energy pathways.

Materials and Methods
Site description. A total of 135 sites distributed in the Romanian Carpathians (128 sites) and the Transylvanian Plateau (7 sites) were investigated. They included representative coniferous (n = 27), deciduous (n = 43) and mixed (n = 26) forests, grasslands (n = 31), scrubs (n = 7) and also a European black pine (Pinus nigra Arnold) plantation. Coniferous forests consisted mostly of Norway spruce (Picea abies L.). Most of the deciduous forests included beech (Fagus sylvatica L.), alone or in association with hornbeam (Carpinus betulus L.); other species were sessile oak (Quercus petraea (Matt.) Leibl) and sweet chestnut (Castanea sativa Mill.). Mixed forests were mostly composed of Norvegian spruce (Picea abies L.), beech (Fagus sylvatica L.), and occasionally silver fir (Abies alba Mill.). Mountain dwarf pine (Pinus mugo Turra.) was present in all scrubs. The sites were situated along an altitudinal gradient from 180 to 2,350 m a.s.l., those located higher in the Carpathians were generally less affected by human activities, some of them being located in pristine natural protected areas. Sampling was carried out between 1974 and 1999, during an extensive survey focused on Romanian biodiversity. Sampling, extraction, identification. Three to seven 25 m 2 plots were selected at random from a surface of about 500 m 2 , which was considered to be representative for the ecosystem type. Roughly, the distances between plots were at least 10 m and therefore sample replicates were in fact pseudoreplicates. In each plot, one sample comprising 10 cores (cores of 25 cm 2 in surface) was collected from the organic horizons (O and AO respectively, see explanation below) and A layer (but sometimes also from B layer, if present) of the mineral soil horizon (0-5 cm and 5-10 cm; cores of 3.8 cm 2 in surface), separately. This resulted in a total of 597 samples including organic horizons only, 378 samples containing exclusively the top 5 cm of the A layer of mineral horizon and 291 samples including solely the 5-10 cm of the A (and sometimes B) layer of mineral horizon, respectively (Table 2). According to the Romanian System of Soil Classification 48 , the organic horizon (O) included litter, fermentation and humification layers (which characterized most forests), while the sod (or turf) horizon (AO) was the shallow surface soil horizon with high proportion of matted roots in grassland soils. Nematodes in these samples were extracted using the centrifugation method 49 , after which they were counted and fixed with 4% formaldehyde solution heated at 65°C. Temporary mass-slides were prepared from each sample and used for examination. At least 150 individuals were randomly identified at the microscope at 400 × magnification and were used to estimate the relative abundance. Adults were identified to genus level whereas juveniles were generally assigned to families. The classification system used was according to Andrássy 50 , Siddiqi 51 and Bongers 52 . The contribution of trophic groups was established according to Yeates et al. 19 A separate set of samples from the A layer of mineral horizon was taken for assessing soil physico-chemical properties, as this horizon was considered less exposed to climatic factors and with more stable and uniform conditions. Standard methods were used for the determination of soil pH (in water), available phosphorus (P 2 O 5 ) (ppm) and potassium (K 2 O; ppm), humus content (%), total nitrogen content (%), total exchangeable bases (SB; (meq 100 g −1 soil), total hydrolytic acidity (SH; meq 100 g −1 soil), cation-exchange capacity (CEC; %), and base saturation (V; %) 53,54 .
Data on annual average temperature and precipitation was compiled from the CARPATCLIM-Climate of the Carpathian Region online database (CARPATCLIM Database, 2013; www.carpatclim-eu.org) for each sampling site, based on its geographical coordinates and using a running mean (e.g. sampling year − 11 up to sampling year − 1), as sampling was done in several field campaigns from 1974-1999. Data analysis. Three trophic groups of nematodes (plant feeding, bacterial feeding and fungal feeding) were taken into account to evaluate soil energy channels 25 . Omnivore and predatory nematodes were excluded because they are located in higher trophic levels in soil food webs than bacterivores, fungivores and herbivores and are not the most relevant or direct groups to assess the relative magnitudes of the herbivory, bacterial-and fungal-mediated channels 5,25 . Percentages of abundance, biomass and metabolic footprint were calculated for each of these trophic groups in each sample. Mean nematode biomass of a genus was assigned as mean biomass of the respective taxonomic family 55 . When mean biomass of the respective taxonomic family for a given genus was not provided in Ferris 55 , biomass of the genus given in Zhao et al. 56,57 was used. Calculation of metabolic footprint for each trophic group was as follows: where Nt is abundance of the t taxa in a sample, Wt is the biomass of the t taxa, and mt is the cp-value of the t taxa. Percentages of the abundance, biomass and metabolic footprint of plant feeding, bacterial feeding and fungal feeding nematodes in the five ecosystem types, at different soil depths were analyzed separately by one-way ANOVA. Briefly, the differences of a given energy channel among ecosystems at each soil depth (27 variables, 3 trophic groups × 3 soil depths × 3 methods) and the differences among the three energy channels for a given ecosystem at a given soil depth (45 variables, 5 ecosystems × 3 soil depths × 3 methods) were analyzed. The percentage data were firstly arcsine transformed to meet assumptions of normality and homogeneity of variance. Alternatively, the data were natural log, square root, or rank transformed to meet the assumptions above. Canonical correspondence analysis (CCA) 58 was performed using CANOCO 4.5 software (Ithaca, NY, USA) to determine the relationships between soil energy pathways and environmental variables, and between soil energy pathways and ecosystem types, with environmental variables as co-variables. Statistical significance was determined at P < 0.05. ANOVAs were performed using SPSS software version 16.0 (SPSS Inc., Chicago, IL, USA). LSD was used to test differences among treatment means; and Tamhane's T2 was used to test differences among treatments when variances of transformed data were unequal. Forward selection and test of significance for the first axis were based on Monte Carlo permutation test (n = 499).