Pedodiversity and Organic Matter Stock of Soils Developed on Sandstone Formations in the Northern Apennines (Italy)

: Pedodiversity is considered the cornerstone of biodiversity. This work aimed to (1) assess pedodiversity according to vegetation, topographic factors, and lithology and to (2) identify the major soil-forming factors on soil organic matter (SOM) stock at a 0–30 cm depth. These goals were reached using data from 147 georeferenced soil profiles distributed along 400–1000 m ( ≤ 1000) and 1000–2134 m (>1000) altitudinal gradients in the northern part of the Apennine chain in Italy. Soils showed mainly weak or incipient development (i.e., Entisols and Inceptisols), which could be attributed to sand-based lithology, high slope gradients, and low SOM accumulation rates, which promote soil erosion processes. However, higher pedodiversity was observed at >1000 m than at ≤ 1000 m, likely due to the higher vegetation cover diversity and climate variability; Spodosols and Mollisols were also found. A greater SOM stock was found at >1000 than ≤ 1000 m, and vegetation seemed to not affect SOM amounts, suggesting a greater influence of climate on SOM content compared to vegetation. Considering ecosystem conservation, the observed spatial pedodiversity could be considered a critical basis for the protection of soil resources and pedodiversity itself in mountain regions.


Introduction
Biodiversity evaluation, also referred to as spatial diversity, is a major topic in ecosystem studies and is mostly applied to aboveground flora and fauna [1,2] but also recently to belowground biodiversity [3,4]. Conversely, few investigations concerning pedodiversity were carried out [5][6][7], although pedodiversity is considered the cornerstone of biodiversity given the strong relationship between soil and flora and fauna [8,9]. Pedodiversity is extremely complex and driven by factors involved in several soil formation processes [10,11], making it difficult to quantify the influence of each factor on soil properties. Nowadays, because of the major role of soil in storing organic matter [12], several studies focused on identifying the factors that mostly significantly affect soil organic carbon (SOC) [13][14][15]. Indeed, the spatial changes in SOC amounts are dependent on biotic and abiotic factors such as climate, the amount and quality of plant debris, vegetation cover, and land use [16][17][18][19]. In addition, SOC stock is affected by soil development according to the different weight of each pedogenesis factor [20,21]. In this regard, recent papers highlighted the importance of soil type knowledge in terms of lithology, the physicochemical characterization of organic debris, and pedogenetic horizons to determine soil forming processes, in order to understand the distribution and quality of SOC stock [17,22].
In recent decades, there has been growing interest in the knowledge of mountain soils because montane ecosystems are characterized both to have high pedodiversity [23][24][25] and

Parent Material
The sites under investigation were characterized by diverse sandstone types as parent rock materials ( Figure 1) whose chemical features [68][69][70][71] are reported in

Climate
The climate of the Northern Apennines is characterized by cold winters with frequent snowfalls and hot summers. The mean annual temperature of the study areas ranges from 14.6 • C at the lowest altitude to 3.3 • C at the highest (Table S2 of  (−3.3 • C), and the warmest ones were recorded in July and August. The mean annual cumulative precipitation was 1004 mm and ranged between 753 and 1703 mm (Table S1 of the Supplementary Materials) with an uneven distribution between the months. In addition, above 2000 m a.s.l., the snow covers the soils for about 6 months per year (from November to April), which also happens at Mount Cimone [64]. The mean air temperature and mean annual cumulative precipitation data were used to determine soil moisture and temperature regimes according to the Newhall Simulation Method [72]. With the approximations allowed by the Soil Taxonomy [52], it was possible to identify a mesic temperature regime between 200 and 1300 m a.s.l. and a frigid temperature regime between 1300 to 2000 m. The soil moisture regime was found to be udic between 200 and 1800 m a.s.l. and perudic over 1800 m of altitude (Table 1). Table 1. Altitudinal range, soil temperature regime, and soil moisture regime of the phytoclimatic zones and subzones of the northern Apennine chain, Italy, according to De Philippis [73,74] and Soil Survey Staff [52] and the location of soil profiles within the altitudinal ranges. The capital letters refer to soil sampling areas, as shown in Figure 1. The climate of the Northern Apennines is characterized by cold winters with frequent snowfalls and hot summers. The mean annual temperature of the study areas ranges from 14.6 °C at the lowest altitude to 3.3 °C at the highest (Table S2 of the Supplementary Materials). The lowest temperatures were recorded in January (1.18 °C on average) with the exception of the highest altitude, where the lowest temperatures were recorded in February (−3.3 °C), and the warmest ones were recorded in July and August. The mean annual cumulative precipitation was 1004 mm and ranged between 753 and 1703 mm (Table S1 of the Supplementary Materials) with an uneven distribution between the months. In addition, above 2000 m a.s.l., the snow covers the soils for about 6 months per year (from November to April), which also happens at Mount Cimone [64]. The mean air temperature and mean annual cumulative precipitation data were used to determine soil moisture and temperature regimes according to the Newhall Simulation Method [72]. With the approximations allowed by the Soil Taxonomy [52], it was possible to identify a mesic temperature regime between 200 and 1300 m a.s.l. and a frigid temperature regime between 1300 to 2000 m. The soil moisture regime was found to be udic between 200 and 1800 m a.s.l. and perudic over 1800 m of altitude (Table 1). Table 1. Altitudinal range, soil temperature regime, and soil moisture regime of the phytoclimatic zones and subzones of the northern Apennine chain, Italy, according to De Philippis [73,74] and Soil Survey Staff [52] and the location of soil profiles within the altitudinal ranges. The capital letters refer to soil sampling areas, as shown in Figure 1.

Soil Profiles
To evaluate the pedodiversity, soil profiles were dug until the parent materials were sampled and then analyzed via genetic horizons and classified according to Soil Survey Staff [52]. Due to the large spatial scale of the area considered in the present investigation, the soil profiles showed large morphological variability, with differences in the genetic horizon sequence, horizon thickness, and soil profile depth. To cope with these issues, in this study, the upper 30 cm of soil profiles was considered for evaluation of the soil physical-chemical properties, and the soil horizons included within this soil depth interval were grouped into three layers: layer A, which included the A horizons; layer B, which included the AB, BA, EA, EB, Bw, and BC horizons; and layer C, which included the CB and C horizons. As a whole, the minimum, maximum, and mean thickness was 1.0, 21.0,

Soil Profiles
To evaluate the pedodiversity, soil profiles were dug until the parent materials were sampled and then analyzed via genetic horizons and classified according to Soil Survey Staff [52]. Due to the large spatial scale of the area considered in the present investigation, the soil profiles showed large morphological variability, with differences in the genetic horizon sequence, horizon thickness, and soil profile depth. To cope with these issues, in this study, the upper 30 cm of soil profiles was considered for evaluation of the soil physical-chemical properties, and the soil horizons included within this soil depth interval were grouped into three layers: layer A, which included the A horizons; layer B, which included the AB, BA, EA, EB, Bw, and BC horizons; and layer C, which included the CB and C horizons. As a whole, the minimum, maximum, and mean thickness was 1.0, 21.0, and 5.7 cm for layer A; 2.0, 23.0, and 10.1 cm for layer B; and 1.0, 27.0, and 14.1 cm for layer C, respectively.

Laboratory Analyses
The pH was determined potentiometrically in deionized water with a 1:2.5 soil:water ratio using a pH-meter electrode (Crison, Germany). The particle size distribution was determined via the pipette method described by Gee and Bauder [75], after the dispersion of soil samples with sodium hexametaphosphate solution. Organic C and total nitrogen (TOC and TN, respectively) contents were determined using a Thermo Scientific Lab EA-1110 dry combustion analyzer. Exchangeable Ca and other cations (Na, K, and Mg) and cation exchange capacity were determined through extraction with a 1 M NH 4 -acetate solution at pH 7 [76]. Cations were measured by inductively coupled plasma optical emission spectrometry (ICP-OES, Ametek, Arcos Spectro). Exchangeable Ca (Ca exch ) and the percentage of base saturation (BS) are reported in this study because of the well-known role of Ca exchange in promoting SOC storage [77] and the importance of BS for soil taxonomic classification and soil fertility.

Soil Organic Matter Stock Calculation
Organic C and TN stocks (SOC stock and TN stock, respectively) within a 0-30 cm soil depth were calculated by the sum of the SOC stock and TN stock of each horizon included within that interval depth using the equation in [78] and expressed as Mg ha −1 : where TOC or TN represents the contents of total organic carbon (%) and total nitrogen (%), BD is the bulk density (g cm −3 ), and sk is the coarse fragment content (Ø > 2 mm; %) by volume. BD was estimated according to the pedotransfer function proposed by Hollis et al. [79] and broadly used for soils located within the Italian landscapes [80][81][82][83][84]. The presence of coarse fragments was recorded during the field survey [85].

Statistical Method
To better detect the influence of the soil formation factors related to altitude, the dataset was divided into two distinct ones: a dataset that grouped the soils located between 400 and 1000 m a.s.l. (hereafter ≤ 1000 m a.s.l.) and another including soils from 1000 to 2134 m a.s.l. (hereafter > 1000 m a.s.l.). Further, a clear vegetation type change was observed at 1000 m a.s.l. Indeed, below 1000 of altitude, CHE and MF generally shared the same altitudinal range. Similarly, the vegetation types found between 1000 and 2000 m of altitude (i.e., CON, BEF, BLB, and GS), in most of the cases, shared the same altitude.
Within each altitudinal range and without considering the soil types (pedodiversity), a Kruskal-Wallis one-way ANOVA [86] was used to compare the effects of land-use and sandstone formation type on the investigated soil parameters. The Kruskal-Wallis test [86] was also used to verify the effect of soil types (pedodiversity) on soil parameters without considering vegetation cover.
Hierarchical cluster analysis (CA) using Euclidean distance and Ward grouping methods [87] was performed on the experimental data of the two altitudinal ranges (i.e., ≤1000 m a.s.l. and >1000 m a.s.l.) to cluster the collected soils in groups sharing similarities. The k-means method was then employed [87], dividing the soil samples into a predetermined number of clusters, with the smallest difference observed between elements within each cluster and the largest one between the elements of different clusters. The groups obtained by CA were processed using the Kruskal-Wallis test to identify differences in the investigated soil parameters. To identify the main drivers affecting the groups obtained by CA, principal component analysis (PCA) was carried out.
A Kruskal-Wallis test was used to evaluate the statistical diversity of OC and TN stocks among the different uses of soil and between the two altitude ranges (i.e., ≤1000 m a.s.l. and >1000 m a.s.l.).
The statistical analysis was conducted using the R software (version 4.2.1) [88].

Pedodiversity and Physicochemical Properties at Lower and Higher Altitude Range
Most of the investigated soils were classified as Inceptisols and Entisols (Table S1 of the Supplementary Materials) [52]. Alfisols at low altitude (200 m a.s.l.), developed on the ADO sandstone type (Figure 1), were also found under CHE. However, because of the scarce number of replicates for soils from 200 to 400 m a.s.l., these samples were deleted from dataset and, therefore, their data were not further elaborated. Within the ≤1000 m a.s.l. altitudinal range, soils were classified as Entisols (18%) and Inceptisols (82%). Most of the Entisols were Typic Udorthents (61%), which were mainly found within CHE, while Lithic Udorthents were less frequently found (39%) and were mainly observed under MF and within fields with slopes ≥ 30 • . Most of the Inceptisols were classified as Typic Dystrudepts (74%) and mainly recorded under CHE. About 18% of Inceptisols were Lithic Dystrudepts, with most of them under CHE, and about 8% were Oxyaquic Dystrudepts, which were found under both CHE and MF. Eutrudepts were also found (18%), all of them under MF.
Taking into account the vegetation cover, CHE presented lower pH, BS values, and contents of Ca exch , sand, and clay than MF ( Table 2). No significant differences for the amounts of TOC and TN were observed ( Table 2).
Comparing the soil subgroups, most of the investigated soils had a loam soil texture class (Table 3), with a trend toward a higher amount of sand observed among the Lithic Udorthents and Lithic Eutrudepts (sandy loam), with a higher amount of silt and lower amount of sand in Lithic Dystrudepts (silty loam) and a higher clay content and lower sand content in Typic Eutrudepts (sandy clay loam).
The soil subgroups identified at a ≤1000 m a.s.l. altitudinal range did not show differences in TOC and TN, although a higher C:N ratio was observed in Typic Dystrudepts compared to Lithic Dystrudepts and Lithic Udorthents. For the Ca exch , BS, and pH, their values were generally lower in Oxyaquic and Typic Dystrudepts than in Eutrudepts and Entisols.

Higher Altitude Range
At >1000 m a.s.l., we observed higher pedodiversity compared to that of the lower altitudinal range dataset (≤1000 m a.s.l.) ( Table 2). Soils were classified as Entisols, Inceptisols, and Spodosols (24, 53, and 8%, respectively). Mollisols were also found under mountain grassland (GS) and represented 4% of soils investigated within the >1000 m a.s.l. altitudinal range (Table 3). Specifically, Lithic Hapludolls were found under GS, showing only a mollic epipedon and no other diagnostic horizons due to the very shallow soil. Entisols were mainly represented by soils with udic and perudic moisture regimes and a temperature regime from frigid to cryic. Lithic Udorthents (27%) were found under CON and BEF, while Typic Udorthents (9%) were observed under all vegetation cover types ( Table 2). A cryic soil temperature regime was observed under GS, where Lithic Cryorthents (4%) were identified.
Most of the Inceptisols had a dystric feature, namely, a base saturation of less than 60 percent at all subsurface horizons between 25 and 75 cm below the soil surface. Furthermore, about 25% of Inceptisols belonged to Typic Dystrudepts, and about 7% were Dystrudepts featuring a lithic contact within 50 cm of the soil surface and/or with an umbric or mollic epipedon less than 50 cm thick (Lithic Dystrudepts, Humic Dystrudepts and Humic Lithic Dystrudepts). Under BEF and CON, Eutrudepts (Humic Lithic Eutrudepts and Lithic Eutrudepts) were also found.
Taking into consideration the vegetation cover, the differences observed for both soil particle size distribution and TOC content were negligible. Notably, the lowest TN content and, as a consequence, the highest C:N ratio were found under BEF. BLB, instead, presented the lowest values of both Ca exch content and BS. Finally, the highest pH values were observed under GS (Table 3).
Soils located at >1000 m a.s.l. belonged mostly to the sandy loam texture class (Table 4). Like at the ≤1000 m a.s.l. altitudinal range, no differences in TOC and TN were observed. The C:N ratio showed some differences, but no clear trend related to soil type was evident. However, as expected, Spodosols presented lower pH, Ca exch , and BS than Mollisols. Table 2. Mean ± standard deviation of pH; base saturation (BS) and C:N ratio values; and concentrations of total organic carbon (TOC), total nitrogen (TN), exchangeable Ca (Ca exch ), sand, silt, and clay at a 0-30 cm soil depth, as well as identified soil types according to Soil Survey Staff [52] of chestnut groves (CHE), mixed forests (MF), beech forests (BEF), blueberry (BLB), coniferous forests (CON), and grasslands (GS) located within 400-1000 m and 1000-2134 m altitudinal ranges of the northeastern Apennine chain (Italy). Within each altitudinal range, different letters indicate significant differences (Kruskal-Wallis test, p < 0.05) among the means.

400-1000 m Altitudinal Range 1000-2134 m Altitudinal Range
28.9 ± 28.   CA was used to cluster the soil samples collected from each horizon of the study sites located at ≤1000 m a.s.l. into three clusters ( Figure S1 of the Supplementary Materials). One cluster (hereafter called Cluster 1) included 45 statistical units mostly belonging to layer A. The largest cluster (hereafter called Cluster 2) had 145 statistical units including soil samples from layer B and layer C. The smallest cluster (hereafter called Cluster 3) mainly included soil samples collected under CHE, although the A horizons of CHE were also found. The Kruskal-Wallis test carried-out on the identified clusters showed the highest values of pH (5.9), TOC content (58.3 g kg −1 ), TN content (3.7 g kg −1 ), Ca exch content (13.6 cmol (+) kg −1 ), BS content (66.4%), and clay content (15%) for the soil samples of Cluster 1. Between Cluster 2 and 3, some further differences were observed: TOC and TN contents were higher in Cluster 3 than in Cluster 2, while Cluster 2 had higher contents of sand and silt than Cluster 3 ( Table 5). The first two principal components of PCA, performed on the dataset at below 1000 m a.s.l. explained 51.0% of the total variance ( Figure 2). PCA allowed us to identity pH, OC, TN, sand, silt, Ca exch, and BS as the main soil parameters in the separation of the clusters (Figure 2).

Lower Altitude Range
CA was used to cluster the soil samples collected from each horizon of the study sites located at ≤1000 m a.s.l. into three clusters ( Figure S1 of the Supplementary Materials). One cluster (hereafter called Cluster 1) included 45 statistical units mostly belonging to layer A. The largest cluster (hereafter called Cluster 2) had 145 statistical units including soil samples from layer B and layer C. The smallest cluster (hereafter called Cluster 3) mainly included soil samples collected under CHE, although the A horizons of CHE were also found. The Kruskal-Wallis test carried-out on the identified clusters showed the highest values of pH (5.9), TOC content (58.3 g kg -1 ), TN content (3.7 g kg -1 ), Caexch content (13.6 cmol (+) kg -1 ), BS content (66.4%), and clay content (15%) for the soil samples of Cluster 1. Between Cluster 2 and 3, some further differences were observed: TOC and TN contents were higher in Cluster 3 than in Cluster 2, while Cluster 2 had higher contents of sand and silt than Cluster 3 ( Table 5).  The first two principal components of PCA, performed on the dataset at below 1000 m a.s.l. explained 51.0% of the total variance ( Figure 2). PCA allowed us to identity pH, OC, TN, sand, silt, Caexch, and BS as the main soil parameters in the separation of the clusters ( Figure 2).

Higher Altitude Range
The CA performed on soil horizons of the study sites located at >1000 m a.s.l. produced three clusters ( Figure S2 of the Supplementary Materials). The smallest cluster (hereafter called Cluster I) had 35 statistical units and included mainly horizon A of the soils characterized by MOD and CERV sandstone types. The largest cluster (hereafter called Cluster II) had 130 statistical units and included deeper horizons. The third cluster (hereafter called Cluster III) had 54 statistical units and included the A horizons of soils developed on the other sandstone types. The highest pH value was observed in Cluster I (5.7) and the lowest one in Cluster III (4.3) ( Table 6). Both TOC and TN showed the highest contents in Cluster III (81.4 and 5.4 g kg −1 , respectively) and the lowest in Cluster II. Ca exch and BS showed the highest values in Cluster I (13.2 cmol (+) kg −1 and 88.3%, respectively) and the lowest ones in Cluster II (3.0 cmol (+) kg −1 and 26.6%, respectively). For the particle size distribution, Cluster II showed the lowest sand content and the highest clay content (57 and 11%, respectively). For the above 1000 m a.s.l. altitudinal range, the first two principal components of the PCA explained 47.1% of the total variance (Table 3). PCA showed the great influence of pH, TN, TOC, sand, silt, Ca exch, and BS on the distinctions among the clusters identified by CA (Figure 3).

Higher Altitude Range
The CA performed on soil horizons of the study sites located at >1000 m a.s.l. produced three clusters ( Figure S2 of the Supplementary Materials). The smallest cluster (hereafter called Cluster I) had 35 statistical units and included mainly horizon A of the soils characterized by MOD and CERV sandstone types. The largest cluster (hereafter called Cluster II) had 130 statistical units and included deeper horizons. The third cluster (hereafter called Cluster III) had 54 statistical units and included the A horizons of soils developed on the other sandstone types. The highest pH value was observed in Cluster I (5.7) and the lowest one in Cluster III (4.3) ( Table 6). Both TOC and TN showed the highest contents in Cluster III (81.4 and 5.4 g kg -1 , respectively) and the lowest in Cluster II. Caexch and BS showed the highest values in Cluster I (13.2 cmol (+) kg -1 and 88.3%, respectively) and the lowest ones in Cluster II (3.0 cmol ( + ) kg -1 and 26.6%, respectively). For the particle size distribution, Cluster II showed the lowest sand content and the highest clay content (57 and 11%, respectively). For the above 1000 m a.s.l. altitudinal range, the first two principal components of the PCA explained 47.1% of the total variance (Table 3). PCA showed the great influence of pH, TN, TOC, sand, silt, Caexch, and BS on the distinctions among the clusters identified by CA (Figure 3).

Organic C and Total N Stocks
SOC and TN stocks differed within altitude classes and vegetation covers (Table 7). Both SOC and TN stocks showed the lowest values at the ≤1000 m a.s.l. altitudinal range. These differences were confirmed by the vegetation covers of CHE and MF. The main covers at ≤1000 m a.s.l. showed lower stocks of SOC and TN compared to the vegetation types located at a higher altitudinal range (73.3 and 4.6, and 70.8 and 4.7 Mg ha −1 for CHE and MF, respectively). Within the higher altitudinal interval, the SOC stock did not differ, while a lower TN stock (7.2 Mg ha −1 ) was found under BEF than under BLB, GS, and CON. Table 7. Mean ± standard error of soil organic carbon stock (SOC stock) and total nitrogen stock (TN stock) in chestnut groves (CHE), coniferous forests (CON), beech forests (BEF), blueberry habitat (BLB), mixed forests (MF), and grasslands (GS) and at 400-1000 and 1000-2134 m above sea level (a.s.l.). For the vegetation types and altitudinal ranges, different letters within each column indicated significant differences (Kruskal-Wallis test, p < 0.05).

Factors Affecting Pedodiversity
The investigated soils are located in a large altitudinal range, and, as a consequence, different hydrological and temperature regimes are present (from udic to perudic and from mesic to cryic, respectively), defining the different climates and phytoclimatic zones ( Figure 1 and Table 1).
In the investigated sites, the soils were mainly characterized by weak or incipient development, mainly including Entisols and Inceptisols (Table 2), which could be attributed to the sand-based lithology, high slope gradients, and low organic matter accumulation rates promoting soil erosion processes and, therefore, weak soil development. Specifically, sandstone weathering processes generate soils with low clay and base contents and a weak soil structure [61,67,89], making them vulnerable to erosion. With organic matter serving as the major actor for the formation and stabilization of soil aggregates [90,91], such vulnerability was further boosted by low clay content, which likely limited organic carbon retention within the mineral horizons [22,92,93] and led to thin A horizons (5.5 cm, on average). However, the incorporation of organic matter within the mineral horizons was likely further prevented by acidic conditions resulting in the accumulation of organic matter on the soil surface as organic horizons [94]. The weak soil development could be also attributed to the steepness of the study sites since most of the investigated soils had a slope gradient greater than 11 • , which is considered the threshold for the occurrence of intense soil erosion [95]. Therefore, the sandstone lithology, weak soil organic matter accumulation, and high slope gradients of these mountain areas of the Apennines promoted the formation of young and poorly developed soils (i.e., Entisols and Inceptisols). Conversely, where erosion is limited, and favorable development conditions occur, such as the study sites at the highest elevations with north-northeast exposure, more intense pedogenesis occurs and is mainly driven by podzolization processes (i.e., Spodosols).
Taking into consideration the altitudinal ranges, higher pedodiversity was observed at both the order and subgroup levels in a range of >1000 m a.s.l. compared to that found in the range of ≤1000 m a.s.l. ( Table 2). Because of the role of vegetation type on soil properties and development [96][97][98], the highest pedodiversity found at >1000 m a.s.l. might be attributable to the higher vegetation cover diversity compared to that at ≤1000 m a.s.l., as discussed in the following paragraphs.
However, the highest vegetation cover diversity observed at >1000 m a.s.l. could be, in turn, due to the higher climate variability in terms of the mean air temperature, which is in accordance with several previous studies [99,100], demonstrating the great influence of climate changes related to altitude on vegetation placement within the different phytoclimatic zones. Indeed, while from about 400 to 1000 m a.s.l., the mean air temperature decreased by about 2 • C, within the >1000 m a.s.l. altitudinal range, the temperature change was about 6 • C (Table S2 of  The investigated area ranging between 400 and 1000 m a.s.l. belongs to the Castanetum phytoclimatic zone (Table 1), which has been populated since ancient times by chestnut stands [104]. Here, suitable soils for chestnut stands are present and are formed on sedimentary or siliceous soils [105] where acidic conditions predominate [106]. Acidic soils are, in fact, favorable for the growth of such plant species [22,107,108], with chestnut being susceptible to the presence of the carbonates that typically come from limestone and marl [109]. The investigated CHE soils were sampled from a mountain farm and had suitable physicochemical characteristics to produce chestnut plants with a good vegetative and adaptive state [110]. Indeed, marked dystric characteristics (Lithic, Oxyaquic, and Typic Dystrudepts) characterized by low values of pH, Ca exch content, BS, and amount of the finest soil particles were monitored under chestnut stands (Table 2). Further, CHE was mainly found on soils developed from parent materials, which were mainly poor in Ca and clay particles (i.e., GRA and MOH). Thus, under CHE, very weak soil variability was observed, mainly because of the ancient human-based selection of soils for plant chestnut stands. Unlike CHE, natural or semi-natural MFs were more evenly distributed among the different investigated sandstone formations. Thus, even if sandstone lithology limited diversity at the level of soil order, the presence of marly pelitic intercalations within the sandstone formation could generate differences at the soil subgroup level. Besides dystric features, under MF, eutric features were found (Lithic and Typic Eutrudepts) in soil developed on marl-enriched sandstones. The presence of such CaCO 3 -rich intercalation promoted the development of soils with high pH values, Ca exch content, and BS percentages (Table 3).
Independently from sandstone formations, nonparametric statistical investigations (i.e., CA and PCA) highlighted that the highest values of BS, TOC, TN, and Ca exch were observed in the A horizon, (Table 5 and Figure 2). This would suggest intense processes within this horizon related to the organic C cycle. Specifically, despite the parent material, the higher values of Ca exch content and BS found in the A horizon could be attributed to the degradation of plant residues [111], as supported by the higher TOC content [112,113]. Thus, organic C cycle occurring in the A horizon allowed us to enhance nutrient (i.e., Ca) concentrations at the soil surface in soils that developed from nutrient-poor parent materials. The A horizons, however, were quite thin, with their mean values of thickness being 5.5 cm (from 1 to 12 cm). The key role of the A horizons in improving soil properties coupled with their thin thickness highlighted the vulnerability of these soils to degradation in areas such as mountainous regions, where soil erosion processes can be intense.

Properties of Soils at the 1000-2134 m a.s.l. Altitudinal Range
The altitudinal range from 1000 to 2134 m a.s.l. included Fagetum, Picetum, and Alpinetum phytoclimatic zones. Different moisture and temperature regimes were recorded from udic to perudic and from mesic to frigid and cryic, respectively ( Table 1). The large variability in phytoclimatic zones, soil moisture, and temperature regimes further promoted pedodiversity within the >1000 m a.s.l. altitudinal range. In the northern Apennines, above the timberline and close to the top of the mountains, the primary grasslands occurred as the end of the elevation sequence of both primary heaths and secondary grasslands [114]. The highest land here includes ridges and steep slopes lashed by winds with high speed that carry away the snowpack accumulating in other areas, increasing the risk of frozen of soils with a consequent effect on the vegetation [115]. Three soil profiles were studied in the Alpinetum phytoclimatic zone of Mount Cimone under Nardetum GS [116]. These profiles were classified as Lithic Cryorthents, characterized by BS > 70% due to the presence of Ca exch on the soil exchange complex. Under Nardetum GS located at the lower altitudes of Corno alle Scale, Lithic Mollisols were found. In Mollisols, it was possible to highlight a greater thickness of the A horizon compared to that in Entisols, but the physicochemical features of these soils were very similar (Table 4). Above the tree line, shrub communities with a predominance of Ericaceae dominated the mountain heathland (BLB, Vaccinium myrtillus and Vaccinium gaultherioides Bijelow) [117,118]. BLB soils showed marked dystric characteristics that were affected by strong leaching (low Ca exch and BS), which, in some cases, resulted in spodic characteristics (Spodic Lithic Dystrudepts) at Corno alle Scale, while Spodosols (Typic Haplorthod) were found at Mount Prado and Mount Cusna. At Mount Cimone, Vittori Antisari et al. [64] investigated paleosols containing, at varying depths (but never exceeding 30 cm), a carbon-enriched horizon (Aub) that separated the top of the soil profile, called modern, from the buried soil, called ancient, showing well-expressed spodic properties. Therefore, in the present study, only modern-developed soils were taken into consideration, thereby providing a rather simple sequence of horizons (A-Bw-C).
It was unexpected to observe the highest C:N ratio of soils under BEF (Table 2). Indeed, the relatively low C:N ratio of beech residues is commonly recognized [119,120]. The high C:N ratio due to the lowest TN content found in the soils under BEF might be attributable to the fast degradation rate of such residues [121,122]. Specifically, the high rates of degradation mobilize a large amount of N, which, if not promptly uptaken by plants, can be lost through leaching.
Similar to the comparison between vegetation types, the soil subgroups identified at >1000 m a.s.l. did not show differences in soil texture (Table 4) likely due to the even distribution of such soil types on the different lithologies. Taking into consideration the investigated soil chemical properties (i.e., TOC, TN, C:N, pH, BS, and Ca exch ), although some differences could be observed among the soil subgroups, a clear trend was not identified (Table 4) due to the wide distribution of these subgroups within the study area in terms of vegetation, topography, and climate.
As observed for soils at the ≤1000 m a.s.l. altitudinal range, the CA performed for the soils at higher altitudes showed differentiation between the A horizons and the subsurface ones ( Table 6). As already observed for the lower altitude, in A horizons at higher altitudes (Cluster I and III), the C cycling facilitated improvement of soil nutrient content, increasing Ca and BS values compared to deeper horizons (Cluster II). However, unlike the comparison between the soil subgroup vegetation types, the CA highlighted the influence of lithology on the A horizon's development. Indeed, the A horizons of soils developed on MOD and CERV showed the highest values of Ca exch content, BS, and pH in accordance with the chemical composition of MOD and CERV (Table A1) compared to that of the other lithologies found at >1000 m a.s.l.

Organic C and Total N Stocks within 0-30 CM Soil Depth
According to previous studies [22,123], within the study area, greater SOC and TN stocks were found at >1000 m a.s.l. compared to ≤1000 (Table 7). Such results were confirmed by vegetation type; indeed, the lowest SOC and TN stocks were observed under CHE and MF, both located at ≤1000 m a.s.l. These findings indicate the greater influence of climate on SOC and TN contents compared to vegetation, as observed by previous studies [47,124]. Indeed, although the vegetation affected the organic matter quality [125,126], until now, the lower mean air temperatures occurring at higher elevations strongly reduced organic matter mineralization by the microbial community [127,128]. Overall, despite the higher SOC and TN stocks at >1000 m a.s.l. compared to those at ≤1000 m, it is important to highlight that our study sites featured soils highly vulnerable to organic C loss due to the aforementioned topographic and edaphic factors, namely, the low amount of clay particles and the generally high slope gradient [129][130][131][132].

Conclusions
The present study highlighted relatively young soils (mainly Entisols and Inceptisols) in the study area. The wide distribution of Entisols and Inceptisols could be attributable to the area's sand-based lithology, high slope gradients, and low organic matter accumulation rates, which make soils vulnerable to erosion processes and, therefore, do not promote soil development. Taking into consideration the subgroup level of the Soil Taxonomy classification, higher pedodiversity was found at altitudes >1000 m a.s.l. compared to those at ≤1000 m, likely due to the higher vegetation diversity and larger mean air temperature ranges at the former altitude compared to the latter.
Within the ≤1000 m a.s.l. altitudinal range, chestnut stands were present in the Castanetum phytoclimatic zone and were characterized by soils with dystric features (Lithic, Oxyaquic, and Typic Dystrudepts) and low values of pH, Ca exch content, BS, and the finest soil particles. Unlike chestnut, the mixed forests contained soils with both dystric and eutric features because such soils also developed on marl-enriched sandstone.
The lower mean air temperatures recorded at >1000 m a.s.l. compared to the lower altitudinal range promoted the accumulation of soil organic matter, which likely reduced the erosion processes and, therefore, promoted soil development. Indeed, at >1000 m, Mollisols and Spodosols were found in addition to Entisols and Inceptisols.
Taking into consideration the SOM stock, the present study demonstrated the greater role of climate (i.e., temperature) on soil organic matter compared to vegetation. Indeed, a clear difference was observed between the investigated altitudinal ranges, while the SOM stock did not show differences between the vegetation types included within each altitudinal range. However, we are aware that the results of the present study need to be interpreted with caution, as we did not compare similar vegetation at different altitudes since vegetation is affected by climate change with respect to altitude.
The spatial pedodiversity observed in the northern part of the Apennine chain in Italy could be considered a critical basis for the protection of soil resources and pedodiversity itself in mountain regions, which are characterized by high complexity related to topography, lithology, and vegetation.