Spatial distribution and risk area assessment of Aphelenchoides besseyi using geostatistical approaches in Giridih district of Jharkhand, India

This study was conducted in Giridih district of Jharkhand (Fig. 1), India during kharif season (rainy season, July-November) in the year 2015. There are 12 administrative blocks in the Giridih district. Among the three major agro-ecological zones of Jharkhand, Giridih belongs to the Central and North Eastern Plateau Zone with an average altitude of 968 feet and annual rainfall of 1,128 mm. Majority of the growing areas are predominantly irrigated or rainfed upland in nature, moreover there are no flood prone areas in the study region. Rice is the main cereal crop covering 68.3% of total agricultural land. (http://agricoop.nic.in/sites/default/files/JKD18-Giridh-31.03.2013.pdf).

Figure 1:

Nematode sampling was conducted across the Giridih during months of October to November, coinciding with the ripening phase of rice. A minimum distance of 1 km was maintained between two sample points. From each block 12 to 18 samples were randomly collected following ‘W’ pattern sample walk method (Karuri et al., 2017). At each site, 30 matured rice panicles were randomly sampled from three adjoining rice fields. Each sampling site was geopositioned and collected panicles were labeled and stored in paper packets before bringing back to the laboratory. In total, 163 samples were collected across the district (Fig. 1C).

Nematodes were isolated using the modified Baermann Funnel method (Schindler, 1961). Briefly, from each sampling site, 100 grains were randomly separated from the collected panicles and pounded with motor and pastel. The grinded material was placed over tissue paper wire gauge assembly placed on a Petri plate filled with water. This assembly was covered with another Petri plate to minimize evaporation loss and kept for 24 hr at room temperature (around 30°C). Extracted nematodes were killed in hot water bath at 65°C for 5 min and fixed in TAF fixative for further analysis (Shepherd, 1970). For identification, fixed nematode samples were processed by Seinhorst’s glycerol-ethanol method (Seinhorst, 1959) and finally mounted on glass slides according to De Grisse (1969). Species identification was done based on morphology and morphometrics key parameters by Khan et al. (2012) using Zeiss Axio ScopeA1. Population density (no. of nematode/100 grains) was measured from 5 ml aliquots of fixed nematode suspension with three replications (Merny and Luc, 1969). Final nematode count was made by multiplying average count (n =3) with the volume of nematode suspension (Mondal et al., 2018).

Population densities of RWTN collected across different blocks were subjected to the Kruskal–Wallis test in R statistical software (version R-3.6.1) to examine if any significant variation exists between administrative blocks of Jharkhand state. Box and whisker plot representing population densities (no. of nematode/100 grains) of A. besseyi across different blocks of Giridih was prepared using ‘SigmaPlot 14.0’ software. Agglomerative hierarchical cluster analysis using average linkage method was performed based on population density of RWTN to infer the distance among the blocks (Kaufman and Rousseeuw, 2009). Normalization of data followed by cluster analysis was performed using ‘hclust’ function in R. In average linkage hierarchical clustering, the distance (L) between two cluster (r, s) is defined as the average distance between each point of a cluster to every points of the other cluster and can be expressed as follows: (1) (1) L ( r , s ) = 1 n r n s ∑ i = 1 n r ∑ j = 1 n s D ( x r i , x s j , \text{(1)} & L(r,s)=\frac{1}{{n}_{r}{n}_{s}}\sum _{i=1}^{{n}_{r}}\sum _{j=1}^{{n}_{s}}D({x}_{ri},\,{x}_{sj},\end{array}\end{eqnarray}?> (1)where x and y are the observations from clusters r and s, respectively. Silhouette analysis was performed using ‘cluster’ package in R. Silhouette width (S _i ) was used to validate the number of clusters formed. Silhouette width was interpreted as follows: S _i  = almost 1; very well clustered, S _i  = around 0; observation lies between two cluster and S _i  = 0; probably placed into wrong cluster.

Spatial pattern of RWTN distribution in Giridih district was characterized by utilizing two different geospatial statistical techniques, point pattern and surface interpolation analyses. First, to identify if any significant hotspots of RWTN infestation exists in the district, point pattern optimized hotspot analysis was carried out. Second, inverse distance weighting (IDW) and Kriging were used to generate continuous map of predicted surface of nematode distribution across Giridih. Finally, risk areas of RWTN (population density > ETL) in the district was delineated by utilizing the indicator kriging (IK) tool, which is also a surface interpolation method. All these operations were performed using ArcMap 10.2, ESRI.

As nematode density data was positively skewed, log(x + 1) transformation was done prior to further analysis. Spatial autocorrelation analysis was conducted using Moran’s I statistic followed by optimized hotspot analysis. Moran’s I, inferred using z-score and p value, determines the presence of spatially significant clusters. If z-score value is greater than 2, it indicates highly significant (p < 0.01) clustering; z < −2 indicates significant dispersed (p < 0.01) and z around 0 indicates random distribution (Yavuzaslanoglu et al., 2012). Spatial relationship was conceptualized using fixed distance band where each feature is analyzed within the context of neighboring features. If global spatial autocorrelation is detected, optimized hotspot analysis tool is used to create statistically significant hot and cold spots using Getis-Ord Gi* local spatial statistic (Ord and Getis, 1995). For a given set of weighted data, Getis-Ord Gi* statistics looks at each feature within the context of its neighboring features and it can be mathematically explained by the following equations: (2) G i * = ∑ j = 1 n w i j x j − X ¯ ∑ j = 1 n w i j n ∑ j = 1 n w i j 2 − ( ∑ j = 1 n x w i j ) n − 1 S , \text{(2)} & {G}_{i}^{* }=\frac{\sum _{j=1}^{n}{w}_{ij}{x}_{j}-\overline{X\,}\sum _{j=1}^{n}{w}_{ij}}{\sqrt[S]{\frac{n{\sum }_{j=1}^{n}{w}_{ij}^{2}-({\sum }_{j=1}^{n}x{w}_{ij})}{n-1}}},\end{array}\end{eqnarray}?> (2)where x _j is the attribute value for feature j. w _ij is the spatial weight between feature i and j, n is equal to the total number features and: (3) (3) X ¯ = ∑ j = 1 n x j n , \text{(3)} & \overline{X\,}=\frac{{\sum }_{j=1}^{n}{x}_{j}}{n},\end{array}\end{eqnarray}?> (3) (4) (4) S = ∑ j = 1 n x j 2 n − ( X ) ¯ 2 . \text{(4)} & S=\,\sqrt{\frac{{\sum }_{j=1}^{n}{x}_{j}^{2}}{n}}-({\overline{X)}}^{2}.\end{array}\end{eqnarray}?> (4)

G_i* is the z-score, which represents statistical significance of clustering for a particular distance. Positive z value signifies high values of clusters similarly negative z represents clusters of low value, while z = 0 represents no statistical significance. This tool also works by looking at each feature within the context of its neighboring feature. So, in a particular location if nematode population density is high and neighboring sites also possesses high nematode density, then that particular location is considered to be part of a ‘hotspot’ and vice versa (for ‘cold spots’). After optimized hotspot analysis, IDW (spatial analyst) tool was used to create an interpolated surface based on z score which effectively visualize the results of hotspot analysis.

Spatial interpolation is a tool to predict the values of a spatial phenomenon at unsampled locations, like if population density of RWTN at locations (x ₁, x ₂, …, x _n ) are (z ₁, z ₂, …, z _n ), the purpose of spatial interpolation is to estimate the value z at some new point of x. Inverse distance weighting (IDW) and ordinary kriging (OK) are two important tools to generate estimated surface across the region of interest.

IDW is one of the simplest interpolation techniques where weighted mean of the neighboring observations are taken into consideration. Weights are usually inversely proportional to a power of distance (Burrough, 1986; Watson, 1992) which, at an unsampled location r, it estimates can be expressed as follows: (5) (5) F ( r ) = ∑ i = 1 m w i z ( r i ) = ∑ i = 1 m z ( r i ) / | r − r i | p ∑ i = 1 m 1 / | r − r j | p , \text{(5)} & F(r)=\,\sum _{i=1}^{m}{w}_{i}z({r}_{i})=\frac{{\sum }_{i=1}^{m}z({r}_{i})/\,|r-{r}_{i}{|}^{p}}{{\sum }_{i=1}^{m}1/\,|r-{r}_{j}{|}^{p}},\end{array}\end{eqnarray}?> (5)where p is a parameter and m is the number of neighbors taken into consideration at a certain cut-off distance. IDW (spatial analyst) tool was used for this purpose with variable search radius and surface map was generated based on z value.

From theoretical point of view, kriging is the optimal interpolation technique to estimate a random variable z at one or more unsampled locations (Santra et al., 2008). In kriging, experimental variogram is used to compute spatial correlation of the random function z(x ₀ ) and can be defined by calculating the semivariance as a function of distance (Eq. (6)): (6) (6) γ ( d ) = 1 2 ∑ { [ z ˆ ( x 1 ) − z ˆ ( x 2 ) ] 2 } . \text{(6)} & \gamma (d)=\,\frac{1}{2}{\sum }^{}\{{[\hat{z}({x}_{1})-\hat{z}({x}_{2})]}^{2}\}.\end{array}\end{eqnarray}?> (6)

Surface map of population density of RWTN was prepared using ordinary kriging (OK) tool in ArcMap 10.2. As population density data showed strong positive skewness, log(x + 1) transformation was carried out before further analysis. Ordinary kriging can be mathematically explained as follows: (7) (7) z ˆ ( x 0 ) = ∑ i = 1 n λ i z ( x i ) , \text{(7)} & \hat{z}({x}_{0})=\,\sum _{i=1}^{n}{\lambda }_{i}z\,({x}_{i}),\end{array}\end{eqnarray}?> (7)where z is the variable of interest at location x _i and x ₀; n is the number of neighbors taken into consideration; λ _i are the weights and z(x _i ) is the value of z at x _i (Armstrong and Boufassa, 1988; Stein, 2012). Before performing kriging, population density data were explored using histogram and normal QQ plot tool in geostatistical analyst, ArcMap 10.2. Thereafter, trend analysis and semivariogram/covariance cloud were performed to identify the presence of global trends and directional influence within the data set. If any trends were observed then trend removal function was used and if the data set was found to have directional influence, then anisotropy was taken into account for interpolation model. To find the appropriate lag size, average nearest neighbor under spatial statistics tool was used. It calculates a nearest neighbor index based on the average distance of each feature from its neighboring feature. Spatial variability of population density of RWTN was expressed by a semivariogram (Warrick et al., 1986; Goovaerts, 1998) which computes the average dissimilarity between data separated by a vector h and can be explained as follows: (8) (8) γ ˆ ( h ) = 1 2 N ( h ) ∑ i = 1 N ( h ) [ z ( x i ) − z ( x i + h ) ] 2 , \text{(8)} & \hat{\gamma }(h)=\,\frac{1}{2N(h)}\sum _{i=1}^{N(h)}{[z({x}_{i})-z({x}_{i}+h)]}^{2},\end{array}\end{eqnarray}?> (8)where N(h) is the number of data pairs within a given class of distance and direction; z(x _i ) is the value of variable at location x _i , z (x _i  + h) is the value of variable at a lag of h from x _i . Semivariogram values were fitted with different theoretical models like exponential, circular, Gaussian, and hole effect model. Mathematical expression of these semivariogram models are as follows.

Spherical model: (9) (9) γ ˆ ( h ) = C 0 + C [ 1.5 h a − 0.5 ( h a ) 3 ] , if 0 ⩽ h ⩽ a . \text{(9)} & \hat{\gamma }(h)={C}_{0}+C[1.5\frac{h}{a}-0.5{(\frac{h}{a})}^{3}],\,{\rm{if}}\,0\leqslant h\leqslant a.\end{array}\end{eqnarray}?> (9)

Exponential model: (10) (10) γ ˆ ( h ) = C 0 + C [ 1 − exp { − h a } ] for h ⩾ 0. \text{(10)} & \hat{\gamma }(h)=\,{C}_{0}+C[1-\exp \{-\frac{h}{a}\}]\,{\rm{for}}\,h\geqslant 0.\end{array}\end{eqnarray}?> (10)

Gaussian model: (11) (11) γ ˆ ( h ) = C 0 + C [ 1 − exp { − h 2 a 2 } ] for h ⩾ 0. \text{(11)} & \hat{\gamma }(h)={C}_{0}+C[1-\exp \{-\frac{{h}^{2}}{{a}^{2}}\}]\,{\rm{for}}\,h\geqslant 0.\end{array}\end{eqnarray}?> (11)

Hole effect model: (12) (12) γ ˆ ( h ) = C 0 + C [ 1 − sin ( π h / a ) π h / a ] for h ⩾ 0. \text{(12)} & \hat{\gamma }(h)={C}_{0}+C[1-\frac{\sin (\pi h/a)}{\pi h/a}]\,{\rm{for}}\,h\geqslant 0.\end{array}\end{eqnarray}?> (12)

C _o, (C + C _o) and a in the above semivariogram models are nugget, sill, and range, respectively; in case of exponential and Gaussian models, a represents the theoretical range. These are the parameter values for omnidirectional semivariogram model. Cross-validation was performed based on three considerations; mean prediction error close to zero; root mean square standardized prediction error close to 1, and root mean square (RMS) and average standard error (ASE) as small as possible. Finally, population density krigged map of RWTN was prepared and symbolized properly.

To identify the areas where population density of RWTN is above the economic threshold level (ETL) level of 300 nematodes/100 seeds, indicator kriging was performed. Indicator kriging was chosen for this purpose because a map showing areas where the population density of RWTN exceeds the ETL would greatly benefit the rice growers. Indicator kriging was performed following steps similar to that mentioned in case of ordinary kriging. Color coded krigged map was generated with contour symbolization delineating high risk areas.

The nematode species was confirmed to be Aphelenchoides besseyi. Infestation of A. besseyi was found throughout the district with varying densities with no significant variation across the blocks (χ ² = 12.9, p > 0.05) (Fig. 2B). Highest population density of the nematode was recorded in the Dumri block followed by Giridih and Bengabad. Lowest RWTN infestation was observed in Deori and Gawan blocks. The agglomerative hierarchical clustering of population density among the 12 blocks of Giridih district identified three main clusters using average linkage method (Fig. 2B). Giridih and Bengabad blocks formed a cluster, while the second cluster was comprised of four blocks, namely, Pirtanr, Gawan, Bagodar, and Tisri. The third cluster was consisted of six blocks with a clear separation between the Dumri block and other members of the cluster due to high population density of RWTN. Silhouette width (S _i ) was found to be 0.19 which represents good quality index for number of clusters found. Average S _i values for three clusters are 0.27, 0.15, and 0.18, respectively (Appendix Fig. A1).

Figure 2:

Moran’s I spatial autocorrelation results suggest presence of significant (z-score 2.36, p < 0.05) spatial clustering pattern among population density of RWTN. Distance threshold of Moran’s I was found to be 9,982.67 meter. Red and green circles show points on the map (Fig. 3) represent hotspots (Gi Bin > 2) and cold spots (Gi Bin < −1), respectively, based on Gi Bin values. Getis-Ord Gi* statistic also suggest presence of significant (z-score > 2) hotspot in the south western parts of Giridih districts comprising of adjoining regions of Dumri and Bagodar blocks (Fig. 3). Significant cold spots (z-score <−1.5) were observed in the adjoining areas of Gawan, Tisri, and Dhanwar blocks (Fig. 3).

Figure 3:

Indicator kriging was performed similarly like ordinary kriging. Among different experimental models (see Supplementary Fig. A4), here also exponential model was found to be best fitted model. Lag size of 3037 meter was used as in case of ordinary kriging. Semivariogram of exponential model has been shown in Figure 5B. Model parameters and cross-validation results were presented in Table 1. Probability distribution map (Fig. 7) was prepared with the threshold value of 300 nematodes/100 grains. Of total interpolated surface, 58.5% area belongs to low probability (20%) risk areas like Gawan, Deori, Birni, Sariya, Bagodar, and parts of other blocks, whereas up to 50% (medium) and 70% (high) probability to cross the ETL comprised of 29.75% and 7.49% of total interpolated surface. Very high probability (75-100%) to cross the ETL boundary was found in Dumri, Giridih, Gande, Bengabad, and some other fragmented parts of the district which comprised 4.17% of the total interpolated surface.

Figure 7:

Rice is the predominant crop of the Giridih district where most of the fields remain fallow after the harvest of rice till next season. The shattered rice grains and regrowth of rice stubbles as ratoon rice maintain and carry over the nematode population from season to season. Furthermore, RWTN present in rice stubbles in the field after rice harvest and the fungi (Fusarium sp., Curvularia sp.) found in the straw were considered probable alternate hosts for the nematode and helped the inoculum to carry over in the rice field (Sivakumar, 1987). The nematodes live in field conditions along with infected seeds might act as the initial inoculum of RWTN for the new rice season. The incidence of RWTN in kharif and boro rice (summer rice) has been recorded from the different districts of West Bengal, India (Pashi et al., 2017a). Moreover, ordinary kriging has been used in the analysis of spatio-temporal distribution of different insects like Helicoverpa armigera (García, 2006). Apart from this, integrated pest management approaches have been designed against pest of orchard and vineyards using geostatistical approaches like ordinary kriging (Sciarretta and Trematerra, 2014). So, spatial interpolation tools like IDW and OK might be useful techniques to investigate the distribution pattern of RWTN across the district.

Being a seed-borne pathogen, presence of significant spatial clustering of RWTN infested fields, as shown by the Moran’s I spatial autocorrelation, indicates that exchanges of infected seeds and plant materials like rice husk and plant debris or infected straw might have resulted in the spread of the nematode to the nearby fields from the initial infestation foci (Sivakumar, 1987; Bridge et al., 2005). Infestation hotspots in the Dumri and adjacent blocks could act as sources of future spread of the nematode disease in other areas of the Jharkhand State. Repeated use of infested seed material could lead to high population build-up of nematodes as observed in these locations. Such a scenario of spread of nematode infestation has been referred to as the contagion effect scenario and was previous observed in case of spread of the golden nematode of potato (Globodera pallida) in the uninfected potato fields of Idaho (Contina et al., 2018). They showed that spread of G. pallida grew in diameter from the original center of infestation toward south west of Idaho might be because of contaminated soil with agricultural equipment or tubers. Several examples exist in literature where both point pattern and surface analysis were utilized to determine nematode infestation pattern. Point pattern or hotspot analysis helps identify the population hotspots among the sampling points, in contrast, surface interpolation helps to understand the probable population density in unsampled locations by generating smoothed surface maps. Presence of significant RWTN infestation hotspot and cold spot in Dumri and adjacent blocks of Gawan, Tisri, and Dhanwar, respectively, might help in future to take management decisions. In another example, point pattern geostatistical approach was employed to identify infestation hotspots of the cereal cyst nematode (Heterodera filipjevi) and root lesion nematode (Pratylenchus neglectus, P. thornei) infesting wheat and barley in Turkey (Yavuzaslanoglu et al., 2012). They also found that population density of these nematodes were positively correlated with electrical conductivity (EC) and percentage sand content of soil. While the primary mode of RWTN infestation is via spread of infected seeds, however, like other phytoparasitic nematodes, environmental factors could also influence population density of RWTN both spatially and temporally (e.g., see Huang et al., 1972; Sun et al., 2009). In this study, we have used the population density information generated during our surveys to generate the spatial distribution pattern of RWTN in Giridih district. As discussed earlier, the population density of the nematode could change over time resulting in changes in the spatial patterns of distribution as well. Therefore, repeated survey of RWTN for several years might help to establish the trend in population density of this important nematode.

Both the surface interpolation techniques, IDW and ordinary kriging, yielded similar prediction surfaces of nematode population density across the district. Experimental variogram for population density of RWTN showed relatively strong spatial correlation indicating the presence of spatial dependency. As spatial dependency is present, kriging is considered as better interpolation tool than IDW due to problems associated with distance-based interpolation methods. The use of ordinary kriging generated map is therefore recommended (Yao et al., 2013; Gong et al., 2014; Shahbeik et al., 2014). In Brazil, a Rotylenchulus reniformis infested cotton monoculture field was sampled to develop a risk map for nematicide application for small crop area using ordinary kriging (Farias et al., 2002).

One of the major drawbacks of using ordinary kriging is that the smoothed map it generates does not capture the extreme population density values (Farias et al., 2002). To avoid this issue, indicator kriging approach was implemented to generate probability risk map of RWTN infestation. This approach corroborates with the study of southern root knot nematode (Meloidogyne incognita) infestation in cotton where risk maps were prepared based on population density, correlation between electrical conductivity of soil and population density or a combination of both using indicator kriging (Ortiz et al., 2010). The risk map generated in our study identifies areas with high probabilities were RWTN population density will cross the economic threshold limit of 300 nematodes/100 grains. As discussed in the introduction section that there is an ambiguity regarding the economic threshold level of A. besseyi. We have used it as 300 nematodes/100 grains, if ETL is <300, the analysis of indicator kriging will be changed accordingly.

The spatial distribution maps generated in our study could be utilized by the farmers and extension agents in devising precautionary measures and formulating management strategies to prevent further spread of this important plant parasitic nematodes in subsequent growing seasons in the district. Management of RWTN could consist of two major strategies: curative and preventative. Disinfection of seed could be a possible means of complete elimination of white tip nematode disease. In high risk areas, curative control measures like sun drying of seeds and pre-soaking seeds in cold water for 18 to 24 hr followed by hot water treatment (51-53°C) for 15 min prior to sowing can be effectively used to minimize the population density build up (Kuriyan, 1995; Bridge et al., 2005; Pashi et al., 2017b). In addition, agronomic practices like early seed bed preparation, irrigation of seed beds immediately after seeding to asynchronies nematode hatching with rice germination, cultivating leguminous crops after rice and destroying crop residues after harvesting could help in decreasing field populations of RWTN (Yamada et al., 1953; Kim et al., 1996).

Like other seed-borne nematodes, establishment of quarantine measures to prevent further spread of RWTN infected seed material from high risk areas could be an effective strategy in stopping the spread of this nematode in low risk areas. This will not only help in reducing the population density load in high risk areas subsequent years, but also reduce the dissemination of the disease in the nearby areas. Use of resistant variety is another eco-friendly management strategy against nematodes. For example, in Russia, resistance screening was conducted with 1,003 rice cultivar against RWTN in glasshouse and three were found to be immune, 10 highly resistant, 164 moderately resistant and rest were susceptible (Popova and Subbotin, 1994; Tülek et al., 2015). However, no studies were conducted so far in India to identify rice cultivars resistant against RWTN, to our knowledge.

Lack of knowledge among farmers regarding this nematode is a crucial factor for wide spread distribution of plant parasitic nematodes in general and RWTN in particular. As nematodes are microscopic and generally do not generate any specific symptoms, farmers often confuse nematode infestation with nutrient or water deficiency. To improve farmers’ perception about RWTN infestation, farmers training workshops can be conducted in agricultural extension centers. Increasing awareness among farmers will reduce the use of contaminated seeds which will help in prevention of spread of RWTN in the district.

In future, the geostatistical methods implemented in our study could be used to generate spatial distribution maps of other agriculturally important plant parasitic nematodes. Beside rice, Aphelenchoides besseyi is also considered as a serious pest of tuberose (Polianthes tuberosa L.) in West Bengal, India causing the ‘floral malady’ disease which inflicts huge economic loss to the farmers (Pathak and Khan, 2009). Similar to our study, identification of A. besseyi infestation hotpots in tuberose infesting areas of West Bengal will be highly beneficial to the farmers by preventing its spread in non-infested areas.