The experimental areas considered in this study belong to the Sugarcane Breeding Program of the Federal University of São Carlos (UFSCar), one of the ten federal university members of the Interuniversity Network for the Development of the Sugarcane Industry in Brazil (RIDESA). Here, we considered two experiments, named Experimental Areas 1 and 2 (EA1 and EA2, respectively), corresponding to the first breeding stage, each located in strategic regions of the state of São Paulo, as detailed here. EA1 represents a traditional cultivation region with high soil fertility for sugarcane cultivation, while EA2 represents a region of crop expansion with low soil fertility. Other environmental conditions are available from the RIDESA/UFSCar coverage area, but they are usually considered only in the final assessment trials where a reduced number of genotypes are available. This allows us to account for genotype × environment (G × E) interactions; i.e., the genotypes are tested in multienvironmental trials (MET) across multiple crop-years and seasons. Before the experiments were planted, a high density of soil sampling was performed to determine the effect of the spatial variability. Fertilization and cultural treatments were carried out as recommended for sugarcane.

Before the implementation of the experiments, both areas were georeferenced. Soil sampling was carried out mainly in a regular grid, where samples were collected at equally spaced points and homogeneously distributed throughout the experimental region.

In EA1, 56 points (8.6 samples/ha) were sampled at depths of 0–20 cm. Each sample was composed of six subsamples collected from an average radius of 3 to 5 m around the central point (Supplementary Figure 1A). The same procedures occurred for EA2; however, since EA2 was larger than EA1, 68 sample points (7 samples/ha) were collected (Supplementary Figure 1B).

The soil samples were sent to the soil laboratory of the Department of Natural Resources and Environmental Protection at UFSCar for analysis and determination of chemical (macronutrients) and physical attributes of the soil.

Thirty days after the experiment was planted, apparent electrical conductivity (ECa) mapping was performed by using the EM38-MK2^® sensor (Geonics, Mississauga, Ontario, Canada), which operates with the principle of electromagnetic induction (EMI). In both experimental areas, the sensor, which was connected to a GPS receiver [type L1 (Trimble)], was used to record the geographic coordinates, and a data collector (Juniper Archer Field PC) was used to store information. Additionally, no vehicle was utilized to transport the sensor in either experimental area, i.e., it was manually operated. Readings were collected in all the interrow lines between the crops.

According to the manufacturer’s information, the EM38-MK2^® equipment simultaneously provides ECa measurements for 0.75 and 0.375 m in the horizontal dipole orientation. This depth reading was considered for use based on Jung et al. (2005), wherein it was demonstrated that it provides the best relationship between the ECa and the soil properties in layers to 30 cm, which is the depth that was contemplated through soil sampling.

The data were subjected to statistical treatments associated with the use of graphical tools, i.e., histograms and boxplots, to assess the shape and dispersion of the data set.

To reduce or eliminate overlap and to select the most representative forms of data from linear combinations of the variables from the soil analysis, the dimensionality of the data was reduced through principal component analysis (PCA) using R software (R Core Team, 2020). The variables were standardized for mean zero and unity variance. Here, we considered only the first two components since a single biplot could be obtained to summarize the results.

The package geoR version 1.8-1 (Ribeiro and Diggle, 2001) was selected for the geostatistical modeling of the principal components. A theoretical semivariogram was modeled by restricted maximum likelihood, or REML (Kerry and Oliver, 2007), and the tested models were spherical, exponential, and Gaussian. The model’s selection was based on the root mean square error (RMSE), coefficient of determination (R²), and mean squared deviation rate (MSDR) via cross-validation. Inferences about spatial dependence were based on the classification proposed by Seidel and de Oliveira (2016). This classification was considered to ensure an understanding of the degree of spatial dependence of the different principal components and to identify the intensity of the variability present in the experimental areas.

In this study, we utilized the inverse distance weighting (IDW) statistical interpolator (Shepard, 1968) since this type of interpolation is suitable when dealing with a high-density data set. After obtaining thematic maps from the information of the different principal components, the position and exact location of each experimental plot were determined (Supplementary Figure 2). After the plots on the thematic maps were overlaid, the values of each pixel (interpolation over 0.5 m × 0.5 m pixels) were extracted within each experimental plot, and after undergoing descriptive statistical analysis, the average value was considered and applied as a covariate in the genetic-statistical model, as detailed in equations 1.2 and 1.4.

Regarding the ECa sensor, due to the high density of data collected, all the raw data in both experimental areas were considered and directly plotted on the segmented plots; therefore, no interpolation procedure was required. Afterward, the average ECa of the soil was calculated for all the points obtained within each experimental plot. Each plot-specific mean ECa value was considered and applied as a covariate in the genetic-statistical model, as detailed in equations 1.2 and 1.4.

Both experiments were previously evaluated by following the same criteria routinely adopted in the early breeding stages of RIDESA/UFSCar; i.e., each family was measured for cane yield (tons of cane per hectare, TCH) when considering the plant cane stage after 12 months in May 2019. The TCH was obtained through mechanized harvesting and total plot weighing (family) with the support mobile truck-mounted weighing equipment. The estimates were calculated using the following equation: TCH = (TW × 10)/PS, where TW is the total weight of the plot (in kg) and PS is the total plot size in m², i.e., 75.6 m².

We used the linear mixed model approach where TCH was the response variable, and four different statistical models were assumed (from 1.1 to 1.4):

The four models have an intercept (μ), a random effect for i-th family (F_i), a fixed effect for the j-th block (b_j), and an error term.

The error term in models 1.1 and 1.2 assumes homogeneity over all the plots and no associations between them, statistically indicated as e_ij∼N(0,σ²). For models 1.3 and 1.4, the error term takes the spatial dependence over plots, i.e., the first-order autoregressive structure for rows and columns, or ε_ij∼N(0,AR(1)_uAR(1)_vσ²). The comparisons between the absence vs. the presence of spatial dependence (models 1.1 vs. 1.3 and 1.2 vs. 1.4) were performed using the Akaike information criterion (AIC) (Akaike, 1974) and the Bayesian information criteria (BIC) (Schwarz, 1978).

Models 1.2 and 1.4 included the environmental information (∑k=1wβk⁢di⁢j) captured using either soil sample information or ECa values, where β_k is the effect of the k-th environmental covariate, and d_ij is the covariate data for the plot containing the i-th family on the j-th block. The covariates inclusion was based on the forward selection approach (Kutner et al., 2004), i.e., (i) each variable was independently tested in the statistical model using the F-test; (ii) the covariates were ordered according to the F-test; (iii) if the highest was significant under 5%, the covariate was added to the model and the process was repeated to include the next one; otherwise, the process was stopped.

The residual variation coefficient (CVe%) was computed as the proportion of the residual error over the experimental mean for the TCH. The broad-sense heritability (H²) and the accuracy of selection (AC) were calculated according to the following equations:

where H²is the broad-sense heritability at the family mean level; σg2 is the variance in genetic effects between families; σe2 is the variance in residual effects (environmental); R is the number of repetitions for families; and AC is the accuracy of selection between families.

Before applying mixed models for data analysis, we verified the spatial distribution pattern of the raw data for each plot/family. We selected the R spatstat package version 1.63-3 (Baddeley et al., 2015), which has exploratory data analysis, model adjustment and simulation functionalities.

To perform the statistical analyses, GenStat software (Payne et al., 2011) was used to predict the genotypic values of the treatments (families) considering the different models, and R software (R Core Team, 2020) was selected for the geostatistics and graphical analyses.

PCA was performed for the 17 physical and chemical variables of the soil. The first two components captured 56.4% of the data variability (Supplementary Figure 3A). The first principal component (PC1) was mainly explained by variables related to soil acidity and base saturation in relation to the cation exchange capacity (CEC), i.e., base saturation, calcium, magnesium, percentage base saturation, percentage aluminum saturation, aluminum, and potential acidity; the only exception was the clay content (Supplementary Figure 3B). The second component (PC2) included a more diverse group of variables, such as pH, phosphorus, silt and the contents of different fractions of sand (total, coarse, and fine sand). The geostatistical analysis on the two principal components showed that, for both fields, PC1 present weak spatial dependence, while PC2 showed strong spatial dependence (Seidel and de Oliveira, 2016; Table 1). This might be explained by the main variables of each PC: most of the variables of PC1 are related to soil chemical properties (acidity), which tend to be more variable within the fields; in contrast, most of the variables of PC2 are related to soil texture (silt and sand contents), which often shows higher spatial dependence.

According to the map obtained by applying the IDW interpolation method (Figure 1), different variability patterns were observed for PC1 (Figure 1A). Negative values tended to represent regions where some parameters of PC1 had high positive expression, such as on the south side of the experiment. In this case, the area was most likely to have high levels of base saturation, calcium, magnesium, percentage base saturation, percentage aluminum saturation, aluminum, and potential acidity. Similar results were obtained when considering PC2 (Figure 1B), which was most associated with pH, phosphorus, silt and the contents of different fractions of sand. Each plot’s means for PC1 and PC2 values were obtained before proceeding with the subsequent analyses.

For the soil ECa, high-density data were obtained (11,722 reading points), which eliminated the need for interpolation. We used the mean value of ECa at the two depths (0.375 and 0.75 m) in each plot for the analyses (Figures 2A,B). Moreover, it is possible to note different patterns of variability, which shows that the soil is not uniform throughout the fields.

The exploratory analysis of the raw TCH data from EA1 indicated variations in the TCH values for the different experimental plots (families) throughout the experiment (Figure 3A). The highest yields are shown in dark blue (maximum = 160 tons/ha), and the lowest yields are indicated in red (minimum = 0 tons/ha). Only three experimental plots showed null results due to data loss (dark red). Figure 3B details the spatial patterns; for example, the yellowish regions located on the south side indicated higher productivity, whereas the bluish concentrated areas on the north side showed lower productivity, thus indicating a spatial variability tendency.

Considering the inclusion of environmental covariates (in models 1.2 and 1.4) in the genetic-statistical models, only PC1 was included using forward selection because, in the first round, it presented the highest significance (p = 0.011) (Supplementary Table 1). In the second round, no variable was significant (Supplementary Table 2).

When considering the four models, some insights can be highlighted. First, the models that took into account the spatial dependence between rows and columns (1.3 and 1.4) allowed better results than the models that considered only a homogeneous variance over plots (1.1 and 1.2). For example, models 1.3 and 1.4 showed the highest H² (0.60 and 0.62, respectively), the highest AC (0.77 and 0.79, respectively), the lowest CVe% (12.44 for model 1.3; 12.11 for model 1.4) (Table 2), and the lowest values of the AIC (4817.24 for model 1.3; 4810.63 for model 1.4) and BIC (4838.76 for model 1.3; 4822.14 for model 1.4) (Supplementary Table 3). Second, the inclusion of environmental covariates was similar for both models (1.2 and 1.4), i.e., the environmental effects (CVe%) were slightly reduced (13.40 for model 1.2; 12.11 for model 1.4). When considering the four models, the complete model (1.4) exhibited the best structure due to its capacity to model the H² environmental effects and increase the AC estimates (Table 2). When considering the different VCOV structures for EA1, the lowest values of the AIC and BIC (Akaike, 1974; Schwarz, 1978) were obtained for model 1.4 (Supplementary Table 3). Therefore, this model is the most recommended model for predicting the genotypic values of the study population.

The PCA performed with the 17 physical and chemical variables of the soil resulted in 44.1% of the variation being explained by the two first principal components (Supplementary Figure 4A). PC1 was essentially defined by the high positive values of aluminum saturation, aluminum, and clay variables, in contrast to the high negative values for the base saturation, calcium, percentage of base saturation, magnesium, pH, silt and cation exchange capacity variables. PC2 was mainly explained by the high positive values of the fine sand, percentage of base saturation, organic matter and pH variables, in contrast to the high negative values of the coarse sand, potential acidity and capacity to exchange cation variables (Supplementary Figure 4B). In contrast to what was observed in EA1, in this field, there was no clear division of the types of soil variables between PCs 1 and 2.

For EA2, the theoretical model of the semivariogram with the best fit to the data set was exponential for both PCs (Table 1). As observed for EA1, PC1 showed weak spatial dependence and PC2 showed strong spatial dependence. These results show that PC1 tends to present the soil variability in short-range distances, while PC2 often represents the long-range distance variability.

According to the map obtained by applying the IDW interpolation method (Figure 4), different variability patterns were observed for PC1 and PC2 (Figures 4A,B); i.e., the variability patterns are more likely to be associated with the previously cited variables.

For the soil ECa, 16,617 reading points were obtained within EA2 at depths of both 0.375 and 0.75 m (Figure 5). According to the data obtained for both depths, different patterns of variability were observed (Figures 5A,B).

The exploratory analysis suggested the presence of variations in the TCH throughout the experimental area (Figure 6A); the highest yields are indicated by the dark blue color (maximum = 350 tons/ha), and the smallest yields are indicated by the reddish color (minimum = 50 tons/ha). Sugarcane presented a higher yield potential in this field than in EA1. In this experiment, there were no missing data. Figure 6B details the spatial patterns.

When considering forward variable selection for genetic modeling, only the covariate ECa at 0.375 m (ECa05) was included in models 1.2 and 1.4, because in the first round, it presented a significance level below 0.05 (0.001 and 0.002, respectively); the other variables showed non-significant Pearson’s chi-square values (Supplementary Table 4). In the second round, the ECa1 variable became explanatory for both models, presenting chi-square values of 0.031 and 0.011 (Supplementary Table 5); the other variables were not explanatory, with values higher than 0.05.

The estimates of the variance components and genetic parameters of models 1.1, 1.2, 1.3, and 1.4 are shown in Table 3. For all models, the CVe presented values between 10 and 22%. A significant decrease was obtained when model 1.3 was considered. In addition, the inclusion of both covariates in this model (model 1.4) allowed for an expressive reduction in this parameter (CVe = 10.66%). For the total genetic variation coefficient (CVg), all models showed results above 10% for TCH. Nevertheless, the significant increase in this parameter was highlighted when models 1.3 and 1.4 were considered (Table 3). The broad-sense heritability increased for TCH when models 1.3 and 1.4 were considered (0.72 and 0.76, respectively). For the latter model, the highest accuracy of the selection value (AC = 0.87) was also observed (Table 3). Considering the different VCOV structures, with or without environmental variables, the lowest values obtained for the AIC and BIC were estimated for model 1.4 (Supplementary Table 3). Therefore, this model is the recommended model for predicting the genotypic values of the study population.