Spatiotemporal analysis of distribution of pest and predator in corn crops

Crop pests have negative impacts on yield. This paper proposes a spatiotemporal geostatistical modeling to compare data of adult Syrphidae fly count and corn leaf aphid Rhopalosiphum maidis (Fitch, 1856) colonies in corn crops. The use of a geostatistical model that allows the space-time variation makes the approach more interesting because it is a more complete model. Multiple regression was used to model the trend component for the variable response adult Syrphidae fly count and corn leaf aphid colonies, with the coordinates serving as covariates and the spatiotemporal variations around the deviation are described by a random spacetime residual field. Finally, the prediction map obtained by kriging may be a biological indicator of possible corn leaf aphid colonies in the corn crop. It was possible to verify that the occurrence of the pest provided a significant increase in adult predators and seminatural habitats may favor populations of natural enemies.


INTRODUCTION
Invasive species are among the most important impacts on humanity, less controlled and less reversible in the ecosystems of the world, with negative consequences, affecting its sustainability, biodiversity and economic systems (SILVA et al., 2017). In the maize crop, for example, the spike attack and stalk can lead to grain yield losses of 15 to 34% (LIMA et al., 2010).
These pests are usually controlled with the use of agrochemicals, however some pests like the corn leaf aphid Rhopalosiphum maidis (Fitch, 1856) present natural enemies that decrease their population. This can be used as a strategy in integrated pest management. BORTOLOTTO et al. (2016) report that the occurrence of populations of Syrphidae flies accompanies the aphid colonies, increasing or decreasing proportionally. This is because the flies from the Syrphidae family are adapted to lay their eggs near the corn leaf aphid colonies, since they feed on aphids in the larval phase and on pollen and nectar in the adult phase. Therefore, the presence of adult Syrphidae flies in the field may be a probable biological indicator of corn leaf aphid colony. PARK;OBRYCHI (2004) say that one way to investigate the spatiotemporal synchronization of predator and prey distributions is to generate and compare distribution maps in temporal sequences.
Tools that can describe the behavior of predators are very useful for indicating and monitoring infestations and pests. These tools include sampling and spatial analysis. Pest sampling works can be seen in FARIAS et al. (2001a). The use of spatial models is already commonly used in the literature, as can be seen in FARIAS et al. (2001b), DUARTE et al. (2015), GARCIA et al. (2017), MEHRJARDI et al. (2008), among others. A more complete approach is time-space analysis, which uses robust models to describe and predict pest incidence by studying space-time dynamics together. M.N. Mello et al. Spacetime is seen as a continuous spatial arrangement combined with a temporal order of events. This union of space and time is defined in terms of its Cartesian product (CHRISTAKOS; VYAS, 1998). This extension of spatial models is a logical evolution of statistical models for mapping spatially and temporally correlated data, and facilitates the distinction between purely temporal, purely spatial or spatiotemporal variability components (KILIBARDA et al., 2014).
The objective of this paper is to model the spatiotemporal variability of corn aphid and Syrphidae family colonies in corn crops to investigate and monitor the dynamics of prey and predator infestation at different phenological stages, serving as a possible biological indicator in the field from time maps.

Study area and sampling
The study was conducted at the Experimental Farm of Igarapé-Açu, FEIGA (01°07'24"S and 47°37'22"W), located in the municipality of Igarapé-Açu, in the Bragantina microregion, in the mesoregion of northeastern Pará, Brazil. The site has an altitude of 39 m above sea level. This farm is recognized regionally for developing studies to improve animal and agricultural production techniques. The farm belongs to the Universidade Federal Rural da Amazônia, but it is also used by Embrapa and Adepará. The climate is defined as Ami, with an average annual rainfall of 2500 mm and an average annual temperature of 25 °C, according to the Köppen classification.
The experimental area covers 10,000 m (1.0 ha). A mango agroecosystem (1.7 ha), a pasture area with a predominance of Megathyrsus maximus (Jacq.) B.K. Simon and S.W.L. Jacobs (0.6 ha), and fragments of secondary forest are adjacent to the main site. The area was divided into 100 plots of 100 m (10 × 10 m). Typical planning of studies in spatial and temporal distribution must take into account spatial grid sampling in the same area and within a regular time interval (HENDERSON;SOUTHWOOD, 2016).
Prior to the corn crop, the area was cultivated with cowpea beans (Vigna unguiculata (L.) Walp), and, after harvesting, it was kept in a fallow period. On April 4, hybrid corn Priorizi M274 Morumbi was sown, with an early cycle that presents rusticity and tolerance to the main corn diseases. Both were sown with a spacing of 0.90 m between rows and 0.15 m between plants, without chemical control intervention during the cultivation period.
Sampling occurred on April 16 (12 days after sowing) to June 25, 2016, totaling 11 observations. The samples that provided the highest incidence of leafhoppers were analyzed using spatiotemporal methodology.
Ten plants were randomly chosen per plot, resulting in a total of 1,000 plants in each sampling. In order to detect the leafhoppers, a visual analysis of all parts of the plants (stem, cartridge, leaves, and ears) was performed. Subsequently, the number of insects per plot was counted. Some sharpshooters were captured with an entomological aspirator and sent to the Biodiversity Laboratory of the Universidade Federal Rural da Amazônia (UFRA), Capanema campus, to be morphotyped. Subsequently, they were sent for identification at the Entomology Department of the National Museum, Universidade Federal do Rio de Janeiro (UFRJ).

Spatiotemporal models
It is a Gaussian spatiotemporal process Z defined over a spatial domain S and a temporal domain T contained in the set of real numbers R, {Z(s,t): (s,t) ∈ (S × T)}, where S ⊆ R d , d = 2 represented the two dimensions in space (longitude and latitude) and T ⊆ R represented only one dimension, namely time. A statistical spatiotemporal process, observations are modeled as a partial realization of a random spacetime function (XU; SHU, 2015). The spatiotemporal variation of Z can be decomposed by the components of the trend m(s,t) and a stochastic residual ε(s,t) (YANG et al., 2015), as can be seen in Equation 1.
In Equation 1 was assumed that Z has the moments from the first to second order. The mean component m(s,t) = E[Z(s,t)] which is the expectation of variable Z. Stochastic residual incorporates three components: spatial, temporal and interaction (DE IACO et al., 2015) and it is given by Equation 2.
where x 1 , x 2 , x 3 and x 4 are the predictor variables of the regression model and β 0 , β 1 , β 2 , β 3 and β 4 are the parameters to be estimated.
For modeling, it was declared that these components are second-order stationary, mutually independent and spatially isotropic. As a binding function, it is g[m(s,t)] = log [m(s,t)]. Due to the high numbers of zero in the response variable, the trend component was modeled by a zero-inflated Poisson (ZIP) distribution. Were considered as covariates the geographic coordinates (latitude and longitude) and a linear and quadratic temporal index. The adjusted trend model is given by Equation 3. " ( , ) = + ! + + " " + + # # + + $ $ + + ! !
where the variables x 1 , x 2 represent latitude and longitude, respectively. The variables x 3 and x 4 represent the linear and quadratic temporal effects, respectively. After the trend modeling, the next step was to model the spatiotemporal variogram on the residuals (MARTÍNEZ et al., 2017). This spatiotemporal variation not explained by the covariates is modeled by the residual component ε(s,t) = Z(s,t) -m(s,t) using the spatiotemporal variogram (γ s,t ) given by Equation 4.
where hs = s -s' and h t = t -t' is for both (s,t) and (s',t') in the space-time domain. The next step in the analysis is to determine a theoretical model based on covariance functions that fit the data of the spatiotemporal sample variogram. To fit a theoretical model to the empirical variogram presented in Equation 4, the generalized product-sum model was used (GRÄLER et al., 2016). Variogram, half the difference in variance, is generally more useful than the covariance function because of its weaker assumptions (XU; SHU, 2015). Thus, the relationship of covariance and semivariogram in second-order stationary processes (constant mean and stationary covariance function) is given by Equation 5. (ℎ ! , ℎ " ) = (ℎ ! , 0) + (0, ℎ " ) − (ℎ ! , 0) (0, ℎ " ) where γ st (h s ,0) and γ st (0,h t ) are valid spatial and temporal marginal variograms, respectively. Parameter k in Equation 5 is related to spatial and temporal sill.

Spatiotemporal kriging
For HEUVELINK et al. (2012) the spatiotemporal kriging formulas do not differ fundamentally, in a mathematical or statistical sense, from those of spatial kriging. According to KILIBARDA et al. (2014), the variogram model is crucial in spatiotemporal kriging to calculate the best unbiased linear predictor, which is given by Equation 7.
where (s 0 ,t 0 ) is the linear predictor, obtained by kriging, of the location (s 0 ,t 0 ), C n is the n × n covariance matrix of the order n × n of the residuals for the n points observed in spacetime, derived from the spatiotemporal variogram, c 0 is a vector of residual covariance for the observed and predicted points and ε is a vector of the residuals at the n observed points.
where m(s 0 ,t 0 ) is the estimated value for the location (s 0 ,t 0 ) obtained from Equation 3 and (s 0 ,t 0 ) as defined in Equation 7. To determine the best theoretical variogram model that fits the sample variogram, the root mean square error (RMSE) was calculated, which was obtained by taking the difference in the point clouds between the sampling variogram and the theoretical variogram.

RESULTS
A gradual increase in the number of corn leaf aphid colonies can be observed over the weeks, with its apex in the last sampling (Fig. 1 a). In Figure 1 Figure 2 present descriptive measures and boxplots as exploratory graphs of (a) adult Syrphidae fly count and (b) corn leaf aphid colony (> 15 individuals) by weekly sampling intervals on a hectare located in the northeastern of the state of Pará. It is possible to verify the high variability of aphid colonies and an approximation of the average and standard deviation in adult Syrphidae flies. It is possible to observe that there are some periods with no observation of species. In the graphs, it is possible to see the great existence of zeros in both species and a possible high variability in each sample. This provides the idea that models that can handle such phenomena should be used to model the trend, such as the ZIP.  To model data trend the regression model was used as indicated in Equation.
3. Table 2 shows the parameter estimates, standard error, t-statistics, and the p-value for corn leaf aphid colonies and adult Syrphidae flies. It is possible to verify that at 0.05 significance level, covariates latitude and time are associated with the counting of corn leaf aphid colony and latitude, longitude and the linear quadratic temporal effects are associated with the adult Syrphidae flies count. The residual produced by this model will be used for spatial dependence analysis. After removing the spatiotemporal trend of the adult Syrphidae flies count and corn leaf aphid colony count, the variogram on the residual is calculated. Figure 3 shows the empirical spatiotemporal variogram model (a) and (c) and the adjusted product-sum model (b) and (d). As can be seen, the residues have a clear correlation in both space and time and the total variation of these residues is explained by spatial as well as temporal components. Note that the spatial structure becomes weaker as time differences increase and the temporal structure becomes low as spatial differences increase. The sample variogram (a) and (c) can be viewed in terms of their marginal variograms, one for space and one for time. The increasing trend in spatial and temporal dimensions in the sample variogram indicates the presence of a strong spatiotemporal correlation in both. Table 3 presents the parameter estimates of the generalized product-sum model considering the different structures for the spatial and temporal components of adult Syrphidae flies counts and corn leaf aphid colony counts. In order to compare the models, the RMSE was used. Works such as by DE IACO et al. (2015),    In this paper, a broader analysis using different structures in the formation of the model was proposed. For the generalized product-sum variogram model, the exponential model for the spatial component and the linear model for the temporal component for adult Syrphidae and Gaussian and linear for corn leaf aphid were considered. The spatial range of adult Syrphidae was 46 m, suggesting that the spatial correlation became insignificant after this threshold. On the other hand, the time interval of 23 days indicated that the temporal correlation became insignificant after this interval. For the aphid colony, the spatial amplitude was 45 m and the temporal interval was 20 days.
Figure 4 (a) shows an increase in corn leaf aphid colonies weekly. In the samplings, the corn leaf aphid emergence in the corn tearing stage until the last day of sampling with the tooth formation stage. In Figure 4 (b), it is possible to observe that there is a tendency of colonization that starts from the west to east. In the first sampling (May 6, 2016), there are a very small number of flies. From the second sampling (May 14, 2016) is already clear the increase of flies in the plantation. This is possibly due to the food supply that must be growing (corn leaf aphid colony). The sampling of May 28, 2016 presents a high rate of infestation. This phase is the corn tearing phase and, according to WEISMANN (2008), it is at this stage that the plant is most susceptible to attacks due to the greater exposure of the tassel and all the leaves. In the last two samples, it is already possible to verify a small decrease in the fly infestation. The errors were estimated and the highest values were at the edge of the area.

DISCUSSION
Natural aphid infestations can produce a range of population densities through different stages of plant growth and development. AL-ERYAN; EL-TABBAKH (2004) report this fact in the stages of growth V10, tearing (VT) and ripening/maturation. They also report that infestation by aphid induces yield loss of 28.14% at growth stages V10-VT, while infestation at maturity growth stages R2 to R4 caused yield loss of 16.28%. CRUZ (2004) comments that infestation begins in isolated plants, mainly in the vegetative period, near the VT stage. This study showed the emergence of the aphid at the VT stage, confirming the hypothesis and higher occurrence at the reproductive stage of tooth formation (R5).
Studies such as YANG et al. (2019) showed that uncultivated habitats, such as forest plots and vegetation around dwellings or wetlands in the landscape, all had positive correlations with corn leaf aphid abundance in wheat fields because they provided sources for corn leaf aphid colonization. Besides the experiment there was an area of hose, pasture and forest fragment. These were shown to influence the colonization of corn leaf aphids and adult Syrphidae flies also in maize, as shown in the maps.
Research such as DAMICONE et al. (2007) and KLINGAUF (1987) report that wind may influence the distribution of other types of winged corn leaf aphids. The corn leaf aphids emerged near the pasture area, and wind can be a possible factor for its displacement. As corn leaf aphid colonies began to appear and grow, there was a significant increase in the number of adult Syrphidae flies in the area (May 28, 2016) that emerged from adjacent areas of forest fragments and entered the corn crop.
Probably, a forest fragment served as a shelter for the flies and provided food, since in adulthood flies feed on pollen and nectar (HOLLOWAY, 1976;MÜLLER, 1883). These seminatural habitats favor populations of natural enemies and improve their efficiency as control agents (ALIGNIER et al., 2014).
Adult flies lay their eggs near aphid colonies (WHITE et al., 1995) for the larvae to feed on them (BELLIURE; MICHAUD, 2001); so, observing adult flies is a possible indication of colonies. Although we do not have samples at the end periods of flies, it is possible to infer that the decrease in flies led to the increase in colonies or that the occurrence of this predator alone was not sufficient to decrease the aphid population. A likely natural barrier to colony control in to experiment in adjacent areas using the flowers to contribute to the fast colonization of adult flies in the presence of aphids.
Many studies, such as GASCH et al. (2015) and MEDEIROS et al. (2019) among others, show that inserting more covariates can improve the accuracy of the model. Perhaps the insertion of one or more covariates could improve the prediction in this study and better explain the influence of the pest versus predator relationship. Some advantages of using this approach are that it provides flexibility for estimating the trend component and it can consider both continuous and categorical variables, as well as applying different binding functions on the predictor, not requiring a Box-Cox transformation, as in HUSSAIN et al. (2010), besides modeling space and time together, not always used in applied research. Some examples of marginal analyses of these two components can be verified in papers from PELISSARI et al. (2017), ROJO; PÉREZ-BADIA (2015) and SCIARRETTA et al. (2018). A good fit can provide Gaussian errors that can be easily interpolated by keeping kriging assumptions and diagnostic tools available to evaluate the model fit (POGGIO et al., 2012). However, the use of this approach has the disadvantage that, in order to model the trend, it assumes independence and the data has known spatiotemporal dependence.

CONCLUSION
A current and flexible approach to estimating the spatiotemporal distribution of pests and natural enemies with a ZIP model and weekly temporal sampling was presented in this paper. The relationship of pests and predators may be a biological indicator in crop control to mitigate the use of agrochemicals in production crops. Seminatural habitats may favor populations of natural enemies. It was possible to verify that the occurrence of the pest provided a significant increase in adult predators. Finally, it is possible to infer that natural barriers can be used to control aphid colonies, favoring their natural enemies.