Multivariate Geostatistical Modeling of Phytophthora rubi and Pratylenchus penetrans in Red Raspberry Fields

Four commercial red raspberry fields were sampled for this study. One field was sampled in each Gresham, Oregon (OR-1) and Lynden, Washington (WA-1) in March 2014, and a second pair of fields was sampled in April 2015 in Portland, Oregon (OR-2) and Lynden, Washington (WA-2). The fields were chosen because they represented a diversity of raspberry growing environments in the region and also had areas that varied in plant vigor. All the fields were planted with the same raspberry variety, “Meeker”, and were naturally infested with P. rubi and P. penetrans. The ages of each field were OR-1 (23 years), OR-2 (12 years), WA-1 (5 years), and WA-2 (4 years). The soil type associated with each field was (USDA Soil Survey, 2024) OR-1 (Powell silt loam), OR-2 (Sauvie silt loam), WA-1 (Lynden sandy loam), and WA-2 (Whatcom silt loam).

An area of approximately 4 ha that contained plants with and without symptoms of Phytophthora root rot was delineated in each field. Within this area, a 6 × 10 uniform grid of 60 sampling points was established, and root and soil samples were collected at each grid point (Miller et al., 1997). GPS coordinates were recorded at each sampling point using a GeoXT 2005 series GPS unit (Trimble, Sunnyvale, CA). The data collected are represented as a set of marked points illustrated in algebraic terms: y=xi,mi,…,xn,mn y = \left\{ {\left( {{x_i},{m_i}} \right), \ldots ,\left( {{x_n},{m_n}} \right)} \right\}

Where x_i is the sampling point locations (Easting and Northing), and m_i is the attribute values (disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation) associated with each sampling point.

Due to differences in the field layout, the distance between each sampling point varied by field. Sampling points at OR-1 and OR-2 were 25 m × 35 m, while those at WA-1 and WA-2 were 35 m × 20 m and 45 m × 20 m, respectively (Fig. 1). Root samples were collected from a single plant at each sample site. Due to differences in plant health, not all samples had the same amount of collected root material, but generally, 5–10 g of root material was collected from each location. The next nearest plant was sampled if no plant was present at a designated sampling point. Approximately 500 g of soil from around the root system was collected using a shovel along with root material. Roots and soil from each sampling point were stored in the same bag and processed in the laboratory. Samples were passed over a 2-mm sieve to separate roots from the soil, and the soil sample was partitioned to facilitate different analyses.

Figure 1:

Layout schematic of a 6 × 10 grid sampling of contiguous quadrats (20 m × 45 m) in the commercial red raspberry production field WA-2 located in Lynden, WA. From the center of each quadrat (red dot), above-ground disease severity was visually assessed and rated; soil and root samples were collected to quantify Phytophthora rubi, Pratylenchus penetrans, and soil texture (sand, silt, clay); and surface field elevations were also recorded at each sampling point.

A 50-g subsample of soil was placed on a Baermann funnel to extract P. penetrans (Walters et al., 2017). Nematodes were collected after 5 days and enumerated using a dissecting microscope. A separate 50-g subsample of soil was placed in a 65°C oven for 5 days before being weighed to obtain soil moisture. Roots from each sample were rinsed with water to remove excess soil and then patted dry with a paper towel. Fine roots (<1 mm) were separated from larger roots, and approximately 1 g of fine roots was placed in a 20-ml glass scintillation vial and then stored at −80°C until further processing for P. rubi. The remaining fine roots were used for P. penetrans extractions, and larger roots were discarded. Fine root material was placed under intermittent mist for 5 days to extract P. penetrans (Walters et al., 2017). After nematode extraction, the root material was placed in a 65°C oven, and after 5 days, the dry weight was determined. P. penetrans were enumerated using a dissecting microscope.

Frozen fine root samples for P. rubi were processed by flash freezing in liquid nitrogen and then pulverized using a tissue pulverizer (Garcia Manufacturing, Visalia, CA) for 20 s. Pulverized roots were transferred into a 1.7 ml centrifuge tube and stored at −80°C until used for DNA extractions. DNA was then extracted from 100 mg of pulverized fine root tissue using a FastDNA SPIN Kit for Soil (MP Biomedicals, Santa Ana, CA) following the manufacturer’s instructions. Extracted DNA was then used to conduct a Phytophthora species-specific quantitative real-time polymerase chain reaction (qPCR) assay described by Bilodeau et al. (2014) with the following modifications. Assays were run on the StepOnePlus machine and software (Thermo Fisher Scientific, Waltham, MA) using PerfeCta Multiplex qPCR ToughMix (Quanta Biosciences, Gaithersburg, MD). Only probes for P. rubi and the plant internal positive control (IPC) were used. The target for the species-specific qPCR assay was the atp9-nad9 region of the mitochondrial genome. The amplification master mix contained the genus-specific primer pair at 0.5 μM and the species-specific P. rubi probe at 0.1 μM. The primers for the plant internal control were present at 0.0125 μM and the probe at 0.01 μM. The cycling conditions for the assay were set at 95°C for 2 min, 50 cycles at 95°C for 15 s, and 60°C annealing temperature for 1 min 30 s, in a reaction volume of 25 μl. Each qPCR run also contained a set of P. rubi DNA standards in the quantities of 2 ng, 0.2 ng, and 0.02 ng, and each standard or sample was run in triplicate (Bilodeau et al., 2014). Following the initial sampling of roots and soil in spring, each field was revisited in summer, and a visual disease rating assessment was performed for each plant at each sampling point. The ratings were based on a 1–10 scale, with 1 representing a healthy and vigorous plant with uniform growth, 5 representing a plant with moderate reduction in growth, and 10 representing an entirely dead plant (T. W. Walters, personal communication).

Soil texture was conducted by A & L Western Labs (Portland, OR). A subsample of soil from each sampling point was analyzed for sand, silt, and clay content.

Global Moran’s I was used to assess spatial autocorrelation in the distribution of disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation in the field (Rossi and Quénéhervé, 1998; Holguin et al., 2015). It was developed to test the null hypothesis of spatial independence or zero spatial autocorrelation (P > 0.05). The test statistic is a weighted product–moment correlation coefficient reflecting geographic proximity (Moran, 1950): I=nS0∑i=1n∑j=1nwi,jyi−y¯yj−y¯∑i=1nyi−y¯2 I = {n \over {{S_0}}}{{\sum\nolimits_{i = 1}^n {\sum\nolimits_{j = 1}^n {{w_{i,j}}\left( {{y_i} - \bar y} \right)\left( {{y_j} - \bar y} \right)} } } \over {\sum\nolimits_{i = 1}^n {{{\left( {{y_i} - \bar y} \right)}^2}} }}

Where y_i are the observations, w_i,j represents the distance weight (distance matrix), n is the number of observations, and S0=∑i=1n∑j=1nwi,j {S_0} = \sum\limits_{i = 1}^n {\sum\nolimits_{j = 1}^n {{w_{i,j}}} } is the sum of all w_i,j.

Global Moran’s I coefficient values range from −1 to 1, where a value <0 indicates negative spatial autocorrelation and a value >0 indicates positive spatial autocorrelation. A positive spatial autocorrelation indicates the presence of clustering, i.e., nearby areas have similar attribute values, whereas a negative spatial autocorrelation indicates dissimilarities in neighborhood data attribute values.

Local Moran’s I, also known as local indicators of spatial association (LISA), was used to identify local spatial clusters and outliers in the distribution of disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation for each field under the assumption of spatial heterogeneity (Anselin, 1995). The local spatial statistic Moran’s I is calculated for each sampling point based on a predefined spatial weight and is expressed as: Ii=yi−y¯∑k=1nyk−y¯2/n−1∑j=1nwijyj−y¯ {I_i} = {{\left( {{y_i} - \bar y} \right)} \over {\sum\nolimits_{k = 1}^n {{{\left( {{y_k} - \bar y} \right)}^2}/\left( {n - 1} \right)} }}\sum\limits_{j = 1}^n {{w_{ij}}\left( {{y_j} - \bar y} \right)}

Where a positive value for I_i indicates that a sampling point has a similar attribute value to the neighboring points (clusters), and a negative value indicates dissimilar attribute values in the local neighborhood (outliers).

Variogram modeling was used to determine and quantify spatial variability of disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation (Larkin et al., 1995; Miller et al., 1997; Madden et al., 2007). A variogram is the expected squared difference between two data values separated by a distance vector (Cressie, 1993; Wackernagel, 2003). The semivariogram represents half the variogram and was first introduced by Matheron (1963). γd=12∑zxi−zxj2 \gamma \left( d \right) = {1 \over 2}\sum {\left\{ {{{\left[ {z\left( {{x_i}} \right) - z\left( {{x_j}} \right)} \right]}^2}} \right\}}

Where γ(d) is the semivariance values as a function of lag distance and z(x_ij) is the point locations associated with their attribute values.

Parameters associated with the variogram include nugget (the short-range variability in the data), sill (the total variance of the empirical variogram), and range (the distance after which the data are no longer correlated). Variogram points above the sill indicate negative spatial autocorrelation and points below the sill indicate positive spatial autocorrelation. The variogram is fitted with a specific known mathematical function such as the spherical, exponential, Gaussian, or Matern models.

The ordinary kriging was used as a probabilistic interpolator to predict disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation at one or more unsampled points. It is considered the best unbiased linear interpolation method, where the expectation of errors equals zero, and the predicted value is a linear combination of weighted known values (Cressie, 1990): z^x0=∑i=1nλizxi \hat z\left( {{x_0}} \right) = \sum\nolimits_{i = 1}^n {{\lambda _i}z\left( {{x_i}} \right)}

Where λ_i are the weights and z(x_i) is the known value of the variable z at sampling point x_i.

Kriging assumes spatial variation to be statistically consistent across the field and considers that an attribute value is neither completely random nor completely determined. Prediction is codetermined by spatial autocorrelation factors, offsets, and random error (Cressie, 1988, 1990).

Ordinary least squares (OLS) multiple regression is a classical statistical method routinely used to predict dependent variables measured at the interval or ratio level. The OLS multiple regression model is expressed as (Ott and Longnecker, 2016): y=β0+β1X1+β2X2+…+βkXk+ε y = {\beta _0} + {\beta _1}{X_1} + {\beta _2}{X_2} + \ldots + {\beta _k}{X_k} + \varepsilon

Where y is the predicted value of the dependent variable, X₁ through X_k are k distinct independent variables, β₀ is the intercept, the value of y when all the independent variables are equal to zero, and β₁ through β_k are the estimated regression coefficients.

In OLS, the estimates and standard errors (SE) are unbiased if the required assumptions are fulfilled—linearity, homoscedasticity, independence of errors, normality, and independent variables (Ott and Longnecker, 2016; James et al., 2021). With spatial dependency, using OLS will produce erroneous predictions due to spatial autocorrelation in the predictor or the OLS regression model residuals, leading to biased parameter estimates (Lopez-Nicora et al., 2020). Therefore, spatial regression models are appropriate for accommodating spatial dependency in either the response variable or the regression model errors (Bivand et al., 2013). The spatial lag model (SLM), spatial Durbin model (SDM), spatial error model (SEM), conditional autoregressive (CAR) model, and spatial autoregressive moving average (SARMA) are widely used to account for spatial interactions between autocorrelated dependent and independent variables (Anselin, 1988; Madden et al., 1988; Bivand et al., 2013; Plant, 2019) and were therefore used in this study.

Spatial lag model, also known as the spatial autoregressive model, only considers the effect of the spatial lag on the dependent variable while focusing on the interaction among observations (Whittle, 1954; Anselin, 1988, 1992) and is defined as: y=β0+ρwiyj+β1X1+β2X2+…+βkXk+ε y = {\beta _0} + \rho {w_i}{y_j} + {\beta _1}{X_1} + {\beta _2}{X_2} + \ldots + {\beta _k}{X_k} + \varepsilon

Where ρ is the spatial coefficient and w_iy_i is the spatially lagged dependent variable for the weights matrix w_i. The parameter ρ also specifies the degree of spatial dependence; when it is positive, a positive correlation exists; and when negative, a negative correlation exists. If there is no spatial dependence on the independent variables, y does not depend on y_i, then ρ = 0, and the equation reverts to a classical OLS multiple regression.

Spatial Durbin model is a modification of SLM and considers the spatial lag’s effect on the independent and dependent variables (LeSage and Pace, 2009) and is defined as: y=β0+ρwiyj+w1Xβ2+β1X1+β2X2+…+βkXk+ε y = {\beta _0} + \rho {w_i}{y_j} + {w_1}X{\beta _2} + {\beta _1}{X_1} + {\beta _2}{X_2} + \ldots + {\beta _k}{X_k} + \varepsilon

Where w₁Xβ₂ is the spatial lag effect of the independent variables.

Spatial error model considers the spatial lag’s effect in the error term and does not include lagged dependent or independent variables (Anselin, 1992). The error is split into random error and spatially structured error and is defined as: y=β0+λwiεi+β1X1+β2X2+…+βkXk+ε y = {\beta _0} + \lambda {w_i}{\varepsilon _i} + {\beta _1}{X_1} + {\beta _2}{X_2} + \ldots + {\beta _k}{X_k} + \varepsilon

Where λ is the spatial error coefficient and w_iɛ_i is the spatially lagged error term for weights matrix w_i. If there is no spatial correlation between the random and spatially structured errors, then λ = 0.

Conditional autoregressive focuses on the conditional distribution of the spatial error terms (Cliff and Ord, 1981; Schabenberger and Gotway, 2005; Bivand et al., 2013) and is defined as: ei|ej~i~N∑j~icijej∑j~icij,σei2∑j~icij {e_i}|{e_{j\sim i}}\sim N\left( {\sum\limits_{j\sim i} {{{{c_{ij}}{e_j}} \over {\sum\nolimits_{j\sim i} {{c_{ij}}} }},} {{\sigma _{{e_i}}^2} \over {\sum\nolimits_{j\sim i} {{c_{ij}}} }}} \right)

Where e_i is the conditional distribution, N denotes a normal distribution with the mean and variance parameters, σei2 \sigma _{{e_i}}^2 is the conditional variance, and c_i,j is the spatial dependence between point locations.

Spatial autoregressive moving average is used for modeling local autocorrelation and can account for spatial dependence among the error terms, like SEM, and the dependent variable, like SLM (LeSage and Pace, 2009), an d is defined as: y=β0+ρwiyj+λwiεi+β1X1+β2X2+…+βkXk+ε y = {\beta _0} + \rho {w_i}{y_j} + \lambda {w_i}{\varepsilon _i} + {\beta _1}{X_1} + {\beta _2}{X_2} + \ldots + {\beta _k}{X_k} + \varepsilon

This study used the software R version 4.4.0 as a modeling language environment for geostatistical analysis and disease mapping (R Core Team, 2018). Exploratory data analysis was performed, and values of disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation were summarized for each field. P. rubi DNA values and P. penetrans in soil and root were normalized with Log₁₀-transformation. Arcsine transformations were applied to disease rating, sand, silt, and clay, and elevation values were square root transformed. The data were transformed into a spatial data frame containing point locations and associated attribute values (disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation) using the function “SpatialPointsDataFrame” within the package “sp” (Pebesma and Bivand, 2005). Spatial matrices and weights were created based on nearest neighbors—with 15 links on average—using the functions “knearneigh,” “knn2nb,” and “nb2listw” from the package “spdep” (Bivand et al., 2024a, 2024b).

Global Moran’s I test for spatial autocorrelation was computed using the function “moran.test” within the package “spdep” (Bivand et al., 2024a, 2024b), and the Global Moran’s I coefficient and P-value were extracted and compiled in a summary table for each field. Local Moran’s I test for spatial association was computed using the function “localmoran” from the package “spdep” (Bivand et al., 2024a, 2024b). Variogram analysis was performed using the function “autofitVariogram.” The function is available within the package “automap” (Hiemstra and Skoien, 2023). Directional variogram analysis was performed using the function “variogram” from the package “gstat” to assess for anisotropy in the spatial distribution of disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation across the field (Pebesma and Graeler, 2023). The ordinary kriging was performed to map disease rating, P. rubi DNA in root, P. penetrans in root and soil, soil texture, and elevation using the function “autoKrige” in the package “automap” (Hiemstra and Skoien, 2023).

Ordinary least squares multiple regression was used to assess the influence of soil texture and elevation on the distribution of P. rubi DNA in root and P. penetrans in root and soil for each field using the function “lm.” Additionally, the effects of P. rubi DNA in root and P. penetrans in root and soil on disease severity rating and the relationship between P. penetrans in soil and in roots were assessed in regression modeling. Homoscedasticity was assessed using the residuals from the fitted models and plotted against the predicted values, and multicollinearity was assessed using the model’s variance inflation factor (VIF) (James et al., 2021). Moran’s I test for spatial autocorrelation in the residuals was performed from the estimated linear model using the function “lm.morantest” from the package “spdep” (Bivand et al., 2024a, 2024b). Finally, SLM and SDM were computed using the function “lagsarlm,” SEM with the function “errorsarlm,” and CAR and SARMA with the function “spautolm” within the package “spatialreg” (Bivand et al., 2024a, 2024b). Akaike information criterion (AIC) was used to compare OLS, SLM, SDM, SEM, CAR, and SARMA models and to determine the best model that fits the data for each field (James et al., 2021). Subsequently, the model with the lowest AIC value was considered the best.

Disease severity was worse in OR-1 and WA-2 than in OR-2 and WA-1 (P < 0.05). The variance-to-mean ratio (VMR) values showed that the distribution of disease severity values was under-dispersed for all the fields. P. rubi DNA concentrations in roots varied significantly across all four fields (P < 0.05) (Table 1). WA-1 had the highest average P. rubi DNA concentrations in roots, followed by WA-2, OR-1, and OR-2, and VMR values were highest in WA-1 and WA-2, indicating the presence of considerable variations in the distribution of P. rubi DNA concentrations in roots (Table 1). P. penetrans densities in roots were significantly higher in WA-2 and OR-2 and were the lowest in WA-1 and OR-1 (P < 0.05), with VMR values indicating strong variations in their distributions for each field (Table 1). P. penetrans densities in soil were significantly higher in WA-2 and more dispersed compared with the other fields (P < 0.05) (Table 1). Soil texture varied greatly across all four fields; WA-1 had the highest content of sand, while OR-1 and WA-2 had the highest silt content, and OR-2 had the highest content of clay (Table 1). For field elevation, only OR-1 had a significant elevation level compared to the other fields (P < 0.05) (Table 1).

Global Moran’s I test was used to assess spatial autocorrelation in the distribution of the measured variable (Table 2). Disease severity and P. rubi DNA concentrations in roots had significant positive spatial autocorrelations in OR-1, WA-1, and WA-2 (P < 0.05). P. penetrans densities in roots had significant positive spatial autocorrelation only for WA-2 (P < 0.05), while P. penetrans densities in soil had significant positive spatial autocorrelation for OR-2, WA-1, and WA-2 (P < 0.05). Soil texture was positively significant for all fields except for clay content in OR-2 (P < 0.05), and field elevations were positively significant for all fields (P < 0.05).

Local Moran’s I test was performed to identify relationships and spatial associations between each observation of the measured variable with their surrounding neighbors. There were significant similarities or positive spatial associations (PSA) to neighboring sampling points had in 37% of the sampled grids for OR-1, 5% for OR-2, 30% for WA-1, and 72% for WA-2 for disease severity (P < 0.05) (Fig. 2). The PSA was significant in 12% of the sampled grids for OR-1, 0% for OR-2, 8% for WA-1, and 45% for WA-2 for P. rubi DNA in roots (P < 0.05) (Fig. 2). For P. penetrans in roots, the PSA was significant in 5% of the sampled grids for OR-1, 3% for OR-2, 0% for WA-1, and 15% for WA-2 (P < 0.05) (Fig. 2). The PSA was significant for P. penetrans in soil in 14% of the sampled grids for OR-1, 12% for OR-2, 44% for WA-1, and 10% for WA-2 (P < 0.05) (Fig. 2). The Moran scatterplot showed positive slopes for the regression lines for disease severity in all the fields; P. rubi DNA in roots in OR-1, WA-1, and WA-2; P. penetrans in roots in OR-1, WA-1, and WA-2; and P. penetrans in soil in all the fields (Fig. 3). The PSA was significant for soil texture, sand contents in 22% of the sampled grids for OR-1, 7% for OR-2, 52% for WA-1, and 67% for WA-2 (P < 0.05) (Fig. 4). Silt contents had significant PSA in 15% of the sampled grids for OR-1 and OR-2, 43% for WA-1, and 32% for WA-2 (P < 0.05) (Fig. 4). The PSA was significant for clay content in 31% of the sampled grids for OR-1, 3% for OR-2, and 53% for WA-1 and WA-2 (P < 0.05) (Fig. 4). Field elevation had significant PSA in 80% of the sampled grids for OR-1, 24% for OR-2, 68% for WA-1, and 63% for WA-2 (P < 0.05) (Fig. 4). There were positive slopes for the regression lines for soil texture and field elevation across all the fields according to the Moran scatterplot (Fig. 5).

Figure 2:

LISA for disease rating, Phytophthora rubi DNA concentrations in roots, and Pratylenchus penetrans densities in roots and soil in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). A positive value for the local Moran’s I indicates that a sampling point has a similar attribute value to the neighboring points (clusters), and a negative value indicates dissimilar attribute values in the local neighborhood (outliers). Solid-filled black shapes indicate the significance of the local Moran’s I coefficient at P < 0.05. LISA, local indicators of spatial association.

Figure 3:

Moran scatterplot for disease rating, Phytophthora rubi DNA concentrations in roots, and Pratylenchus penetrans densities in roots and soil in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The scatterplot is characterized by four quadrants corresponding to the four types of spatial association and illustrates relationships between an individual sampling value (standardized observation) and the average values at neighboring sampling points (standardized spatial lags). The lower left and upper right quadrants indicate spatial clustering or PSA. The lower left quadrant is characterized by low values for both the spatial lags and the observations (Low–Low), and in the upper right quadrant, there are high values for both the spatial lags and the observations (High–High). The upper left (Low–High) and lower right (High–Low) quadrants indicate spatial outliers or NSA. The slope of each regression line corresponds to the Moran’s I coefficients listed in Table 2. NSA, negative spatial association; PSA, positive spatial association.

Figure 4:

LISA for soil texture (sand, silt, and clay) and field elevation in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). A positive value for the local Moran’s I indicates that a sampling point has a similar attribute value to the neighboring points (clusters), and a negative value indicates dissimilar attribute values in the local neighborhood (outliers). Solid-filled black shapes indicate the significance of the local Moran’s I coefficient at P < 0.05. LISA, local indicators of spatial association.

Figure 5:

Moran scatterplot for soil texture (sand, silt, and clay) and field elevation in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The scatterplot is characterized by four quadrants corresponding to the four types of spatial association and illustrates relationships between an individual sampling value (standardized observation) and the average values at neighboring sampling points (standardized spatial lags). The lower left and upper right quadrants indicate spatial clustering or PSA. The lower left quadrant is characterized by low values for both the spatial lags and the observations (Low–Low), and in the upper right quadrant, there are high values for both the spatial lags and the observations (High–High). The upper left (Low–High) and lower right (High–Low) quadrants indicate spatial outliers or NSA. The slope of each regression line corresponds to the Moran’s I coefficients listed in Table 2. NSA, negative spatial association; PSA, positive spatial association.

Experimental variogram analyses were performed to evaluate the spatial variability and fit a spatial model to the distribution of the measured variable. Disease severity displayed some increased levels in the spatial continuity when plotting the distance lags against the semivariance for OR-1 and WA-2 (Fig. 6), and Gaussian and exponential models were fitted to their respective variograms (Table 3). P. rubi DNA concentrations and P. penetrans densities in roots were mostly characterized by steady and constant trend lines when assessing the semivariance over lag distances of pairwise points. P. rubi spatial continuity was more pronounced in WA-1, and P. penetrans in roots was highly defined in WA-2 (Fig. 6), and both were fitted with an exponential model (Table 3). P. penetrans densities in soil had a defined spatial continuity in WA-1 and WA-2 (Fig. 6); spherical and exponential models were fitted to their respective variograms (Table 3). There was strong spatial dependency for sand and silt contents in WA-1 and WA-2, for clay content in OR-2, and for field elevation in OR-1, WA-1, and WA-2 (Fig. 7).

Figure 6:

Experimental variogram graph for assessing the spatial structure of disease rating, Phytophthora rubi DNA concentrations in roots, and Pratylenchus penetrans densities in roots and soil in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The x-axis shows the distance lag, which is created by plotting the average variability between pairs of sampling points, and the y-axis shows the calculated value of the semivariance, where a greater value indicates less spatial correlation between pairs of points. The variogram model type and fitting parameters for each commercial field are listed in Table 3.

Figure 7:

Experimental variogram graph for assessing the spatial structure of soil texture (sand, silt, and clay) and field elevation in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The x-axis shows the distance lag, which is created by plotting the average variability between pairs of sampling points, and the y-axis shows the calculated value of the semivariance, where a greater value indicates less spatial correlation between pairs of points. The variogram model type and fitting parameters for each commercial field are listed in Table 3.

Directional variogram analyses were performed to test for the potential presence of anisotropy in the distribution of the measured variable. Overall, there were no major directional variations in the distribution of disease severity, P. penetrans densities in soil, and sand and silt contents. However, the presence of moderate anisotropy in the distribution of P. rubi DNA in roots in OR-2 and WA-2 and with P. penetrans densities in roots for WA-1 was detected (Fig. 8). There was moderate anisotropy for clay content in OR-2, WA-1, and WA-2. Field elevation had high levels of anisotropy for all the locations (Fig. 9).

Figure 8:

Directional experimental variogram graph for assessing potential directional variations in the spatial structure of disease rating, Phytophthora rubi DNA concentrations in roots, and Pratylenchus penetrans densities in roots and soil in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). Directional angles used for the graph include 0° (North–South), 45° (Northeast–Southwest), 90° (East–West), and 135° (Southeast–Northwest). The overlapping of all the angle point lines indicates no directional variations in the spatial structure of the attribute values (isotropy), while major gaps between the angle point lines indicate the presence of directional trends (anisotropy).

Figure 9:

Directional experimental variogram graph for assessing potential directional variations in the spatial structure of soil texture (sand, silt, and clay) and field elevation in four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). Directional angles used for the graph include 0° (North–South), 45° (Northeast–Southwest), 90° (East–West), and 135° (Southeast–Northwest). The overlapping of all the angle point lines indicates no directional variations in the spatial structure of the attribute values (isotropy), while major gaps between the angle point lines indicate the presence of directional trends (anisotropy).

Ordinary kriging was used as a probabilistic interpolator to estimate the values of the measured variable at unsampled locations in the fields. Before mapping the kriging predictions, density plots were used to capture and compare the distribution of the predictive outcomes. Disease severity values were skewed toward high disease levels for OR-1, restricted between lower and middle values for OR-2 and WA-1, and a bimodal distribution of two distinct peaks of low and high disease values for WA-2 (Fig. 10). P. rubi DNA concentrations in roots had a trimodal distribution centered on the average DNA values for OR-1; a single distribution peak of low DNA values for OR-2; a prominent positive skewness ranging from low to high DNA values for WA-1; and finally, a bimodal distribution of DNA values for WA-2 (Fig. 10). P. penetrans density in roots were restricted to low values dominated by a single distributional and overlapping peak for OR-1 and WA-1; a trimodal distribution centered slightly above the average values of nematodes in roots was recorded for OR-2; and WA-2 was dominated by a major and multiple minor peaks of high nematode values in roots (Fig. 10). P. penetrans densities in soil were dominated by a single major peak of distribution for OR-1 and OR-2, where OR-1 encompassed lower values while OR-2 higher values of nematodes in soil; an extensive multimodal distribution and skewed to lower values of nematodes in soil for WA-1; and finally, WA-2 was dominated by a single peak of high distribution values of nematodes in soil (Fig. 10).

Figure 10:

Kriging density graph for comparing the distribution of interpolated values of disease rating, Phytophthora rubi DNA concentrations in roots, and Pratylenchus penetrans densities in roots and soil across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The x-axis shows the attribute values of the variable, and the y-axis shows the proportion value assigned to each range of the attribute values.

Sand content was characterized by low distribution values for OR-1, ranging between 25% and 30%, and OR-2, between 23% and 35%; WA-1 had the highest range of sand content (>55%) with a bimodal distribution; and finally, WA-2 encompassed a bimodal distribution centered to lower range values and an extended skewness reaching 50% of sand content (Fig. 11). Silt content was characterized by a cluster of multimodal density distributions in the upper level of high values for OR-1, OR-2, and WA-2, and WA-1 had the lowest range of values for silt content with a prominent bimodal distribution (Fig. 11). Clay content was characterized by a clear segmentation of OR-1 and OR-2, where OR-2 had the highest distribution levels of clay content; WA-1 dominated the lower range values for clay content; and WA-2 had a trimodal distribution and was skewed toward high values of clay content (Fig. 11). Field elevation was characterized by a cluster of low distribution values for OR-2, WA-1, and WA-2, and OR-1 had the highest distribution levels of elevation values with an extensive skewness toward low values (Fig. 11).

Figure 11:

Kriging density graph for comparing the distribution of interpolated values of soil texture (sand, silt, and clay) and field elevation across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The x-axis shows the attribute values of the variable, and the y-axis shows the proportion value assigned to each range of the attribute values.

Predictive maps generated by ordinary kriging analysis highlighted the distribution of the measured variables across the fields. Disease severity distribution was characterized by a single focus located westside of OR-1 with a concentric distributional ripple effect expanding eastbound. In OR-2 a monotonous distribution of disease severity values was observed. WA-1 had a descending wave of disease severity values from north to south. WA-2 was characterized by an elongated clump of high disease severity values, stretching from northeast to southeast with radiating waves toward both west and east sides of the field (Fig. 12). P. rubi DNA concentrations in roots followed similar distribution patterns to that of disease severity for OR-1 with overlapping ripple effect expanding eastbound. In OR-2 a quasi-monotonous distribution of low fungal DNA values in roots was observed. WA-1 was characterized by a major focus of high values of fungal DNA in roots and was located slightly toward the northwest side of the field and surrounded by several satellite clusters located toward the southeast, southwest, and west corner of WA-1. In WA-2, quasi-similar distribution patterns similar to disease severity were observed with an elongated focus of high fungal DNA values, stretching from the northeast to the southeast sides of WA-2 (Fig. 13).

Figure 12:

Spatial distribution of disease rating across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict disease ratings at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Figure 13:

Spatial distribution of Phytophthora rubi DNA concentrations in roots across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict P. rubi DNA in roots at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

In OR-1, P. penetrans densities in roots had a quasi-monotonous distribution of low nematode values with an elongated focus of low densities stretching from west to east. OR-2 is dominated by a large, centered clump of low P. penetrans densities in roots. WA-1 were characterized by a disorganized and nonstructured distribution pattern of low and high nematode densities in roots. In WA-2, quasi-similar distribution patterns similar to disease severity and P. rubi DNA concenrtrations in roots were observed, where high P. penetrans densities in roots were mostly found stretching from the northeast to the southeast sides of WA-2 (Fig. 14). In OR-1 and OR-2, P. penetrans in soil had the presence of multiple foci with low to median densities. WA-1 was characterized by a single clump of high nematode densities in soil with a rippling effect in all directions; however, large areas of extremely low nematode densities predominated northwest and northeast sides of WA-1. WA-2 was characterized by high P. penetrans densities in soil and predominantly circumscribed at the northeast, east, and southeast edges of the field, and a defined rippling wave of low P. penetrans densities was located slightly toward the center of WA-2 (Fig. 15).

Figure 14:

Spatial distribution of Pratylenchus penetrans densities in roots across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict P. penetrans densities in roots at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Figure 15:

Spatial distribution of Pratylenchus penetrans densities in soil across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict P. penetrans in the soil at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Sand content had a quasi-monotonous distribution with a limited focus of low sand values for OR-1. OR-2 was characterized by a large clump of low sand values stretching from the westside to the center of the field, and median high values of sand were observed toward the northwest and southeast sides of OR-2. In WA-1, there was a massive area of extremely high sand content located southside of the field and radiating northward. In WA-2, quasi-similar distribution patterns of sand content similar to that of disease severity, P. rubi DNA concentration in roots, and P. penetrans densities in roots were observed with a stretching area of high sand content from northeast to southeast sides of WA-2 (Fig. 16). Silt content was predominated by a vast area amounting to over 50% of silt distribution for OR-1. OR-2 had multiple large swaths of areas encompassing silt content values ranging between 40% and 50%. WA-1 was characterized by mostly low values for silt content (22%–28%) located mostly in the southside of WA-1; however, a single clump of high silt content was observed in the northeast corner of the field. In WA-2, there were multiple foci of high silt content located northwest and southwest sides of the field (Fig. 17). Clay content was characterized by a large focus of high clay values located northwest corner of OR-1 and radiating southward. OR-2 had the highest level of clay content compared to the other fields and was characterized by a large centered clump spanning from north to south. WA-1 had a three-layered area encompassing low and median values for clay content. In WA-2, there was a focus of low values in the east and northeast sides and median values in the west and southwest corners of the field (Fig. 18).

Figure 16:

Spatial distribution of the percentage of sand across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict the percentage of sand at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Figure 17:

Spatial distribution of the percentage of silt across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict the percentage of silt at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Figure 18:

Spatial distribution of the percentage of clay across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict the percentage of clay at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Field elevation was predominated by a focus of high elevation values located southwest corner of OR-1 and radiating northeastward. In OR-2, there was a quasi-monotonous elevation distribution that was characterized by a single focus located north side of the field. WA-1 was represented by a four-layered area spanning from north to south with increasing elevation values. WA-2 followed similar patterns with a three-layered area spanning from northwest to southeast with increasing elevation values (Fig. 19).

Figure 19:

Spatial distribution of field elevation across four commercial red raspberry production fields in Oregon (OR-1 and OR-2) and Washington (WA-1 and WA-2). The ordinary kriging was used as a probabilistic interpolator to predict field elevation at one or more unsampled points. Northing divides the map from north to south and Easting from west to east. Northing and Easting are grid references and are expressed in meters.

Ordinary least squares multivariate regression analyses were performed for each field location before performing Moran’s I test to detect spatial autocorrelation in the residuals. P. rubi DNA concentrations in roots was not significantly influenced by sand, silt, clay, and elevation for OR-1, OR-2, and WA-2 (P > 0.05), and a severe multicollinearity effect was observed for WA-1 due to overinflated variance (Table 4). Spatial autocorrelation in the residuals was significant only for WA-2 (P < 0.05) (Table S1 in Supplementary Materials). Based on a comparative regression model analysis, the SARMA model was able to detect significant effects of sand, silt, and elevation on P. rubi DNA concentrations in roots for OR-1 (P < 0.05; AIC = 179.81, R² = 0.27). Additionally, the SARMA model revealed a significant influence of sand, silt, clay, and elevation on P. rubi DNA concentrations in roots for WA-1 (P < 0.05; AIC = 192.91, R² = 0.33). OR-2 and WA-2 had the lowest AIC values based on the SDM model; however, there were no significant effects of soil texture and elevation on P. rubi DNA in roots (P > 0.05) (Table 4).

Pratylenchus penetrans densities in roots were not significantly influenced by soil texture and field elevation in any of the fields when using OLS multivariate regression analysis (P > 0.05) (Table 5). Severe multicollinearity effects were observed for soil texture across all fields except OR-2. No significant differences were observed for the spatial autocorrelation in the residuals for all fields (P > 0.05) (Table S1 in Supplementary Materials). For OR-1 and WA-2, the SARMA model revealed significant effects of field elevation on the distribution of P. penetrans densities in roots (P < 0.05; AIC = 0.31, R² = 0.42). For OR-2 and WA-1, the OLS model had the lowest AIC values compared to the other models; however, soil texture and elevation had no significant influence on P. penetrans densities in roots (P > 0.05) (Table 5). P. penetrans densities in soil was not significantly influenced by soil texture and field elevation in none of the fields when using OLS (P > 0.05), and severe multicollinearity effects were observed across all fields except OR-2 (Table 6). Spatial autocorrelation in the residuals was significant for OR-2, WA-1, and WA-2 (P < 0.05) (Table S1 in Supplementary Materials). WA-1 was the only field in which sand, silt, and clay had a significant effect on P. penetrans densities in soil using the SDM model (P < 0.05; AIC = 111.71, R² = 0.59) (Table 6).

Disease severity was significantly influenced by P. rubi and P. penetrans in soil for OR-1, WA-1, and WA-2 (P < 0.05) in the OLS model, and no severe multicollinearity effects were observed (Table 7). However, spatial autocorrelation in the residuals was significant for OR-1, WA-1, and WA-2 (P < 0.05) (Table S2 in Supplementary Materials). For OR-1 and WA-2, the SDM model had the lowest AIC values, −96.37 and −98.44, respectively, and P. rubi significantly affected disease severity (P < 0.05) with R² values, 0.61 and 0.73, respectively. Additionally, P. penetrans in soil had a significant effect on disease severity in WA-2 (P < 0.05). For OR-2, the OLS model had the lowest AIC values, −153.46, compared to the other models. However, no significant effects of P. rubi and P. penetrans were observed (P > 0.05). Finally, for WA-1, the SLM model had the lowest AIC values, −56.27, and P. rubi had a significant effect on disease severity (P < 0.05) with an R² value of 0.43 (Table 7).

Pratylenchus penetrans densities in roots were significantly influenced by P. penetrans densitiesin soil for OR-1, OR-2, and WA-2 (P < 0.05) when using the OLS model, and no severe multicollinearity effects were observed (Table 8). No significant differences were observed for the spatial autocorrelation in the residuals across all fields (P > 0.05) (Table S3 in Supplementary Materials). For OR-1, the SLM and SARMA models had the same lowest AIC values, 128.26, and a similar R² value of 0.23, with P. penetrans densities in soil significantly affecting P. penetrans in the roots (P < 0.05). For OR-2, the SARMA model had the lowest AIC value, 189.57, an R² value of 0.33, and P. penetrans densities in the soil significantly affected P. penetrans in the roots (P < 0.05). For WA-1, the OLS model had the lowest AIC value, 134.05. However, P. penetrans densities in soil and P. rubi DNA concentrations in roots had no significant effect on P. penetrans densitiesin the roots (P > 0.05). For WA-2, the SDM model had the lowest AIC value, 183.53. However, P. penetrans densities in soil and P. rubi DNA concentrations in roots had no significant effect on P. penetrans densities in the roots (P > 0.05) (Table 8).