Cadmium in Soils Cultivated with Avocado in Peru: Interaction with Soil Variables and Its Relationship with Parent Material

1. Introduction

Heavy metals, such as cadmium (Cd), are naturally present in soils, with concentrations varying depending on their geological origin [1]. However, anthropogenic activities, including mining, intensive agriculture, and industrial processes, have markedly elevated these levels, raising concerns about their impacts on ecosystems and public health [2,3]. Among these elements, Cd is particularly notable for its high toxicity, strong bioaccumulation potential, and environmental persistence [4], posing a significant risk when it enters the food chain [5].

In the Peruvian context, avocado (Persea americana) cultivation has become a consolidated and important agricultural export product [6]. The presence of Cd in soils is a growing concern. Recent studies in other producing areas have reported high concentrations of this metal, even in areas with no obvious anthropogenic sources [7], suggesting that the parent material could be the main determining factor. This is particularly relevant, considering that plants can absorb Cd and accumulate it in their fruits [8], which affects both food safety and the commercial value of the crop.

The spatial distribution of heavy metals in agricultural systems is influenced by a complex set of factors, including edaphic properties (pH, texture, cation exchange capacity), geological features, and management practices [9,10]. To analyze these patterns, methodological approaches have evolved from classical spatial interpolation techniques [11] to advanced tools, such as the geodetector, which allows for the assessment of interactions between multiple environmental variables [12]. Nevertheless, in tropical soils cultivated with avocado, significant knowledge gaps remain, particularly regarding the mechanisms controlling Cd accumulation in soils with basic pH, the role of underlying lithology and geological age on its spatial distribution, and the interactions between Cd and other elements such as Zn, Pb, and Na under tropical soil conditions [12,13].

This study aims to fill these gaps through an innovative approach that combines geochemical analysis through Inductively Coupled Plasma Mass Spectrometry (ICP-MS), multivariate statistics (PCA, Spearman correlations), and advanced spatial modeling (geodetector). It represents the first systematic effort to understand Cd accumulation patterns in avocado soils in Peru, integrating for the first time geochemical, edaphic and spatial approaches in a tropical context, contributing to the broader understanding of Cd sources, distribution, and associated risks, which has gained significant relevance for the scientific community and policy makers [14,15].

The specific objectives include characterizing the distribution patterns of Cd in soils from avocado plantations across different regions of Peru, identifying the most influential edaphic and geological factors in its accumulation, and establishing criteria for risk zoning based on the parent material. The relevance of this research lies in its potential to inform sustainable soil management strategies for high-value commercial crops, thereby contributing to both food security and the competitiveness of Peruvian agricultural exports. The findings will be particularly valuable for designing more efficient sampling protocols, prioritizing areas for remediation programmes, and developing technical guidelines for site selection in future avocado plantations.

2. Materials and Methods

2.1. Study Area

A total of 410 soil samples were collected from agricultural areas dedicated to avocado (Persea americana) cultivation in the departments of Ancash, Apurímac, Ayacucho, Cusco, Huancavelica, Ica, La Libertad, and Lima. The samples represent a wide diversity of soil and climatic conditions, vegetation coverage, soil taxonomic classes, types of parent material, including residual igneous, sedimentary, and transported soils (alluvial, eolian, fluvial, glacial).

2.2. Soil Sampling

We adopted the soil sampling model proposed by Havlin et al. [16], in which each soil sample represents an agricultural unit with similar crop, topographic, and agronomic management characteristics. The sampling units (SU) were defined according to different boundaries. The slope boundary refers to the distinction between hillside and flat areas. The soil boundary involves identifying areas with different textural classes and soil colours. The crop boundary refers to areas with different rootstocks and planting ages. Finally, the area boundary sets a maximum of 10 hectares per soil sample. After identifying the SU, the sampling design was implemented, considering five avocado trees per homogeneous lot, with the trees being representative and randomly distributed throughout the plantation. Finally, the collection of soil samples consisted of one sub-sample of soil per tree taken within the canopy projection at a depth of 30 cm. A total of 410 soil samples were collected from different avocado plantations in Peru.

2.3. Soil Analysis

Soil samples were analyzed at the Soil, Water, and Foliar Laboratory Network of the National Institute for Agrarian Innovation (LABSAF-INIA). Prior to physicochemical analysis, samples were pre-treated by air-drying at temperatures below 40°C and then sieved to obtain a fraction smaller than 2 mm, following procedures established by the International Organization for Standardization [17]. The characterization analysis included the following parameters and reference methodologies: Particle size distribution (sand, silt, and clay percentages) using the Bouyoucos density meter method [18]; pH in a 1:1 soil-to-water ratio [19]; Electrical conductivity (EC) in a 1:5 soil-to-water aqueous extract [20]; Calcium carbonate equivalent by acid neutralization method [18]; Organic matter, total organic carbon, total carbon and total nitrogen of the soil, determined using a LECO CN828 elemental combustion analyzer; Available phosphorus (P) for neutral and acidic soils using the Bray and Kurtz method [18]; Available potassium (K), determined following Bazán Tapia [21]; Exchangeable bases concentration (Ca²⁺, Mg²⁺, K⁺, and Na⁺) extracted with ammonium acetate, and exchangeable acidity (H⁺, Al³⁺) extracted with potassium chloride [18]. The effective cation exchange capacity (CECe) was calculated as the sum of exchangeable cations. In addition, bulk density (BD) and total CEC at pH 7 were obtained from the global digital soil mapping system SoilGrids [22], with a spatial resolution of 250 m. These data were extracted for the 0–15 cm and 15–30 cm depth intervals by downloading a file in TIFF format and extracting the values in Qgis.

2.4. Extraction and Processing of Geospatial Variables

In addition to the physicochemical parameters determined through laboratory analysis, a set of potential drivers influencing the spatial distribution of cadmium concentrations was selected based on a literature review and relevant research. These include the geoform type (GT), obtained from the 1:250,000-scale geomorphological map [23], which represents broad relief units that influence the processes of element accumulation or runoff, being compatible with the regional scale of the study. Parent material (PM) and geological age (GA), both obtained at a 1:100,000 scale from the national geological map (from sheet 21h to 27l), provide the necessary detail to capture key lithological variations. Soil type (ST), at a 250 m resolution (approximately 1:250,000) from SoilGrids [22], accurately reflects general soil conditions well-suited for regional analysis. The life zone (LZ), at 1: 100,000 from the life zone map [24], summarizes broad ecological factors such as climate and altitude, while the climate classification (CLI) obtained from the climate distribution map of Peru [25], although more general (1: 400,000), provides a macro view of the climatic environment. The altitude (ALT), with 30 m resolution, serves to refine spatial interpretations and compensate for the generalization inherent in other layers.

The selected spatial scales were adjusted to match the territorial dimension of the study, allowing for a clear representation of the environmental processes that may influence cadmium distribution. Data were obtained from reliable sources and at appropriate levels of detail, given the nature of the phenomenon under investigation, thereby minimizing potential spatial inconsistencies. While some variation exists in the accuracy of the geospatial layers used, this is not expected to significantly impact the results, as the chosen resolutions are sufficient to properly identify the general patterns of Cd distribution and the associated explanatory factors.

2.5. Univariate Statistical Analysis

Univariate descriptive statistical analysis was conducted using R software (version 4.4.1). Various tools and libraries were utilized to generate graphs and support the statistical analysis of soil properties. The dataset was initially imported directly from the clipboard using the read.delim() function. Numerical variables were selected using the select(where(is.numeric)) function from the dplyr package. Descriptive statistics were then calculated for each variable, including the mean, standard deviation, variance, coefficient of variation, skewness, kurtosis, minimum and maximum values, 25th (Q1) and 75th (Q3) percentiles, median, and the p-value significance from the Shapiro-Wilk normality test. These statistics were obtained using the summarise(across()) function and subsequently rearranged into long and wide formats using functions from the tidyr package to facilitate presentation and interpretation.

2.6. Multivariate Statistical Analysis

Non-parametric bivariate correlations among numerical variables were evaluated using Spearman’s correlation coefficient, implemented via the rcorr() function from the Hmisc package. Spearman’s correlation coefficient (r) does not assume normality in the distribution of variables or linearity in their relationships, making it particularly suitable for the typically heterogeneous and skewed nature of soil data. The correlation matrix (r) and the associated p-value matrix were generated by transforming the numerical data into matrix format. For visualization, the corrplot function was employed. Only statistically significant correlations (p < 0.01) were displayed, omitting those without significance.

To reduce dimensionality and explore multivariate patterns among the evaluated edaphic variables, Principal Component Analysis (PCA) was performed using the FactoMineR package, with visualization supported by factoextra. Previously, categorical variables were excluded, and only standardized numerical data were retained. Categorical study factors were converted into factor-type variables to facilitate visualization and grouping purposes. The analysis included the assessment of explained variances for each principal component (eigenvalues), as well as the individual contributions of the variables. The quality of representation (cos²) and the distribution of individuals within the factorial space were plotted, considering groupings according to the previously defined factors. Biplots with confidence ellipses were also generated to explore potential latent structures associated with soil classes. All graphical analysis was developed with ggplot2.

2.7. Application of the Geodetector Algorithm Principle

The Geodetector (GD) is a statistical method designed to evaluate whether a geographical factor, such as soil type or parent material, influences the spatial distribution of a dependent variable, such as total cadmium concentration in soil. It operates under a fundamental principle: if two variables (i.e., Cd and geoforms) exhibit similar geographical patterns, they are likely associated. A key advantage of the Geodetector is that it does not rely on strict assumptions about the data, enabling the exploration of interaction effects between two combined factors on the variable of interest. In this study, the Geodetector model developed by Wang et al. [26] was applied to analyze spatial heterogeneity in cadmium concentrations and to understand how certain factors influence this pattern. This approach, based on variance partitioning statistics, quantifies the degree of spatial association between a continuous dependent variable and multiple categorical factors. The model is articulated in four main components: the factor detector, ecological detector, risk detector, and interaction detector. The factor detector assesses whether a given factor is associated with the spatial distribution of a continuous variable by analyzing whether they share a similar spatial pattern that suggests a possible cause-and-effect relationship. For each factor analyzed, a power factor ‘q’ is determined that quantifies its influence on the dependent variable, according to:

$$\mathit{q}=1-{\displaystyle \frac{1}{\mathit{n}{\sigma }_{\mathit{H}}^2}}{\sum }_{\mathit{i}=1}^{\mathit{m}}({\mathit{n}}_{\mathit{D},\mathit{i}}{\sigma }_{{\mathit{H}}_{\mathit{D},\mathit{i}}}^2)$$

(1)

Where D represents a categorised driving factor; m is the number of categories of the factor F; H determines the spatial distribution of the dependent variable; q is the power of the driving factor (D) over H; n represents the total number of samples and ${\sigma }_{\mathit{H}}^2$ is the total variance of H in the study area; ${\mathit{n}}_{\mathit{D},\mathit{i}}$ represents the number of samples in the i-th layer of D and ${\sigma }_{{\mathit{H}}_{\mathit{D},\mathit{i}}}^2$ is the variance of H over the i-th layer of D. The parameter q is quantified in the range of [0;1] [12], where a higher q value indicates a greater influence [26].

The ecological detector enables the assessment and comparison of the influence of determinants to identify which factor has the most significant effect on the spatial distribution of the analyzed variable. This process is performed in a bivariate way using the F-test, according to:

$$\mathit{F}={\displaystyle \frac{{\mathit{n}}_{\mathit{C},\mathit{p}}\left({\mathit{n}}_{\mathit{C},\mathit{p}}-1\right){\sigma }_{\mathit{C},\mathit{m}}^2}{{\mathit{n}}_{\mathit{D},\mathit{p}}\left({\mathit{n}}_{\mathit{D},\mathit{p}}-1\right){\sigma }_{\mathit{D},\mathit{m}}^2}}$$

(2)

Where F is the value of the F-test, $n_{C,p}$ and $n_{D,p}$ determine the number of samples of the determinant factors C and D within a stratum of the categorisation (p). ${\sigma }_{C,m}^2$ and ${\sigma }_{D,m}^2$ denote the variances reflecting the dispersion of the values of factors C and D.

The risk detector evaluates differences in the mean values of a dependent variable across the categories (or strata) of a given factor using a t-test. As the magnitude of these differences increases, so does the factor’s influence on the dependent variable. The interaction detector analysis examines how two factors (x₁ and x₂) jointly influence cadmium accumulation in soils, identifying whether their interaction weakens, strengthens, or operates independently of each other. This is assessed by comparing q(x₁ ∩ x₂) with q(x₁) and q(x₂). If q(x₁ ∩ x₂) is less than min[q(x₁), q(x₂)], the interaction is non-linearly weakening. If it lies between min[q(x₁), q(x₂)] and max[q(x₁), q(x₂)], it is unilaterally weakening. If it exceeds max[q(x₁), q(x₂)], it is considered bi-directionally enhancing. If it is greater than q(x₁) + q(x₂), the two factors are deemed to interact non-linearly and independently.

2.7.1. Factor Data Processing

Based on findings from previous studies and considering the availability and accessibility of relevant data, a total of 18 factors were selected as potential determinants of the ecological risks associated with total cadmium (Cd) content in avocado plantation soils in Peru. The evaluated soil properties included: soil type, particle size distribution (percentages of sand, silt, and clay), pH, organic matter content, available phosphorus (P), exchangeable sodium percentage (ESP), exchangeable calcium percentage (ECP), exchangeable potassium percentage (EPP), exchangeable magnesium percentage (EMP) and total concentrations of lead (Pb), magnesium (Mg), and zinc (Zn).

Geomorphological (GT), geological (GA and PM), edaphic (ST), ecological (LZ), and climatic (climate distribution) variables were also incorporated into the analysis. Soil properties were obtained through laboratory analyses, while the remaining variables were derived from spatial databases in vector and raster formats, which were subsequently converted into polygons. These spatial variables, used in the Geodetector analysis, were processed using geoprocessing tools in QGIS 3.40 (Bratislava). All spatial layers were clipped to the area of interest and reprojected to the WGS 84 / UTM Zone 18S coordinate reference system to ensure spatial consistency. To reduce spatial noise, small or fragmented polygons were removed. Before categorization, the “dissolve” tool was used to group values that shared common attributes.

The driving factors were categorized using expert knowledge, national metrics, and classification algorithms such as Jenks natural breaks and quantile ranking. The application of the classification method was based on enhancing the magnitude of the “q” statistic [27]. Finally, the "Point Sampling Tool" plugin was used to extract the values of each variable at the sampling points, ensuring that each point was assigned a single value per variable, ready for analysis in Geodetector.

2.8. Statistical Analysis of Non-Parametric Comparison

Differences in Cd concentrations between groups defined by the factors ST, GA, and PM were evaluated using non-parametric statistical tests, given the absence of assumptions of normality and homogeneity of variances. The Kruskal-Wallis test was initially used to identify overall differences in the median total Cd content across the groups. When significant differences were detected (p < 0.05), post hoc multiple comparisons were conducted using Dunn’s test with Bonferroni correction to control for Type I error, implemented via the dunnTest() function from the FSA package. These analyses enabled the identification of specific group pairs exhibiting significant differences. To illustrate these differences, significance letters were assigned to each factor level using the multcompLetters() function from the multcompView package, based on the adjusted p-values. The distribution of Cd by group was visualized using box plots enhanced with individual data point dispersion (jitter), generated with the ggplot2 package.

2.9. Cartographic Mapping of Influencing Factors in the Spatial Variability of Cd

Based on the results of the statistical analyses, thematic maps were generated to depict the categories within each factor that exhibited significantly different values in cadmium distribution patterns. Map production was carried out using QGIS 3.40 (Bratislava). Attribute filters were applied to isolate significant classes, and a differential categorized symbology was defined between variables. Ambiguous scales were avoided, and contrast with the background was prioritized to enhance interpretability.

All maps adhered to the principles of thematic cartography, including a graphic scale, explanatory legend, and title. A minimum cartographic unit (MCU) was established to avoid displaying polygons with minimal areas, which would hinder visual interpretation. Contiguous polygons sharing the same category were generalized into a single unit when thematic accuracy was not compromised. Particular attention was given to maintaining a clear visual hierarchy, minimizing background graphic load while emphasizing significant areas. Accordingly, maps were produced for the Lima region using the factors with the highest explanatory power (q > 15% according to Geodetector): parent material (PM), geological age (GA), geoform type (GT), soil type (ST), climate (CLI), and life zone (LZ).

The categories of each factor were weighted according to their respective q-values from the Geodetector analysis and subsequently combined using an additive index, spatially overlaid to identify areas with different levels of risk. This methodology enabled the classification of spatial risk levels for the occurrence of elevated Cd levels in the Department of Lima. Areas with the highest overlap of high-risk categories were identified as zones of most significant concern.

3. Results

3.1. Statistical Distribution of Cd and Edaphic Properties: Identification of Anomalies and Geochemical Heterogeneity

Table 1 and Table 2 present the results of the physicochemical, metal, and metalloid analyses of the 410 soil samples examined. Total Cd content in the soils ranged from 0.00 to 11.98 mg·kg⁻¹, with a mean of 0.54 ± 1.00 mg·kg⁻¹, a variance of 1.01, and a coefficient of variation (CV) of 185.09%. These findings indicate a high degree of spatial heterogeneity in Cd accumulation in soils from avocado plantations in Peru.

The statistical behavior of Cd in the assessed soils exhibited extreme positive skewness (skewness = 7.30), indicating a strongly right-skewed distribution. This distribution suggests that most of the sampling units contain low Cd levels, while a small subset displays anomalously high concentrations.

Regarding the shape of the distribution, the high kurtosis value (kurtosis = 73.46) classifies it as extremely leptokurtic. This statistical behavior indicates not only a higher probability of occurrence of extreme values compared to a normal distribution (kurtosis = 3), but also suggests that such extremes are not the result of random noise or stochastic fluctuations. Instead, they reflect recurrent and systematic patterns in the data.

In addition to Cd, several soil fertility-related variables, such as electrical conductivity (EC), available phosphorus (Olsen P), total nitrogen, exchangeable calcium (Ca²⁺), exchangeable sodium (Na⁺), effective cation exchange capacity (CECe), as well as divalent and monovalent cation ratios (Ca²⁺ + Mg²⁺/K⁺, Ca²⁺/Mg²⁺, and Mg²⁺/K⁺), also exhibited statistically anomalous distributions. These variables exhibited extreme values of positive skewness, leptokurtic distributions, and coefficients of variation. The total metals and metalloids in the soil reported in Table 2, including Mo, Ni, Ag, Pb, Se, Tl, Ca, Na, Zn, Sb, and Cr, displayed distributional characteristics comparable to those observed for Cd.

3.2. Spearman Non-Parametric Correlation Between Total Cd Content and Fertility Parameters and Total Metal Accumulation in Avocado Plantation Soils of Peru

Spearman correlation analysis was conducted to assess non-parametric monotonic relationships between cadmium (Cd) concentration and 26 soil physicochemical properties (Figure 1). The results revealed potential patterns of co-variance between total Cd accumulation and edaphic factors associated with phenomena such as cation exchange and soil sodicity.

Moderate positive correlations (R > 0.30) were observed between Cd and Na⁺, ESP, and Mg²⁺ (0.42, 0.39, and 0.31, respectively), along with low but statistically significant positive correlations with silt content and EMP (0.23 and 0.14, respectively). In contrast, moderate negative correlations were found with ECP (-0.32), and low negative correlations with the Ca²⁺/Mg²⁺ ratio, total cation exchange capacity (CECt), Mg²⁺/K⁺ ratio, sand content, and pH (-0.25, -0.23, -0.18, -0.16, and -0.16, respectively).

Figure 2 presents the Spearman correlation analysis between total Cd content and 22 total metals and metalloids in the soil. The results indicate high positive correlations between Cd and Zn, Mg, Pb, Na, and Cu (0.60, 0.50, 0.46, 0.44, and 0.44, respectively). Moderate positive correlations were also observed with Se, Sr, V, Cr, As, and Ag, with correlation coefficients ranging from 0.32 to 0.37. Additionally, low but statistically significant positive correlations were found with TI, Mn, Ni, Mo, Sb, Fe, and K, with values ranging from 0.13 to 0.28.

These associations, assessed using the non-parametric Spearman approach, reflect a systemic tendency for certain soil fertility and geochemical variables to co-vary with Cd accumulation. The results provide strong evidence of a joint multivariate behavior, in which processes such as cation exchange, sodicity, and geochemical co-accumulation positively regulate Cd concentrations. Conversely, in soils with higher pH and greater sand content, Cd levels appear to be negatively regulated.

3.2. Principal Component Analysis of Fertility Parameters and Total Metals Associated with Cd Accumulation in Avocado Plantation Soils of Peru

Principal component analysis (PCA) was conducted by selecting variables with the highest correlations with total soil Cd, thereby reducing collinearity and transforming the variables to homogenize variances. This approach accounted for 46.5% of the variance across components 1 to 4 (Figure 3). Component 4 was included in the analysis because Cd contributed more significantly to it than to components 2 and 3. The variables contributing most to the variance explained by component 1 were Zn, sand content, silt content, Cd, Cu, exchangeable sodium percentage (ESP), and available phosphorus. In contrast, the variables most strongly associated with component 4 were Mg, sand content, silt content, and Cd. These findings highlight the significant influence of primary mineral weathering processes in agricultural soils under avocado cultivation, resulting in the simultaneous release of metals and metalloids, including Zn, Cu, Mg, and Cd. Furthermore, the analysis underscores the importance of variables that modulate Cd accumulation, particularly sand and silt content, ESP, and pH.

The clustering analysis using ellipses, with the department of origin as a categorical variable, reveals notable patterns. The department of Lima (pink ellipse) shows a pronounced orientation to the right (Dim1) of the graph, influenced by high values of Cd, Zn, Mg, and ESP. In contrast, the departments of Huancavelica (light blue ellipse), Ayacucho (green ellipse), and Ancash (red ellipse) display considerable overlap, with centroids positioned closely together, indicating similar soil profiles among these regions. These three departments exhibit an inclination toward variables negatively associated with component 1, such as pH, sand content, and ECP. However, the presence of both high and low Cd content hotspots cannot be ruled out, as the proximity of their centroids to the origin of the biplot suggests a multivariate profile not determined by compositional extremes, but by an intermediate balance of variables. In the case of Ica (blue ellipse), the ellipse is oval-shaped and parallel to Dim1, indicating low internal dispersion, but with multivariate similarity to the department of Lima, particularly in terms of high contents of Cd, Zn, Mg, and ESP, with silt content also contributing to this similarity. Lastly, the department of La Libertad (purple ellipse) is represented by a compact, circular ellipse of small size, signifying low internal dispersion and high compositional homogeneity among the analyzed samples. Its centroid lies in the negative quadrant of Dim1, opposite to those of Lima and Ica, suggesting soils characterized by low total Cd concentrations and low ESP levels, alongside higher pH values and greater sand content.

3.3. Application of the Geodetector Model for Identifying Determinant Factors in the Spatial Variability of Total Soil Cd in Avocado Plantations of Peru

Table 3 presents the results of the Geodetector analysis, highlighting Zn, Mg, life zone (LZ), soil parent material (SPM), climate, and geological age (GA) as the factors with the highest explanatory power on the spatial variation of total soil Cd. These factors accounted for 27.88%, 26.20%, 25.54%, 23.38%, 20.73%, and 18.21% of the variation, respectively, all with high statistical significance (p < 0.01). However, no statistically significant differences were observed among these variables (Figure 4), suggesting that the spatial variability of total soil Cd is likely driven by the combined effect of these variables, which implies the involvement of interrelated geological, geochemical, and climatic processes.

Additionally, the results indicate that Pb, soil type (ST), geomorphology (GT), exchangeable sodium percentage (ESP), and exchangeable calcium percentage (ECP) also contribute to the spatial variation of soil Cd, with values of 15.53%, 14.81%, 11.60%, 10.28%, and 8.59%, respectively, all statistically significant (p < 0.01). Similar to the most influential variables, these soil characteristics do not exhibit statistically significant differences among themselves (Figure 4), suggesting a combined effect on the spatial variation of Cd in soils under avocado cultivation. It is necessary to specify that these variables integrate characteristics related to the cation balance between Na and Ca, as well as their interaction with Cd accumulation.

Other significant factors (p < 0.01), though with lower explanatory power, included silt content (5.99%), sand content (5.46%), and available phosphorus (3.74%). However, the individual influence of these variables on Cd accumulation is statistically equal. In contrast, variables such as clay content, pH, and soil organic matter, although traditionally considered influential, showed a low q-value between 2.43 and 2.45%, with no significant differences (p > 0.05).

3.4. Non-Parametric Comparison of Total Cd Levels According to Parent Material of Origin, Taxonomic Type (WRB) of the Soil and the Geological Age of the Parent Material, in Soils Cultivated with Avocado

The Kruskal–Wallis test revealed significant overall differences in Cd levels (mg·kg⁻¹) among the different soil types (Figure 5). Post hoc analysis using Dunn’s test showed that Leptosols had significantly higher total soil Cd levels compared to Cambisols (medians: 0.64 vs. 0.20; z = 5.17; p-adj < 0.0001), Calcisols (0.64 vs. 0.20; z = 3.91; p-adj < 0.0001), and Luvisols (0.64 vs. 0.18; z = -3.35; p-adj = 0.03). These differences suggest that Leptosols, characterized by limited pedogenic evolution and potentially low contents of clay and stabilized organic matter, may promote Cd accumulation, in contrast to more developed soils or soils with higher immobilization capacity. Although other soil types, such as Phaeozems and Vertisols, also exhibited high medians, the small sample size in these cases reduced statistical significance, preventing the detection of real differences after conservative Bonferroni adjustment in Dunn's test. These findings highlight the significance of soil type as a crucial modulator of Cd geochemical behavior in avocado agroecosystems, indicating a strong interaction with factors related to parent material and soil mineralogy.

Highly statistically significant differences in total soil Cd concentrations were observed based on the origin of the parent material (Figure 6). Soils developed from residual parent material derived from intermediate igneous rocks exhibited the highest Cd levels (1.33 mg·kg⁻¹), significantly exceeding those from clastic sedimentary rocks (0.26 mg·kg⁻¹; z = -5.41; p-adj < 0.0001) and felsic igneous rocks (0.20 mg·kg⁻¹; z = 5.36; p-adj < 0.0001). Similarly, soils formed from residual parent material derived from intermediate igneous rocks also showed significantly higher values of accumulated Cd than soils with transported parent material of alluvial (0.24 mg·kg⁻¹; z = -6.43; p-adj < 0.0001) and aeolian (0.18 mg·kg⁻¹; z = -5.90; p-adj < 0.0001) origin.

Similarly, soils developed from glacially transported parent material exhibited a high median (1.12 mg·kg⁻¹); however, due to limited sample size and resulting low statistical power, no significant differences were detected compared to most other groups. In contrast, soils derived from aeolian-transported parent material showed the lowest Cd median (0.18 mg·kg⁻¹), with no significant difference to most other groups (p-adj > 0.05). This behavior suggests that soils formed from intermediate igneous rocks, potentially enriched in primary amphibole group minerals such as hornblende, commonly found in intermediate igneous rocks like diorites and tonalites, may represent a significant geogenic source contributing to Cd accumulation in soils.

Regarding the geological age of the parent material, highly significant differences in total Cd concentrations were observed in soils cultivated with avocado (Figure 7). Particularly, soils developed from Cretaceous geological formations exhibited the highest Cd concentrations (1.21 mg·kg⁻¹), significantly exceeding those associated with most other geological periods. The most pronounced differences were found in comparison with soils derived from Quaternary (0.24 mg·kg⁻¹; z = -6.13; p-adj < 0.0001), Permian (0.20 mg·kg⁻¹; z = -5.10; p-adj < 0.0001), and Paleogene (0.12 mg·kg⁻¹; z = -4.29; p-adj < 0.0001) formations, all with strong statistical significance. Additionally, though more moderate, significant differences were observed with soils of Jurassic (0.26 mg·kg⁻¹; z = -3.14; p-adj = 0.05) and Neogene (0.28 mg·kg⁻¹; z = -3.20; p-adj = 0.04) origin.

On the other hand, although soils developed on Devonian, Ordovician, and Paleogene formations exhibited low Cd medians (0.30, 0.32, and 0.12 mg·kg⁻¹, respectively), they did not significantly differ from those developed on Cretaceous formations (p-adj > 0.05). This lack of statistical significance may be attributed to limited statistical power, resulting from both small sample sizes and high intra-group variability.

3.5. Spatial Representation of Influential Factors in Cadmium Distribution Based on Geodetector Analysis

Analysis using the Geodetector (GD) algorithm identified a statistically significant relationship between the spatial distribution of cadmium (Cd) in soils and six categorical factors. By visually representing these factors: parent material (PM), geological age (GA), geomorphological type (GT), soil type (TS), climate (CLI), and life zone (LZ), it was possible to identify regions where specific classes exhibited higher average Cd concentrations compared to other categories within the same variable (Figure 8).

In terms of parent material, soils developed on intermediate igneous rocks (r-ii) exhibited elevated cadmium levels. A similar spatial pattern was observed for certain Cretaceous geological formations, which overlapped with areas of higher Cd concentration. Distinct associations were also found with specific soil types, particularly Leptosols, and ecological zones such as tropical humid montane forest (th-MF), tropical low montane hawthorn steppe (tl-MHS), and low montane tropical desert scrubland (lm-TDS). Additionally, the climatic category "semi-arid with dry-temperate winter" (D(i)B′) and certain geoform types, namely alluvial plains (P-aluv), mountains on sedimentary rock (M-sr), mountains on volcanic rock (M-vr), and undifferentiated terraces (T-u), were spatially associated with areas of increased Cd accumulation.

This visual interpretation is supported by the quantitative results of the q-statistic, which measures the extent to which each factor contributes to the spatial variability of Cd. The variables with the highest explanatory power were life zone (LZ) and parent material (PM), with q-values of 0.26 and 0.23, respectively. These were followed by climate (CLI) with q = 0.21 and geological age (GA) with q = 0.18. In practical terms, this indicates that approximately a quarter of the observed variation in soil Cd levels can be attributed independently to each of these factors. All q-values were statistically significant, providing robust support for their interpretation.

Figure 9 presents the potential risk map for the occurrence of elevated Cd accumulation values in soils across the agricultural areas of the Lima department, encompassing a total of 224,349.27 hectares. This map reflects the integrated effect of the most relevant categories from the key factors influencing the spatial variability of Cd. Four potential risk zones were identified and classified into levels: low, moderate, high, and very high. The majority of the agricultural area falls within the low-risk category (Category 1), covering approximately 175,106.35 hectares, which accounts for approximately 78.1% of the total area. This is followed by Category 2 (moderate risk), with 28,446.39 hectares (12.7%), Category 3 (high risk), with 14,262.29 hectares (6.4%), and Category 4 (very high risk), comprising 6,534.24 hectares (2.9%).

4. Discussion

4.1. Geochemical Implications of Anomalous Cd Distributions and Their Association with Edaphic Variables in Avocado Plantations of Peru

The pronounced positive skewness and extreme leptokurtosis of total Cd content in the soil (median = 0.26 ± 1.00) of avocado plantations in Peru, cultivated in soils with a pH of 7.42 ± 0.43, indicate a high probability of elevated Cd outliers' occurrence. These findings highlight the need for targeted attention in affected areas, as even small portions of the cultivated land may exceed the permissible Cd limits for agricultural soils, which range from 1 to 3 mg·kg⁻¹ [28].

This response is highly relevant for diagnosing areas with elevated total soil Cd, as its occurrence is driven by systemic factors rather than random variation. The statistical evidence supports the existence of edaphic hotspots, induced by localized micro-environmental and geochemical processes, such as selective adsorption onto specific soil fractions or the influence of point sources of contamination [29].

Univariate statistical analysis also revealed that soil properties related to cation ratios (Na⁺, Ca²⁺, CEC, Ca²⁺+Mg²⁺/K⁺, Ca²⁺/Mg²⁺, and Mg²⁺/K⁺), as well as the total accumulation of metals and metalloids (Zn, Pb, Na, Se, Cr, Ag, Tl, Ni, Mo, Sb, and Ca), exhibit heterogeneous behaviour similar to that of Cd [30,31]. These findings suggest that these soil characteristics may be directly or indirectly modulating the spatial distribution of Cd. Consequently, the base cation saturation ratio and the geogenic co-accumulation of metals appear to play a significant role in shaping Cd dynamics in avocado plantation soils of Peru.

4.2. Edaphic Mechanisms of Cd Accumulation: Insights from Multivariate Statistical Analysis

Spearman’s correlation analysis and principal component analysis revealed positive correlations and a strong influence on the explained variance between total Cd and variables related to soil sodicity (Na⁺ and ESP), suggesting that soils with high sodium saturation tend to accumulate more Cd. This finding is consistent with Zahedifar [32], who reported that sodic (0.5–0.96 mg·kg⁻¹) and saline-sodic (2.8–3.6 mg·kg⁻¹) soils exhibit higher total Cd accumulation compared to saline (0.2–0.44 mg·kg⁻¹) and normal soils (0.4 mg·kg⁻¹). In these soils, characterized by high exchangeable Na⁺ content, the phenomenon of colloidal dispersion of the clay–humic complex and soil aggregates predominates [33]. This process increases the exposure of negatively charged sites and reduces cation selectivity at exchange sites, thereby promoting the adsorption of metal cations such as Cd [34].

Among the variables most positively correlated with Cd were Mg²⁺ and EMP, while the Ca²⁺/Mg²⁺ ratio and ECP exhibited significant negative correlations. These results suggest that cationic competition between Na⁺ and Mg²⁺ against exchangeable Ca²⁺ may enhance Cd fixation in soil matrix constituents and increase its mobility into the soil solution [35]. Principal component analysis supports these findings because it identifies the high influence of total Mg on the explained variance of soils with high Cd accumulation. It also reinforces the hypothesis of cationic competition between Ca²⁺ and Cd, given their similar elemental properties, such as ionic radius (Cd²⁺ = 0.95 Å; Ca²⁺ = 0.99 Å) and divalent charge [36,37].

Low but statistically significant negative correlations were also observed between Cd and parameters such as pH, total CEC, and sand content, along with a low positive correlation with silt content. Regarding pH, it is well established that this parameter influences the solubility and mobility of labile Cd in soils [38] and is considered one of the most important predictor variables for Cd bioaccumulation in crops [39]. Specifically, Cd becomes more soluble under acidic soil conditions and tends to precipitate as cadmium carbonate or hydroxide under alkaline soil conditions [40]. However, in the present study, total soil Cd accumulation was determined through chemical digestion followed by analysis using inductively coupled plasma mass spectrometry (ICP-MS). This technique quantifies Cd across all soil phases, including Cd adsorbed on ion exchange complexes, bound to Fe, Mn, and Al oxides, associated with soil organic matter, incorporated into primary and secondary minerals, and precipitated as salts or insoluble complexes [41]. Consequently, this methodology does not measure soluble or labile Cd, nor does it distinguish the chemical forms of Cd speciation. Therefore, the correlation between pH and total Cd is weak; however, this does not preclude the possibility that soils with parent material derived from alkaline rocks may have lower Cd content.

Cadmium accumulation has also been shown to depend on soil particle size distribution [28]. In our principal component analysis, the sand and silt content emerged as prominent variables, explaining a significant proportion of the influence of component 1 on the explained variance of these soils. Our results indicate that soils with a higher silt content tend to accumulate more Cd than those with a higher sand content. Although previous studies have reported that soils with a greater proportion of fine particles (< 2 µm) adsorb 1.2 to 1.4 times more Cd than those with a higher proportion of coarse particles [28]. Our findings in avocado plantation soils of Peru reveal that silt content (particles between 2 and 20 µm) has a greater influence on total Cd accumulation than clay content. This may be attributed to silt functioning as a primary Cd reservoir, potentially through the formation of surface complexes with carbonates originating from the weathering of plagioclase [42]. Conversely, the observed negative correlation between sand content and Cd accumulation can be explained by the enhanced vertical migration of Cd in sandy soils. Higher sand content is typically associated with greater soil aeration and more alkaline conditions, both of which increase Cd mobility, facilitating its leaching into deeper soil layers [43].

Spearman correlation analysis between total Cd concentrations and the total contents of 22 metals and metalloids in the soil revealed systematic patterns of geochemical co-accumulation linked to the parent material. The strongest positive correlations were observed with Zn, Mg, Pb, Na, and Cu. This association can be attributed to the geochemical inheritance of the parent material [44]. In our study, intermediate igneous rocks, specifically Cretaceous diorites and tonalites, were identified as the parent material associated with the highest levels of Cd accumulation. These rocks often coexist with primary minerals, such as plagioclase and biotite, within their structure.

Cadmium (Cd), zinc (Zn), and lead (Pb) share a combination of structural, electronic, and thermodynamic similarities that influence their mobility, adsorption, and partitioning among mineral and organic soil phases [28]. These elements exhibit comparable ionic radii in their most common oxidation states (0.95 Å for Cd²⁺, 0.74 Å for Zn²⁺, and 1.19 Å for Pb²⁺) [45]. This similarity supports their participation in isomorphic substitution processes within crystalline structures, including carbonates, oxides, clay minerals, and primary sulfide-type minerals such as sphalerite (ZnS), galena (PbS), and greenockite (CdS) [28].

Zinc (Zn) showed the strongest correlation with cadmium (R = 0.60) among the 48 soil variables evaluated and was also the variable that contributed most to the explained variance in the soils studied. Both Zn and Cd possess a closed d¹⁰ electronic configuration [42], indicating that their d orbitals are fully occupied. As a result, they exhibit low redox reactivity, that is, they are not readily oxidized or reduced. Instead, these elements tend to form electrostatic complexes by binding to negatively charged functional groups, such as carboxylic (-COO⁻) and phenolic (-OH) groups, present in soil organic matter or on the surfaces of clay minerals [2]. Lead (Pb), on the other hand, has a 6s² electronic configuration, known as an “inert pair,” which refers to an energetically stable electron pair that is not easily ionized or shared [46]. This property makes Pb a chemically heavier element, and like Zn and Cd, it forms electrostatic bonds with oxygenated functional groups rather than strong covalent bonds [44].

The physicochemical characteristics of Cd, Zn, and Pb support their geochemical co-adsorption onto soil colloids, particularly in soils with a high content of humic complexes and/or the presence of Fe and Mn oxides, typical of highly weathered acid soils [42]. However, despite the high variability in clay content (18.56 ± 9.78%) and organic matter content (2.04 ± 1.69%) in the studied soils, no significant correlations with Cd were observed. These findings underscore the relevance of alternative mechanisms, such as the influence of weakly weathered parent material near the soil surface, in controlling Cd accumulation. This characteristic is typical of young soils, such as Leptosols, which are derived from parent materials rich in primary minerals, including plagioclase and biotite, and are contaminated with trace metals or sulfide minerals that have adsorbed metal ions.

4.3. Soil and Geological Factors Controlling Total Cd Accumulation in Avocado-Cultivated Soils: Influence of Taxonomic Type, Parent Material, and Geological Age

Non-parametric tests were conducted to identify specific edaphological characteristics associated with elevated total Cd content in soils. The results indicate that Leptosols accumulate significantly more Cd than Cambisols, Calcisols, and Luvisols. Likewise, parent materials of Cretaceous origin derived from intermediate igneous residual material reached significantly higher Cd levels. To enable a more comprehensive soil analysis of these interactions, it is essential to clarify the key differences among these soil types.

Leptosols are poorly developed soils, typically formed over cemented parent material and characterized by a shallow effective depth (< 50 cm), which can significantly constrain agricultural activity [47]. They commonly exhibit lithic, petrocalcic, or gypsic contact, which causes edaphic restriction in root growth and limited horizonation of the soil profile [48]. Additionally, Leptosols have low water-storage capacity and good drainage, features typical of soils with low organic matter content in the surface horizon [49]. However, due to their shallow profile, metals such as Cd tend to accumulate in the upper soil layers, particularly in soils with a basic pH, where Cd becomes less soluble [50].

The more recent geological age of the parent material in Leptosols indicates a higher proportion of primary minerals in the early stages of weathering [51]. The low degree of alteration of the primary minerals, mainly amphiboles and micas, supports the higher content of trace Cd in their crystalline structure [52]. Thus, in highly technified agricultural soils, the frequent use of acidified irrigation water leads to the dissolution of these minerals and the release of Cd into the solid phase and the ion exchange complex of the soil [53].

Cambisols are also poorly evolved soils, with incipient alteration of the parent material [47]. However, this alteration is sufficient to modify the soil’s physical and chemical properties compared to the original material. These soils are typically characterized by the presence of a cambic subsurface diagnostic horizon [47]. The cambic horizon exhibits colour and structure formation derived from the parent material, without significant illuviation processes; that is, its differentiation is due to the alteration of primary minerals or the incipient formation of secondary minerals.

Thus, in contrast to Leptosols, Cambisols contain higher levels of clay minerals and soil organic matter, reflecting greater pedogenetic maturity [54]. For this reason, it tends to accumulate Cd through organic chelation complexes or metal-specific adsorption mechanisms [54]. However, the higher soil profile horizon of Cambisols likely produces a dilution effect on total Cd concentrations, whereas Leptosols exhibit higher concentrations (>3 mg·kg⁻¹) within a reduced topsoil layer thickness. Nonetheless, Cambisols with high levels of Cd content (>1 mg·kg⁻¹) have been identified in avocado plantations, particularly in alluvial and glacially transported soils of Quaternary geological age. These edaphic hotspots within certain Cambisols are driven by dynamic geomorphological features such as undifferentiated terraces, mountainous areas formed on volcanic rocks, and riverbeds.

Luvisols, unlike Leptosols and Cambisols, exhibit a higher degree of pedogenic development and are characterized by the presence of an argic subsurface horizon (Bt) with significant clay accumulation, typically exceeding 15 cm in thickness [47]. These soils have a fine texture, moderate organic matter content in the A horizon, high CEC (> 20 cmol·kg⁻¹), neutral to alkaline pH, and are dominated by illite-type minerals [55]. In these Luvisols, the highest concentrations of Cd (> 1 mg·kg⁻¹) have been observed in soils of residual origin derived from intermediate igneous rocks (diorite and tonalite) formed during the Cretaceous period. This characteristic coincides with the highest Cd values found in Leptosols (> 3 mg·kg⁻¹). However, it is likely that in Luvisols, Cd released through weathering is present in a higher concentration in the argic horizon, where it is firmly attached to the specific surface of clay minerals, thereby diluting its concentration in the arable layer of the soil [56].

Finally, significant differences were also observed between Leptosols and Calcisols. Calcisols are defined by the presence of a calcic subsurface horizon (Bk), characterized by a carbonate content exceeding 15% within this layer [47]. These soils typically exhibit a neutral to alkaline pH, a moderate CEC lower than that of Luvisols, and a low organic matter content, primarily due to the leaching of soluble organic compounds [57]. Notably, Calcisols presented the lowest levels of Cd accumulation in the arable layer of avocado plantation soils (median = 0.2 mg·kg⁻¹), with pronounced differences compared to Leptosols. The lowest Cd concentrations in these soils were associated with parent material of aeolian-transported origin, which contributes to the formation of sandy layer deposits, influencing the generation of sandier soils with good infiltration.

It is likely that, in the arable soil layer, most cation exchange sites are saturated with calcium, thereby displacing Cd into the soluble phase of the soil [37,58]. Then, in the soil solution, Cd forms complexes with short-chain soluble organic acids [59], which can migrate with the downward flow of water, particularly in permeable horizons (sandy or coluvial). Upon reaching the calcic subsurface horizon, which is rich in carbonates, Cd can precipitate as otavite (CdCO₃) or become incorporated into mixed Ca-Cd carbonate phases [43,60,61].

5. Conclusions

The findings of this study demonstrate that total cadmium (Cd) content in avocado-cultivated soils of Peru exhibits a highly heterogeneous spatial distribution, with notably elevated values (>3 mg·kg-1) that are not attributable to random variability, but rather to the systematic influence of edaphic and geogenic factors. The results provide strong evidence that the main explanatory variables of Cd accumulation in basic pH are of pedo-geochemical nature. Key explanatory variables include soil taxonomic classification, parent material, geological age, total content of Zn, Pb, and Mg, exchangeable sodium percentage (ESP), and specific cation ratios (Ca²⁺/Mg²⁺). Moreover, this research identifies Leptosols formed on Cretaceous intermediate igneous parent materials (diorite and tonalite) as having the highest total Cd concentrations, indicating a pronounced lithogenic contribution. This is likely linked to the weathering of trace-metal-rich primary minerals and the shallow depth of the arable soil layer. Based on the integrated results, a risk map was developed to illustrate the potential occurrence of elevated total Cd levels in soils across the Lima region. The analysis indicates that approximately 49,242.92 hectares, representing 22% of Lima’s agricultural land, have a high probability of occurrence and are at moderate to very high risk of Cd accumulation in soils.

The identification of soil types with lower cadmium accumulation, such as Calcisols, Cambisols, and Luvisols, compared to Leptosols, provides critical insights for developing soil management strategies that aim to mitigate Cd levels in vulnerable agricultural areas. These more pedogenetically evolved soils, characterized by better profile development and the presence of subsurface horizons with chemical immobilization capacity (argic or calcic horizons), offer natural attenuation mechanisms for Cd. Such mechanisms include organic chelation, specific adsorption onto high cation exchange capacity (CEC) clays, and precipitation into insoluble forms, such as mixed carbonates. In contrast, Leptosols, with their limited soil profile differentiation and close proximity to lithogenic parent material, tend to favor surface Cd accumulation. Consequently, Calcisols, Cambisols, and Luvisols may serve as functional reference soils for guiding the sustainable management of avocado-growing areas in Peru.