The study was conducted at the San Ramón Agrarian Experimental Station of the National Institute of Agrarian Innovation (INIA), located in the Yurimaguas district, Alto Amazonas province, Loreto region, Peru (approximately 5.90° S, 76.12° W). Yurimaguas has a humid tropical climate typical of the lower Amazon, with average monthly temperatures of approximately 26–27 °C (temperature normals reported for Yurimaguas) and annual precipitation ranging from ~2,500 to 2,600 mm, characterized by pronounced seasonal rainfall. Intra-annual temperature variation is minimal, and rainfall is persistent throughout the year. Winds are generally weak (around 1–2 m·s^-¹), and relative humidity remains high year-round (> 70%). These climatic conditions have been documented for Yurimaguas and the Huallaga basin in both international and regional climate sources, including compilations of local monthly climate normals (12).

Cocoa (Theobroma cacao L.) seedlings of the CCN-51 clone were used, established in a greenhouse, and evaluated up to four months of age. The seedlings were grown in 2 kg polyethylene bags filled with a local agricultural substrate (one plant per bag). The soil used in the trial was pre-incubated for 45 days after amendment-conditioning with phosphorus sources (phosphate rock) and prior to transplanting to stabilize the soil–amendment reactions. Foliar application of zinc sulphate began once the seedlings developed true leaves, with four applications administered at 15-day intervals using a manual backpack sprayer with a fine-cone nozzle to ensure uniform foliage coverage. Biometric evaluations (plant height, stem diameter, number of leaves, root length, and root volume) were conducted every 30 days during the four-month greenhouse period.

The preliminary soil characteristics prior to treatment application were as follows: 5.29 pH, 1.61% organic matter content (OM), 0.07% total nitrogen (N), 6.95 mg·kg^-¹ available phosphorus (P); 75.97 mg·kg^-¹ available potassium (K), 0.06 mg·kg^-¹ available cadmium (Cd), 0.25 mg·kg^-¹ total cadmium, and a sandy loam texture. A 3 × 4 factorial arrangement was implemented in a completely randomized design, with three replicates per factor combination, resulting in 36 experimental units. The results of each treatment are seen in Table 1. The selection of rock phosphate (RFP) and zinc sulfate (ZnSO₄) levels was based on agronomic, edaphic, and mechanistic considerations. The RFP rates (0, 57.27, and 114.55 mg·kg^-¹) represent a gradient from zero phosphorus input to a high application dose, designed to evaluate the potential of RFP to modify soil pH, increase phosphorus availability, and promote Cd immobilization through the formation of Cd–phosphate complexes, while also accounting for the possible presence of Cd impurities in natural phosphate sources. The ZnSO₄ doses (0, 98.80, 197.50, and 527.80 mg·plant^-¹, applied foliarly) were defined based on reported ranges effective in inducing Zn–Cd physiological antagonism without causing phytotoxic effects. This approach enabled the assessment of dose-dependent responses and interactions with soil phosphorus. Consequently, the factorial experimental design facilitates the identification of management combinations with the greatest potential to mitigate Cd transfer from soil to plant.

The operational database was structured into three measurement compartments: soil, leaf tissue, and seedling, each containing representative variables. Soil variables included pH, organic matter (OM, %), total nitrogen (N, %), available phosphorus (P₂O₅, mg·kg^-¹), available potassium (K₂O, mg·kg^-¹), exchangeable cations (Ca²⁺, Mg²⁺, K⁺, Na⁺), effective cation exchange capacity (ECEC), and trace metals (available Cd²⁺, Zn²⁺, and Cu²⁺). Leaf tissue variables comprised macronutrients (P₂O₅, K₂O, Ca²⁺, Mg²⁺, Na⁺) and micronutrients (Fe³⁺, Mn²⁺, Cu²⁺, Pb²⁺, Cd²⁺, Zn²⁺), expressed in ppm or percentage depending on analytical practice. Seedling variables included substrate pH, plant height (cm), stem diameter (mm), number of leaves, root length (cm), and root volume (mL). Analytical determinations followed INIA’s standard laboratory protocols: pH was measured by potentiometry, macronutrients by extraction and colorimetry or flame photometry depending on the measured analyte, and metals by acid digestion with detection by Atomic Absorption Spectroscopy (AAS) or Inductively Coupled Plasma–Optical Emission Spectrometry (ICP–OES). Treatment coding and sample traceability were maintained following the field logbook and experimental data matrix.

All analyses were conducted in R version 4.x using RStudio/Posit IDE, with fully reproducible scripts in R Markdown or Quarto. The analytical workflow relied on the tidyverse, ggplot2, rstatix, corrplot, FactoMineR/factoextra, MASS, and car packages (13–15).

Descriptive statistics and exploratory data analysis (EDA) included means, standard deviations (SD), 95% confidence intervals (95%CI), coefficients of variation (CV), skewness, and kurtosis to characterize soil, leaf, and seedling variables and verify assumptions before modelling (15, 16). Correlation analyses and heatmaps were used to examine Pearson correlation matrices and visualise covariation (e.g., P₂O₅–Cd²⁺, K₂O–Zn²⁺) and detect collinearity before multivariate analysis (17). This study reported high correlations within each compartment. Principal Component Analysis (PCA) was applied separately to soil, leaf, and seedling datasets using prcomp and FactoMineR to reduce dimensionality, synthesise latent axes (e.g., P₂O₅–(Cd²⁺, Zn²⁺)–K₂O vs pH/Na⁺), and visualise gradients (18, 19). The study details the variance explained by the first two components. Factorial ANOVA/GLM (ZnSO₄ × RFP) was used to assess main and interaction effects, with verification for assumptions and effect sizes (20). A pattern of non-parallel lines in the interaction plots supports a ZnSO₄ × RFP interaction (Table 1). Visual moderation analysis (scatterplots and simple slopes), using dispersion plots by RFP level and simple slope analysis, helps to interpret the direction and magnitude of the interaction (21), where slopes are expected to decrease with higher RFP levels.

Boxplots were used to compare medians, interquartile ranges (IQRs), and outliers robustly among RFP × ZnSO₄ treatments (22), showing clear separations between treatment combinations. Linear Discriminant Analysis (LDA) was applied as a supervised classification method to evaluate treatment separability based on seedling traits, complementing PCA (18). Finally, the relative uptake index (foliar Cd²⁺/soil Cd²⁺) was used as a dimensionless metric to prioritize treatments with lower soil-to-plant transfer efficiency, providing evidence of attenuation under high RFP levels (6, 23).

Prior to multivariate analyses (PCA and LDA), all variables were centered and scaled using z-score standardization, resulting in a mean of zero and a unit variance for each variable. This transformation was required to account for differences in units and magnitudes among variables (e.g., pH, mg·kg^-¹, cmol(+)·kg^-¹, and biometric traits), ensuring that no variable artificially dominated the multivariate structure. No additional logarithmic transformations were applied, as preliminary exploratory analyses did not indicate extreme deviations or skewness that would justify their use.

Overall, the soil exhibited low to moderate acidity (pH: M = 5.52, SD = 0.61, 95% CI [5.32, 5.72]; n = 36), with a right-skewed distribution (skew = 1.01) and mesokurtic kurtosis of 0.25. The pronounced dispersion of Cd²⁺ and P indicates heterogeneous soil nutrient supply, which may influence foliar Cd²⁺ uptake and modulate the plant’s response to ZnSO₄ (Table 2). Foliar Zn²⁺ levels were elevated, with moderate dispersion (M = 78.49 mg·kg^-1, SD = 25.65, 95% CI [69.86, 87.11], CV = 32.7%). Regarding the Zn²⁺–Cd²⁺ ratio, the average foliar Zn²⁺: Cd²⁺ ratio was approximately 7.1:1. This ratio suggests a favourable relative availability of Zn²⁺ despite the high and heterogeneous foliar Cd²⁺ concentrations. Nonetheless, descriptive statistics alone do not demonstrate Zn–Cd antagonism; rather, they contextualize the magnitude and variability on which the treatments operate (Table 3). Seedling biometrics showed that the aerial traits were stable to moderately variable: seedling height (M = 22.76 cm, SD = 3.81, 95% CI [21.47, 24.05], CV = 16.8%) and stem diameter (M = 7.92 mm, SD = 0.79, 95% CI [7.66, 8.17], CV = 10.1%). This pattern indicates that that the greater dispersion observed in root traits suggests that management factors (ZnSO₄ and phosphate rock) and/or the availability of Cd²⁺ and P (PO₄³^-) in the soil exert a more substantial influence on root architecture than on aerial biomass at this stage.

The correlations based on soil analysis results revealed a strong covariation block among available phosphorus (P₂O₅; as PO₄³^- in solution), total cadmium (Cd²⁺), exchangeable cations, and pH. Extremely high correlations were observed between P₂O₅ and exchangeable K⁺ (cmol(+)·kg^-¹; (r(34) = 0.97, p < 0.001), P₂O₅ and exchangeable Na⁺ (cmol(+) ·kg^-¹; r(34) = 0.95, p < 0.001), as well as between P₂O₅ and Cd²⁺ (r(34) = 0.93, p < 0.001). Cadmium also correlated strongly with K⁺ (r(34) = 0.92, p < 0.001) and Na⁺ (r(34) = 0.90, p < 0.001). Soil pH was positively correlated with P₂O₅, Cd²⁺, K⁺, and Na⁺ (r ≈ 0.73–0.75, p < 0.001), and moderately with soil Zn²⁺ (r = 0.45, p < 0.01). In contrast, K₂O (mg·kg^-1) showed weak or non-significant relationships with most variables (|r| ≤ 0.27, ns). These correlations suggest that phosphate rock (RFP) applications simultaneously increased soil P₂O₅ (PO₄³^-), pH, and exchangeable cations (K⁺, Na⁺), while soil Cd²⁺ co-increased with these factors, likely due to cadmium impurities in the RFP and/or co-accumulation in the extracted fractions. The positive pH–Cd²⁺ correlation (rather than the expected negative relationship) supports the hypothesis of a shared source of Cd rather than changes in availability driven by soil chemistry (Figure 1).

Correlations based on leaf concentration data revealed strong covariation between nutrients and cadmium. Very high correlations were observed for P₂O₅–Cd²⁺ (r(34) = 0.88, p < 0.001) and Na⁺–Cd²⁺ (r(34) = 0.84, p < 0.001). High correlations were also found for K₂O–Zn²⁺ (r(34) = 0.72, p < 0.001), P₂O₅–Na⁺ (r(34) = 0.77, p < 0.001), and P₂O₅–Zn²⁺ (r(34) = 0.64, p < 0.001). The Zn²⁺–Cd²⁺ association was moderate to high (r(34) = 0.59, p < 0.001), and K₂O–Cd²⁺ showed a moderate correlation (r(34) = 0.37, p < 0.05). Only K₂O–Na⁺ was not significant (r = 0.12, ns). These results indicate that, in strictly correlational terms, foliar Zn²⁺ is positively associated with Cd²⁺ (r = 0.59). This pattern does not support a simple antagonistic relationship at the leaf level. Instead, it suggests treatment collinearity (e.g., higher ZnSO₄ doses applied alongside greater RFP inputs) or shared uptake drivers, such as increased plant vigour or co-transport under elevated P₂O₅ and Na⁺ supply (Figure 2).

In the biometric correlations of 4-month-old cocoa seedlings, strong and consistent associations were observed among growth traits. Seedling height and root length showed a very strong correlation (r(34) = 0.88, p < 0.001), indicating that taller seedlings also tend to develop longer roots. Seedling height and stem diameter were moderately to highly correlated (r(34) = 0.58, p < 0.001), and stem diameter and root length exhibited a similar moderate-high association (r(34) = 0.57, p < 0.001). Together, these three relationships suggest a coordinated growth pattern between aerial and root structures, supporting the notion that treatments enhance overall vigour (e.g., combinations of RFP and ZnSO₄) promote simultaneous aerial and root development (Figure 3).

Other associations were weak or non-significant: number of leaves with seedling height (r = 0.30, ns) and with root length (r = 0.22, ns). Root volume was also unrelated to the remaining traits (|r| ≤ 0.08, ns). This suggests that root volume was either more variable or less effectively captured by the experimental arrangement compared with root length and seedling diameter.

The PCA of the soil analysis results explained 85.4% of the total variance (Dim1 = 70.5%; Dim2 = 14.9%; Table 4; Figure 4). Dim1 grouped the fertility/chemical status variables, dominated by pH, P₂O₅, Zn²⁺, Cd²⁺, K⁺ (cmol(+)·kg^-¹), and Na⁺ (cmol(+)·kg^-¹), as shown in Table 5. The small angles between these vectors, all oriented toward the positive semi-axis of Dim1 in Figure 4, indicate strong positive correlations within this block, consistent with the patterns observed in the heat maps. Dim2 was clearly dominated by K₂O (the longest vector and the most significant contributor), showing that variability in soil K represents an independent component separate from the main pH–P₂O₅–(Zn²⁺, Cd²⁺)–exchangeable bases (K⁺, Na⁺) block (Figure 5).

The PCA of the foliar analyses conducted on 4-month-old cocoa seedlings accounted for 90.0% of the total variance across the first two dimensions (Dim1 = 67.6%; Dim2 = 22.4%; Table 6, Figure 6). Dim1 represents a clear ion-accumulation axis: K₂O and Zn²⁺ exhibited the most significant contributions (longest vectors and highest chromatic intensity), followed by Na⁺, whereas Cd²⁺ and P₂O₅ contributed to a lesser extent (Table 7; Figure 6). Crucially for the Zn²⁺–Cd²⁺ interaction hypothesis, the Zn²⁺ and Cd²⁺ vectors show opposite orientations along Dim2 (Zn²⁺ with a positive component and Cd²⁺ with a negative component), reflecting an inverse relationship along that axis. This indicates that, at equal levels of the general accumulation axis (Dim1), higher foliar Zn²⁺ levels are associated with lower foliar Cd²⁺ concentrations (Table 7). This pattern is consistent with an antagonistic effect of ZnSO₄ on Cd²⁺ uptake or translocation, aligning with the central premise of the study (Figure 6).

The PCA of the biometric analysis results for cocoa seedlings explained 71.6% of the total variance across the first two dimensions (Dim1 = 50.8%; Dim2 = 20.8%), as shown in Figure 4. The vector configuration indicates that Dim1 represents an aerial-size gradient: seedling height and root length load positively and with intermediate magnitude, while stem diameter and number of leaves also align with Dim1, though with smaller contributions (Table 8). The small angles between these vectors in Figure 5 are suggest the positive correlations among them. Dim2 is dominated by root volume, as indicated by the longest vector and highest contribution (Figure 5; Table 8), suggesting that root volume introduces orthogonal variation to aerial size. The biplot shows no distinct separation among the treatment groups (coloured points), suggesting that the RFP × ZnSO₄ combinations influence the magnitude of the traits rather than their covariation pattern. In the context of the study objective, improved morphological performance is associated with positive values along Dim1 (taller seedlings with longer roots) and Dim2 (greater root volume), which may serve as plausible indicators of enhanced tolerance to Cd²⁺ stress when Zn²⁺ supply is adequate.

A multigroup LDA was conducted using biometric variables as predictors and treatment as the grouping factor. The classification performance shown in the figure corresponds to an apparent accuracy of 81% (≈29 correct classifications out of 36 observations), which is markedly higher than the accuracy expected by chance (≈8.3% with 12 groups). An approximate 95% confidence interval for accuracy ranges from 0.68 to 0.93.

Regarding the separation structure, LD1 accounted for most of the separation among treatments (horizontal distribution, as shown in Figure 15), while LD2 provided secondary differentiation (vertical distribution). In practical terms, treatments with high LD1 values (right side of the plot; e.g., T11, T8, T9) exhibit a relatively stronger biometric profile, characterized by linear combinations reflecting greater height and diameter, and, to a lesser extent, root length, compared with treatments on the negative side of LD1. Treatments with negative LD1 values (left side; e.g., T2, T3) cluster seedlings with comparatively lower biometric values. LD2 adds further separation for specific subsets. For example, T1 and T12 are distinguished from others with similar LD1 values but different LD2 positions, suggesting contrasts along an orthogonal trait combination (e.g., root volume and leaf number versus height/diameter), which helps explain certain overlaps and classification errors. The central cluster, where several intermediate treatments (e.g., T4–T7) overlap, accounts for the misclassification fractions contributing to the overall 81% accuracy. From an experimental interpretation perspective, these intermediate treatments share similar biometric phenotypes, consistent with ZnSO₄ × RFP combinations that do not generate strong contrasts in seedling growth.

Within the framework of the research project “Antagonistic interaction between zinc sulphate and cadmium in cocoa seedlings (Theobroma cacao var. CCN-51) under phosphate rock application”, the observed separation indicates that biometric responses integrate the Zn²⁺–Cd²⁺ interaction signal as modulated by RFP. Treatments combining conditions that favour Zn²⁺↔Cd²⁺ antagonism tend to cluster at one extreme of LD1 (reflecting superior biometric performance), whereas less favourable combinations appear at the opposite extreme (reflecting poorer performance). The pattern of partial overlap suggests that, although the overall signal is strong, there is intragroup heterogeneity (e.g., variation among pots or individual seedlings) that should be more tightly controlled in future experimental designs.

In Figure 15 and Table 2, relative uptake (foliar Cd²⁺/soil Cd²⁺; lower values indicate better performance) varied widely among treatments (range = 33.08; maximum in T4 = 53.12; minimum in T9 = 20.04). The relative difference between the worst and best treatments was 62.3% = 1 − (20.04/53.12), indicating a substantial potential for reducing Cd²⁺ transfer to leaf tissue. Ranked from best to worst (lowest to highest relative uptake), the five treatments with the lowest values were: T9 (20.04), T2 (25.43), T1 (26.11), T12 (26.19), and T11 (26.84), as shown in Table 2. At the opposite end of the distribution, the treatments with the highest relative uptake were T4 (53.12), T6 (34.68), T3 (33.73), and T8 (33.27).

The colour of the bars encodes the dose of phosphate rock (RFP): 0 mg·kg^-1 (red), 57.27 mg·kg^-1 (green), and 114.55 mg·kg^-1 (blue). Based on the values in Table 4, the averages within each RFP level were as follows: RFP = 114.55 mg·kg^-1: ≈ 25.43 (T9, T10, T11, T12); RFP = 57.27 mg·kg^-1: ≈ 31.90 (T5, T6, T7, T8); RFP = 0 mg·kg^-1: ≈ 34.60 (T1, T2, T3, T4). Thus, the highest RFP dose (114.55 mg·kg^-1) was associated with an average reduction of approximately 26.5% in relative Cd uptake compared with the 0 mg·kg^-1 treatments (mean difference = 34.60 – 25.43 = 9.17 units). Moreover, intragroup variability was substantially greater at RFP = 0 mg·kg^-1 (range ≈ 27.69) than at RFP = 57.27 mg·kg^-1 (≈ 7.26) or RFP = 114.55 mg·kg^-1 (≈ 8.62), suggesting more stable control of Cd²⁺ transfer when the phosphorite dose is increased.

The ranking pattern supports the attenuation of foliar Cd²⁺ uptake under management combinations that incorporate higher RFP doses (blue bars in Figure 16), consistent with the antagonism hypothesis (i.e., reduced Cd²⁺ bioavailability through phosphate complexation and/or fixation, and/or synergistic interactions with ZnSO₄). In practical terms, T9 emerges as the most effective treatment (20.04), whereas T4 represents the least favourable scenario (53.12). Overall, the data presented in Figure 16 and Table 3 indicate that increasing RFP, together with selecting appropriate ZnSO₄ combinations (Figures 7–9), consistently reduces the foliar Cd²⁺/soil Cd²⁺ ratio while also decreasing variability. This trend is both statistically and agronomically significant for minimizing Cd²⁺ accumulation in the leaf tissue of cocoa seedlings.