The macrorelief of Romania comprises mountainous units (the Carpathians), hilly units (the Subcarpathians), plateaus, and plains. Within these major landform units, distinct regions have been delineated based on specific morphological characteristics (Posea, 2005), as shown in Figure 1. Surface lithology varies in relation to the major geomorphological units distributed around the Carpathians (A1, A2, A3), which occupy the central part of the country. Mountainous areas are characterized by a diverse assemblage of metamorphic, igneous, and sedimentary rock formations. The Subcarpathians comprise three regions (B1, B2, and B3), whose defining common feature is the alternation of depressions and hills. Lowland plains (E1, E2) are dominated by loess and loess-like deposits, whereas limestone and sandstone formations prevail in hilly (C4) and plateau regions (C1, C2, C3, D) (Figure 1).

Romania has a temperate continental climate, with a mean annual temperature of about 9.5 °C. Temperatures decline sharply with altitude, from below 0 °C in high mountains to 12 °C in the southern and southeastern regions. Mean annual precipitation averages 700 mm yr⁻¹, increasing with elevation and from east to west. The lowest values (<400 mm yr⁻¹) occur in the southeastern part, while high mountain areas receive over 1,000 mm yr⁻¹ (Dumitrescu and Bîrsan, 2015).

The spatial distribution of soils in Romania is illustrated in Figure 2. Given the high degree of correlation between RSST (2012) and WRB (2022), and in order to improve the clarity of the text, soil type names are replaced by references.

The soil cover is highly heterogeneous and exhibits a distinct concentric zonation around the Carpathian Mountains. Mountain regions are predominantly characterized by acidic soils, particularly Cambisols, Podzols, and Umbrisols, with Andosols occurring locally in volcanic areas of the Eastern Carpathians. The surrounding hilly regions and plateaus are mainly dominated by Luvisols and Phaeozems, which gradually transition into an extensive belt of Chernozems in the peripheral lowland plains (Figure 2). Fluvisols are associated with major river valleys and commonly co-occur with Gleysols and saline soils (Solonchaks and Solonetz). The relatively similar areal proportions of Cambisols, Luvisols, and Chernozems are consistent with the spatial distribution of the main landforms across Romania (Patriche et al., 2023).

The Romanian System of Soil Taxonomy (RSST, 2012) is a hierarchical classification system comprising soil classes, types, and subtypes (Florea and Munteanu, 2012). The 12 soil classes include 1–4 soil types, while the 29 soil types correspond to the Reference Soil Groups (RSGs) of the World Reference Base (WRB) (IUSS Working Group WRB, 2022). Each soil type is further divided into 2–3 subtypes, depending on map scale, which correspond to WRB qualifiers (IUSS Working Group WRB, 2022).

Soil classes group soil types and appear only in the legend; thus, soil types constitute the primary cartographic unit on soil maps. A soil type corresponds to real, natural soil bodies occupying distinct landscape areas, shaped over time by the interaction of soil-forming factors and processes (Munteanu and Florea, 2001). At lower taxonomic levels, the qualifier typic is applied when a type conforms to the standard horizon sequence. Increased complexity, such as the presence of a specific horizon, is indicated by additional subtypes.

Each polygon on the Soil Map of Romania (1:200,000) is linked to an attribute table indicating either a single soil type with one subtype or a combination of two to three soil types forming soil complexes. In this study, this structure was retained; however, on the soil type map, complex names were reduced to the dominant soil type to ensure comparability and clarity.

A cartographic scale of 1:1,000,000 has been adopted across European Union (EU) countries for geographic databases and national soil maps. In the present study, a scale of 1:750,000 was selected as a suitable solution for pedodiversity analysis across geographical regions and landform units, with 100 ha established as the minimum cartographic unit (Jones et al., 2005).

In the first stage, the vector database SIGSTAR-200, derived from the Soil Map of Romania (1:200,000), was generalized (upscaled) to 1:750,000 (Table 1). This process yielded soil type–subtype vectors consistent with the Romanian System of Soil Taxonomy (RSST, 2012) (Florea and Munteanu, 2012; Vintilă, 2018). Generalization procedures were designed to preserve, as much as possible, the pedogenetic information that characterizes soil landscapes at the given scale (Constantini and L’Abate, 2016).

Modern digital technologies, such as Remote Sensing and Geographic Information Systems, provide powerful tools for soil identification and mapping. Hence, for refining polygons of soils formed on sandy substrates (Arenosols), spectral indices were applied (Secu et al., 2022). Furthermore, the Copernicus High Resolution Water and Wetness Layers (2009–2015) have proven effective in updating soil maps in floodplain areas, offering detailed insights into hydrological regime dynamics (Copernicus, 2024).

In the second stage, the soil type–subtype vector was further generalized. Insular polygons representing soil subtypes were aggregated with the dominant polygon encompassing them, while maintaining the taxonomic consistency of qualifiers. Based on the hierarchical structure, vector data and attributes were dissolved, resulting in a coarser spatial subdivision (Schiavina et al., 2023). This procedure reduced the number of polygons and simplified soil nomenclature (e.g., Cambic and Luvic Chernozem were aggregated as Chernozem) (Figure 3). However, this process inevitably entails information loss, a limitation also emphasized in the specialized literature concerning the upscaling of spatial data (Kokkonen et al., 2006).

According to conventional approaches, Shannon’s entropy is the most widely applied index for quantifying pedological diversity, as it reflects both the relative abundance and distribution of soil classes (McBratney and Minasny, 2007). In this study, pedodiversity was calculated by considering the number of polygons and the number of taxa at two hierarchical levels, following the Romanian System of Soil Taxonomy (Florea and Munteanu, 2012), and across two spatial entities: geographical regions and landform units (Figure 3). This approach enables comparability among different territorial units and supports the identification of pedological heterogeneity both at the regional scale and at the level of landform units.

The soil type–subtype and soil type vectors were first clipped to the boundaries of geographical regions and then to landform units using ArcGIS 10.8.2 (ESRI, Redlands, CA, USA) (Figure 3). Landform units (e.g., plains) encompass multiple geographical regions (e.g., Romanian Plain, Western Plain), so the number of polygons assigned to each spatial structure varies with the size and complexity of genetic factors. To prevent aggregation errors, the surface area of each taxonomic unit in the soil type–subtype vector was kept identical to that in the soil type vector.

Shannon’s index represents a measure of information about a set of objects (e.g., soil types) with different probabilities of occurrence within a region or cartographic unit (Shannon and Weaver, 1964; Ibáñez et al., 2013). Given that two vectors of different levels of detail were employed, relative entropy, joint distribution, and the connectivity index were calculated and adapted to the requirements of this study (Yabuki et al., 2009) (Figure 3).

Through aggregation, some taxa lose information, while others retain the same number of polygons in both the soil type–subtype and soil type vectors. The latter, typically occurring in small areas, reflect paleogeographic evolution, geological structure, and related factors (Ibáñez et al., 2013). Based on this, taxa within each spatial unit were classified into two groups: variable taxa and constant taxa (Figure 3). Entropy was calculated separately for variable taxa extracted from the soil type vector (Y_A) and the soil type–subtype vector (Y_B), and for constant taxa (Y_C), which are common to both groups (Figure 3).

In this study, the ratio between cartographic units and taxonomic units was quantified using the “relative diversity index” (Pielou, 1966), also referred to as “the diversity index” (Yabuki et al., 2009), calculated using Equation 1. Y A , B , C = − ∑ p = 1 n p i ⁡ ln ⁡ p i ln n (1) where p_i is the proportion of polygons covered by the i-th pedotaxon, and n is the number of taxa within each spatial unit. Y_A represent the spatial distribution of polygons for soil types, Y_B represent the dispersion of polygons for soil type–subtypes, while Y_C denotes the joint distribution of polygons for both soil types and subtypes, generally corresponding to rare and endemic soils. Relative diversity equals 0 when the number of polygons assigned to taxa within a territorial unit is non-uniform and reaches 1 when the distribution is equiprobable (Yabuki et al., 2009; Lo Papa et al., 2011).

The mutual information between two variables (A and B) is associated with joint entropy and was employed to indicate simultaneous distribution (Ibáñez and De Alba, 2000; Yabuki et al., 2009). In this study, the number of polygons corresponding to soil type–subtypes that share the same distribution as soil types constituted the joint entropy, which was calculated using Equation 2. S A , B = − ∑ i = 1 n ∑ k = 1 n p A i , B i ln ⁡ p A i , B i / ⁡ ln n (2) where A i B i   represents the probability of the simultaneous distribution of polygons for soil type–subtypes and soil types, and n is the number of taxa common to both groups.

In the RSST (2012), soil nomenclature distinguishes two groups of qualifiers: the typical qualifier (haplic in the WRB), indicating the core soil concept, and additional qualifiers (principal and supplementary in the WRB), reflecting variability in pedogenetic factors and anthropogenic influences. Polygons defined by the typical qualifier on the higher-resolution map correspond to soil types on the generalized map, so their overlap represents simultaneous spatial distribution, with probabilities ranging from 0 to 1.

To evaluate soil map connectivity at two levels of detail, the connectivity index (R), originally applied to soils and land use by Yabuki et al. (2009), was adapted and calculated as follows: R A , B = 1 2 Y C × S A , B × Y A + Y B (3) where Yc represents the entropy of constant taxa, S_(A,B) denotes the entropy of simultaneous distribution, and Y_A + Y_B is the sum of entropies of variable taxa. The connectivity index R_(A,B) approaches 1 when the number of polygons associated with variable taxa does not differ significantly between the soil type–subtype and soil type vectors and tends toward 0 when the differences are substantial.

The aggregation of taxonomic units from soil type–subtypes to soil types was performed in ArcGIS 10.8. The number of polygons for each of the two spatial structures was obtained using the Patch Analyst module (Rempel et al., 2012), which enables the quantification of the spatial configuration of cartographic units.

Descriptive statistics, Wilcoxon Signed-Rank Test, and the Shapiro–Wilk test for polygon distribution normality within each spatial entity were processed in SPSS 27 (International Business Machines Corporation, IBM). Database organization and index calculations were conducted in Microsoft Excel.

To assess the impact of generalization on the polygon populations across different relief units, Wilcoxon Signed-Rank Test was employed, as it is suitable for comparing related samples. The null hypothesis (H₀) states that the median difference between values before and after generalization is zero, while the alternative hypothesis (H₁) states that the median difference is non-zero, indicating a significant effect of reducing the number of polygons from soil type–subtype to soil type across geographical regions and landform units.

In addition, similarity and differences among diversity indices, joint distribution, and the connectivity index across the 15 geographical regions were assessed using bivariate regression, implemented in OriginPro 2025 (OriginLab Corporation, Northampton, MA, USA). This statistical approach enabled the identification of functional relationships among indices and the assessment of their spatial comparability.

For Romania, the soil type–subtype vector contains 1,878 polygons, whereas the simplified soil type vector includes 1,232 polygons. The hierarchical proportions shift between these taxonomic levels. At the type–subtype level, five major groups—Luvisoluri, Cernoziomuri, Eutri- and Districambosoluri, Aluviosoluri, and Faeoziomuri—account for 62.2% of polygons. After aggregation at the soil type level, the hierarchy changes, with Eutri- and Districambosols, Luvosols, Phaeozems, Luvisols, and Chernozem representing only 51.8% (Figures 2, 4 for WRB correspondence). This demonstrates that taxonomic simplification reduces both the number of polygons and the relative balance among main soil types, affecting pedodiversity assessment and regional-scale interpretation.

Geographical regions and landform units show a unidirectional decrease in the number of polygons from soil type–subtypes to soil types (Table 2). For soil type–subtypes, the number of polygons is higher in lowland areas compared to mountain units, whereas the hierarchy is reversed when considering soil types (Table 2).

The Wilcoxon Signed-Rank Test indicates differences in the polygon populations after aggregation for landform units (Table 3). A similar pattern was observed for geographical regions, with the exception of the Mehedinți Plateau.

When comparing the number of polygons at the soil type–subtype and soil type taxonomic levels to statistical maxima, disparities in the hierarchies formed among landform units and geographical regions are revealed. However, a common trend persists: the top position is consistently occupied either by lowland or by mountain units. For the type–subtype level within landform units, the maximum number of polygons follows the sequence: plains > Carpathians > Subcarpathians > plateaus. By contrast, for geographical regions, the hierarchy differs slightly: plains > Carpathians > plateaus > Subcarpathians (Tables 2, 3).

These differences reflect the distinct ways in which landscape structure and accessibility influence the taxonomic detail of mapping. In plains, high accessibility and favorable fieldwork conditions enable better characterization of soils at the subtype level, thereby increasing the number of polygons. In mountain areas, although accessibility is more limited, the diversity of parent materials and relief variability confer a numerical advantage to soil types compared to plains.

The Shapiro–Wilk test (95% confidence interval) applied to the number of polygons in each region indicated a normal distribution. The number of polygons within the four landform units also followed a normal distribution, as the data originated from relatively homogeneous geographical regions.

Excessive skewness at the soil type–subtype level, however, indicated an asymmetric distribution of polygon counts. In the plains and Subcarpathian units, this imbalance was driven by the large number of polygons belonging to dominant soils, namely, Chernozems (183) and Luvisols (114).

Shannon’s Evenness Index, or diversity index (Y_A,B,C) show similar mean values for geographical regions (0.92, 0.90, 0.86), with very low variance (<0.003), suggesting a limited dispersion of values around the mean. In this context, Y_(A,B,C) does not reasonably differentiate spatial distribution in terms of polygon counts and, when considered individually, is therefore not suitable for ranking spatial units.

By contrast, the entropy of simultaneous distribution exhibited a wider range of values compared to the Evenness Index (0.25–0.88) (Figure 5). Geographical regions with a small number of variable taxa (≤3), such as the Dobrogea Plateau (C2) and transitional units from mountainous to plateau or plain landscapes (e.g., Getic Subcarpathians - B3 and Moldavian Subcarpathians - B1), record high values of simultaneous distribution (0.88, 0.66, 0.59), lying in the upper whisker interval (Figures 1, 5).

Mixture analysis facilitated the separation of connectivity index values for soil type–subtype and soil type taxonomic levels into three categories: weak (<0.25), moderate (0.25–0.45), and strong (>0.45) connected regions. Weak connectivity between two taxonomic levels is characteristic of mountain regions (Southern and Eastern Carpathians), the plains (Western Plain), and the plateau region (Transylvanian Depression) (Figure 6a). Low connectivity indicates spatial diversity, reflecting the participation of numerous subtypes in the formation of the soil cover or the presence of soil complexes, particularly in high mountain areas (Southern Carpathians). The connectivity index (RA,B) reached its highest value in the Dobrogea Plateau (0.72), located in southeastern Romania, the country’s most arid region. This was followed by transitional regions (Getic Subcarpathians, Moldavian Subcarpathians) situated more centrally within the country (Figures 1, 6a). The third group comprises regions with diverse physical–geographical characteristics and sizes (Romanian Plain, Danube Delta, Curvature Subcarpathians, and Western Carpathians), with an average connectivity index of 0.36 (Figure 6a).

The connectivity index produced a hierarchy of landform units that differed from that of geographical regions. Specifically, the sequence of landform units based on the connectivity index was as follows: plains (0.17) < plateaus (0.23) < Carpathians (0.27) < Subcarpathians (0.28) (Figure 6b). A common feature of some regions (e.g., B1, B3) and units (B) with weak connectivity was the small number of variable taxa (B1 = 3, B3 = 3, B = 5) forming the variable group of soil types and subtypes. For most taxa, polygon counts showed minimal differences between the two taxonomic levels, indicating a relatively homogeneous structure, except for Luvisols in the Subcarpathians and Cambisols in the Western Hills.

The relationships among diversity indices, joint entropy, and the connectivity index at the regional level were modeled using bivariate correlation. Joint entropy showed a very strong association with the connectivity index (0.98). The diversity index of variable taxa at the type–subtype level had a weak influence on both joint entropy (0.39) and the connectivity index (0.46) (Figure 7). Constant taxa formed a separate group, exhibiting a weak relationship with the diversity Index of variable taxa, represented by soil types (0.22) (Figure 7).

The database used in this study consists of the set of taxa belonging either to geographical regions or to landform units (Figure 1). Each taxon is represented by a variable number of cartographic polygons (Table 1). Through the upscaling process, this number may remain unchanged or may decrease, depending on the taxonomic level and the complexity of the pedological landscape structure.

The difference in the number of polygons between soil type–subtypes and soil types becomes zero only when the spatial distribution of a taxon is identical at both taxonomic levels. The group of constant taxa (C) was obtained by selecting, for each region and unit, those taxa that maintain an identical number of polygons following aggregation (Figure 3). The number of constant taxa varies across regions, with polygons assigned to each taxon ranging from a minimum of 1 to a maximum of 9. Any taxon classified at the soil type–subtype level consistently exhibits a greater number of polygons than at the soil type level (Table 2; Figure 4). The group of variable taxa was divided into two numerically balanced subroup for each region and unit.

Based on previous observations, data overlap is observed for constant taxa. The concept of overlap is widely discussed in related fields. Some modeling approach disregards data from “overlapping regions” and focuses only on “non-overlapping regions” (Xiong et al., 2010). The term “overlapping region” is used in a generic sense, applied to two-dimensional datasets (Santos et al., 2023), and should not be confused with an actual geographical region. In this context, the “non-overlapping regions” define the group of variable taxa and correspond to the sum of the entropies of the variable taxa (Y_A and Y_B) as expressed in Equation 3. The “overlapping region” corresponds to the group of constant taxa and is represented by the entropy of the constant taxa (Y _c ) in Equation 3 (Figure 3).

Unlike the method proposed by Yabuki et al. (2009), this study omitted the grid approach, since the surface area allocated to each taxon does not change through aggregation. For the same reason, the sum of entropies of variable taxa does not employ a fraction with numerator 1. The spatial diversity of variable taxa (Y_A, Y_B) across geographical regions did not reveal significant differences between the two groups, suggesting limited variability in entropy from one taxon to another (Figure 6). When soil areas remain constant, pedodiversity should not substantially increase when moving down the taxonomic scale (Minasny and McBratney, 2007). The results obtained here are consistent with this hypothesis, confirming that the transition from soil types to subtypes does not necessarily imply a significant increase in spatial diversity.

When diversity is interpreted in the sense proposed by Yabuki et al. (2009), constant taxa appear to significantly influence the overall level of spatial diversity, even if the information is redundant in the aggregation process. The entropy of constant taxa is incorporated in Equation 3 as a factor that equally influences the sum of entropies of the variable taxa.

The rarefaction analysis module in PAST 5.0 was used to estimate the number of taxa based on the number of polygons. Richness (S) was tested for both a small number of individuals and a large number of individuals (Colwell et al., 2012). The analysis focused on geographical regions and relief units, respectively mountainous and plain areas selected for modeling. The number of taxa estimated from polygon counts was equal to that considered in this study for the Southern Carpathians (14) and slightly lower for the Western Plain, with values of 13 for soil types and 12.6 for soil type–subtypes. Non-significant differences between the number of taxa extracted from the maps and the estimated values for each spatial unit support the reliability of the soil maps.

Aggregation considerably reduces the degree of variation within the analyzed area and is inherently associated with information loss (Kokkonen et al., 2006). At the working scale applied here, aggregation could not be performed for the Mehedinți Plateau (C5), which comprises only five taxa, four of which were represented by a single polygon. For this reason, it was merged with the Getic Plateau (C3), based on the principle of equivalent rank, at the level of geographical regions (Figure 6a).

The results indicate that the connectivity index tends to reach high values in small geographical regions, such as the Dobrogea Plateau, Getic Subcarpathians, and Moldavian Subcarpathians, which are characterized by a low number of variable taxa (2, 3, and 2, respectively). Interestingly, the number of variable taxa in a region does not appear to directly influence the value of the connectivity index. For example, although the Getic Subcarpathians and the Danube Delta both contain three variable taxa, joint entropy is lower in the former (0.48) compared to the latter (0.59), which substantially affects the connectivity index. This finding confirms observations from the literature that probability measures alone are not always sufficient to explain the complexity of spatial structures and must be complemented by indices that capture interactions and informational redundancy (Baez et al., 2011). Consequently, the connectivity index emerges as a valuable tool not only for describing pedological diversity but also for understanding the structural stability of spatial relationships in regions with contrasting characteristics.

A strict comparative analysis of connectivity levels between cartographic units across different taxonomic levels is not yet possible due to the lack of similar studies. Nevertheless, the results obtained here allow the identification of common features of entropy applied to the soil cover. Consequently, the connectivity index is lower in Carpathian regions characterized by higher altitudes (A2 and A1) than in lower-altitude regions (Western Carpathians), indicating that altitude and morphological complexity enhance spatial diversity in high-elevation mountain areas (Figure 6a). Similarly, Shannon pedodiversity is higher in high-altitude mountain regions than in lowland areas, reflecting the strong influence of biodiversity, climatic diversity, and topographic variability, which together increase the number of soil types per unit area (Ibáñez and Effland, 2011; Ibáñez et al., 2013; Pascual-Aguilar et al., 2015; Constantini and L'Abate, 2016; Fu et al., 2018). The Eastern Carpathians exhibit an extensive spatial distribution of volcanic and sedimentary rocks associated with Andosols, Eutric Cambisols, and Luvisols. In the Southern Carpathians, crystalline schists dominate, and altitudes above 2000 m, together with glacial landforms, promote the development of Podzols and Umbrisols. These conditions contribute to higher pedodiversity compared to the Western Carpathians, where only a few peaks exceed 1800 m and Cambisols and Luvisols prevail (Figures 2, 6a).

The greatest dispersion of connectivity index values was recorded for the plateau regions. The low connectivity index observed for the Transylvanian Depression (0.30) suggests a low probability that the polygons associated with variable taxa (Y_A and Y_B) display similar spatial distributions, in contrast to the Dobrogea Plateau, which exhibits a higher index (0.88) (Figure 6a). These pedodiversity contrasts are largely driven by the distinct climatic characteristics of the two regions. The Transylvanian Depression is characterized by a temperate-continental climate, whereas Dobrogea is among the driest regions of Romania (<500 mm/year), with extreme temperatures and low precipitation levels, which result in the lowest pedodiversity (Minasny et al., 2010; Manea et al., 2024; Serban and Maftei, 2025).

The connectivity index places mountain regions (A1, A2) and the plain region (D1) in close positions (Figure 6a). The Southern Carpathians, Eastern Carpathians, and Western Plain share lower values of joint entropy (0.25, 0.31, 0.27) and the diversity index for constant taxa (0.87, 0.78, 0.84) compared with other regions. Guo et al. (2003) reported comparable values of taxa richness and Shannon’s diversity index for soils developed on varied parent materials, despite differing climatic conditions in the Western and Central Southern regions of the USA.

More recent studies further confirm that in certain landscapes—particularly where alluvial and aeolian landforms predominate—parent materials exert a stronger influence on pedodiversity than other environmental factors (Luo et al., 2021). In the Western Plain, the genesis of Fluvisols, Vertisols, and Solonetz is primarily influenced by parent materials and these soils are classified among the variable taxa, whereas Arenosols display a more uniform spatial distribution and are assigned to the constant taxa group. This analysis highlights the complex interaction between landforms and parent materials: while mountain regions exhibit higher diversity due to morphological variability, in plains, the control exerted by parent material becomes the determining factor in structuring pedodiversity.

Typically, larger territories are characterized by greater heterogeneity compared to smaller ones, but this feature also depends on geographical location (Ibáñez et al., 2014). Nevertheless, the Romanian Plain shows a higher connectivity index (0.38) for the two taxonomic levels than the Western Plain (0.27), even though the former covers a larger area (Table 2). The surface-connectivity mismatch can be explained by the uniformity of the parent materials in the eastern Romanian Plain, represented by aeolian loess deposits and sandy loess (Gherghina et al., 2006; Jipa, 2014). The lithology associated with the highest degree of aridity (Bondoc and Prăvălie, 2015) limits the evolution of chernozems towards a greater diversity of subtypes.

Entropy represents the average amount of information transmitted from the source and received at the destination (Lombardi et al., 2016). In this study, the information corresponds to the number of polygons in the two spatial entities (units and regions), before and after aggregation. On the other hand, Shannon’s theory is based on the selection of a message from a set of possible messages, where “selection” is a statistical concept reflecting the constraint in choosing a particular situation (Cole, 1993).

Considering the soil type–subtype vector as the source and the soil type vector as the destination, the signal transmitted by the intermediate product of the diversity index with joint entropy indicates higher pedodiversity in the Western Plain (0.11) compared to the Romanian Plain (0.17) (Equation 3). In the first region, both the soil type–subtypes and the soil types encompass three taxa with a higher number of polygons, compared to only two taxa in the Romanian Plain. The diversity index for constant taxa is lower in the Western Plain (0.85), reflecting a more varied data structure than in the Romanian Plain (0.90). Conversely, the sum of the diversity index (Y_A + Y_B) for variable taxa indicates a more uniform polygon distribution in the Western Plain (1.9) than in the Romanian Plain (1.7).

These results indicate that discrepancies may occur between pedological reality and the statistical interpretation of pedodiversity. For instance, describing a region using a single diversity index value limits the ability to capture the spatial complexity of soil landscapes. As shown in other studies, using a single global indicator to characterize the spatial structure of pedodiversity can obscure local variations and scale-dependent relationships (Vasat et al., 2023). These limitations point to a broader conceptual issue: pedodiversity has been regarded by some authors as a “controversial concept” (Jie et al., 2001), due to methodological constraints and the difficulty of simultaneously integrating pedogenetic factors, spatial variability, and taxonomic classifications.

The connectivity index indicates a higher degree of independence among variable taxa in plains (0.17) compared to other landform units, a pattern that is also evident at the geographical region level (Figure 6a). The joint distribution of soils across the four relief units is consistent with the connectivity index results, following the sequence: plains (0.22) < plateaus (0.27) < Carpathians (0.36) < Subcarpathians (0.39), and shows a strong correlation (0.97) between the two measures (Figure 6b). The greater number of taxa and polygons within landform units, relative to geographical regions, leads to a narrower range of connectivity index values.