Potential Bioclimatic Ranges of Crop Pests Zabrus tenebrioides and Harpalus ruﬁpes during Climate Change Conditions

: The ground beetles Zabrus tenebrioides and Harpalus ruﬁpes (Coleoptera, Carabidae) are two of the most prevalent pests of wheat and other grasses. This article presents current data on their distribution and the results of modelling the bioclimatic ranges of these species using the maximum entropy method. To improve the model, we used various RStudio packages including the R script “thin points 4-1-18.R” package spThin and the «Raster» RStudio package. We determined the climatic parameters that promote the dispersal of the species, as well as the optimum conditions for the growth of Z. tenebrioides and H. ruﬁpes . Maps forecasting the distribution of the studied species were generated through the perspective of two climate scenarios: RCP 2.6 and RCP 8.5. For the modelling, we utilised 435 geographic points of Z. tenebrioides occurrence and 653 points of H. ruﬁpes occurrence. Both species have similar bioclimatic ranges, and the most favourable conditions for them are ﬁelds of grain crops. The most signiﬁcant parameters inﬂuencing Z. tenebrioides are those of moisture, whereas H. ruﬁpes is most sensitive to the temperature parameters. According to the generated climatic models for both species, a decrease in the areas of their ranges would occur in their eastern, more continental areas, with a slight shift towards the north.


Introduction
Ground beetles (Carabidae) are one of the most important components of the fauna of natural and anthropogenic ecosystems. The practical significance of ground beetles is first determined by their diverse biology and their abundance in all terrestrial communities [1][2][3]. In most cases, their significance is positive because of the primary role of predatory Carabidae in regulating the numbers of many insects, terrestrial mollusks and other invertebrates, including a number of harmful pests of agriculture and forestry [4][5][6][7][8][9]. At the same time, some species of ground beetles (phytophages and mixophytophages) notably harm agricultural crops, and sometimes harm pastures and planted seeds of forest species [10][11][12][13][14]. A great level of practical importance is acknowledged for representatives of the Zabrini and Harpalini tribes [15][16][17][18]. They mostly belong to mixophytophages that feed on plants and animals and, less commonly, on specialised phytophages [19].
Identifying the patterns through which ranges of ground beetles establish themselves is a fundamental issue closely associated with the history of species in the context of the global and local evolution of their natural environment and biodiversity [20,21]. Having determined these patterns, we are able to estimate the perspectives of the further existence Diversity 2021, 13, 559 2 of 15 of species and changes in their ranges caused by various factors. Such studies are closely related to solving numerous practical tasks of fauna genesis, biogeography, their protection and a rational use of natural resources [22]. The emergence of new methods, particularly of GIS-technologies, has significantly extended the opportunities for ecologists, allowing them to simulate the dynamics of the ranges of species based on their relationships with environmental factors [23].
The present-day ranges of many species and the changes they undergo significantly depend on the climate. Therefore, identifying the climatic niche of species is an important component of biogeographic and ecological studies [24]. The study of climate factors is especially interesting because of the current global changes in the climatic system and the rearrangement of biospheric processes [25,26]. Methods for ecological modelling, based on GIS technologies, allow for the distribution of species to be modelled not only for the current period, but also for the future by taking into account forecasts of climate changes. Furthermore, the distribution of a species depends largely on the trophic needs at the various stages of its development, the details of the spatial distribution of any competing species, the predatory and parasitic species and technogenic and agrogenic pollutants in the environment [27][28][29][30].
Zabrus tenebrioides (Goeze, 1777) is an average-sized beetle with a broad distribution range spanning the area from England and southern Sweden to North Africa and Asia Minor [31][32][33][34]. This species is common in Europe, found in almost all southern regions where grain crops are cultivated, and consumes wild-growing and cultivated grasses [35][36][37]. Zabrus tenebrioides is a dangerous pest of grain crops and may cause great damage to barley, wheat, rye and maize. Both beetles and their larvae are harmful [38]. The long-lasting evolution of the representatives of the Zabrus genus allowed this species to almost completely abandon predation [39][40][41][42][43][44][45]. A much lower adaptivity to a plant-based diet is observed among representatives of the Harpalini tribe [46].
Global climate changes may affect the ranges of Z. tenebrioides and H. rufipes. According to a number of studies, including studies we conducted earlier, the climate changes cause a shift in climatic ranges towards the North. Nonetheless, different species of ground beetles respond to these changes at different rates. For certain species, the impact of global warming is relatively low [68][69][70], however, mesophilic species living under forest canopies would face the greatest impact due to changes in the range of forest communities [71][72][73][74]; whereas hygrophilic [75][76][77][78][79] and solonchak [80][81][82] species would be the least impacted by climate change. Therefore, we studied the ecological niches of these two species in the multidimensional space of climatic factors, determined the main factors limiting these two species of ground beetles, and evaluated the changes in their distribution in the context of global warming.

Materials and Methods
This study is based on both materials collected by the authors and records retrieved from the Global Biodiversity Information Facility (GBIF). Field collections conducted over 1995-2019 in various areas of Ukraine allowed us to assemble 75 points of Z. tenebrioides occurrence and 105 points of H. rufipes occurrence. Studies were carried out in 2008-2018 in the territory of the Czech Republic, covering 22 biotopes, including 14 in the steppe zone of the Republic and 8 in the pre-mountain forest-steppe zone. A total of around 7 thousand specimens was collected using pit-fall traps, soil sampling and manual collection. In each Diversity 2021, 13, 559 3 of 15 biotope, 20 traps were set, and the material was collected once a decade. The dominating species in agrocoenoces was Z. tenebrioides, some individuals of which were also found in the pre-mountain area. Harpalus rufipes was recorded ubiquitously, which were most dominant in the steppe zone of the Republic, particularly in flood-plain forests and forest belts during hot periods.
Therefore, we used 692 points of occurrence of Z. tenebrioides and 1672 points of the occurrence of H. rufipes. For a more even distribution of the points of occurrence, the initial sampling was confirmed through filtration in RStudio software (R Core Team, Boston, MA, USA,), after which there were 435 points remaining for Z. tenebrioides ( Figure 1) and 653 remaining for H. rufipes ( Figure 2). Duplicate records were removed manually, and the thinning was performed, to filter points, using, the R script "thin points 4-1-18.R" package spThin (Functions for Spatial Thinning of Species Occurrence Records for Use in Ecological Models) (https://github.com/mlammens/spThin/issues, accessed on 23 March 2021). thousand specimens was collected using pit-fall traps, soil sampling and manual collection. In each biotope, 20 traps were set, and the material was collected once a decade. The dominating species in agrocoenoces was Z. tenebrioides, some individuals of which were also found in the pre-mountain area. Harpalus rufipes was recorded ubiquitously, which were most dominant in the steppe zone of the Republic, particularly in flood-plain forests and forest belts during hot periods. An additional 617 points of Z. tenebrioides and 1,567 points of H. rufipes were obtained from GBIF (doi:10.15468/dl.djj7bb; doi:10.15468/dl.8eegbx).
Therefore, we used 692 points of occurrence of Z. tenebrioides and 1,672 points of the occurrence of H. rufipes. For a more even distribution of the points of occurrence, the initial sampling was confirmed through filtration in RStudio software (R Core Team, Boston, MA, USA,), after which there were 435 points remaining for Z. tenebrioides ( Figure 1) and 653 remaining for H. rufipes ( Figure 2). Duplicate records were removed manually, and the thinning was performed, to filter points, using, the R script "thin points 4-1-18.R" package spThin (Functions for Spatial Thinning of Species Occurrence Records for Use in Ecological Models) (https://github.com/mlammens/spThin/issues, accessed on 23 March 2021).
.   thousand specimens was collected using pit-fall traps, soil sampling and manual collection. In each biotope, 20 traps were set, and the material was collected once a decade. The dominating species in agrocoenoces was Z. tenebrioides, some individuals of which were also found in the pre-mountain area. Harpalus rufipes was recorded ubiquitously, which were most dominant in the steppe zone of the Republic, particularly in flood-plain forests and forest belts during hot periods. An additional 617 points of Z. tenebrioides and 1,567 points of H. rufipes were obtained from GBIF (doi:10.15468/dl.djj7bb; doi:10.15468/dl.8eegbx).
Therefore, we used 692 points of occurrence of Z. tenebrioides and 1,672 points of the occurrence of H. rufipes. For a more even distribution of the points of occurrence, the initial sampling was confirmed through filtration in RStudio software (R Core Team, Boston, MA, USA,), after which there were 435 points remaining for Z. tenebrioides ( Figure 1) and 653 remaining for H. rufipes ( Figure 2). Duplicate records were removed manually, and the thinning was performed, to filter points, using, the R script "thin points 4-1-18.R" package spThin (Functions for Spatial Thinning of Species Occurrence Records for Use in Ecological Models) (https://github.com/mlammens/spThin/issues, accessed on 23 March 2021).
.   The potential range was predicted using Maxent 3.4.4 software (http://biodiversity informatics.amnh.org/open_source/maxent, accessed on 5 March 2021). The capabilities of this rapidly developing software are constantly improving, and resultantly, the quality of the models it generates. Maxent is now recognised as one of the best algorithms for modeling ranges of species. As mentioned above, we used the RStudio packages and SDMtoolbox applications to improve the model with regard to occurrence data, as well as climatic characteristics. According to a number of authors, a high correlation between variables makes it harder to estimate the contribution of each variable to the developed model of species distribution, and in some cases, it is impossible to do so. Therefore, we used "Raster" RStudio package, which detects and deletes highly correlated variables. In the results correlation file (csv), a matrix was obtained, an analysis of which allowed us to identify all variables above and below 0.75. Those variables scoring higher than 0.75 were removed from the analysis [83][84][85][86][87][88].
To avoid a shift of the predicted range towards territories that were studied to a higher degree (those with higher number of occurrence points), we used the ENMEval package in the R. A bias file, which was developed and used in Maxent modeling.
ENM-eval generates the csv file "enmeval-results" through which the "best" model can be selected, based on one or more evaluation criteria. We used the model settings that resulted in the lowest AICc. A simulation was performed according to the settings that corresponded to the minimum AICc value and to ∆AIC = 0.
To generate the model, we used predictor functions, linear features and quadratic features. To obtain independent datasets, we used spatial cross-validation. We set 10 replicates. To predict the probability occurrence and to model the spatial distributions of Z. tenebrioides and H. rufipes, we used a complimentary log-log (cloglog) function. We set a maximum of 500 iterations for all the cases. When developing the predicted ranges, we turned off the "extrapolation" feature because we did not know how the species would react to new conditions.
To verify the developed models, we used AUC under ROC-curve. The accuracy parameters of the model were evaluated using random test samples. To evaluate the obtained model, we set the random test percentage at 25, i.e., the program randomly selected 25% of the entire range of points so as to test the developed model for predicting performance. Therefore, 75% of the data for the model was used as a training sample and 25% as a test sample. To evaluate the predicting performance of the model, we used a threshold of prediction binarization whereby the conditions were considered favourable for species if the prediction had been higher than a certain threshold value, and unsuitable if lower. A threshold of 10% was used for the modelling in order to exclude individuals living in untypical habitats from the analysis. This means that 90% of the points of occurrence included in the analysis are in the "potential range", while the remaining 10% that are outside of this area are considered as living in atypical conditions and were therefore not taken into account during the identification of their climatic niche.
The AUC measures the property of the model to distinguish between places where the species is present and places where it is absent, alternating between 0 and 1. In our model, the average AUC value for Z. tenebrioides is 0.933 ± 0.011, i.e., the reliability of the generated model is quite high (Figure 1).
For all cases of modelling, we used the climatic model CCSM 4 (Community Climate System Model). We took into account two scenarios of climate change: the RCP2.6 (increase of the temperature on the planet is considered to be 0.9 • C on average by 2100) and the RCP8.5 (growth by 4.1 • C). The work on the layers and the estimation of the areas was performed using QGIS 3.18.0 software (Quantum GIS, www.qgis.org, accessed on 5 March 2021).

Results
According to our analysis, the most significant factors for Z. tenebrioides are the bio 15-the variation coefficient of seasonal pattern of precipitations, bio 17-the precipitations of the driest quarter, bio 18-the precipitations of warmest quarter, bio 19-the precipitations of the coldest year quarter, bio 2-the mean range of daily temperatures, bio 7-the annual temperature range, bio 8-the mean temperature of the most humid quarter, bio 9-the mean temperature of driest quarter.
The most significant factors determining the distribution Z. tenebrioides are, the precipitations of the driest quarter, the variation coefficient of the seasonal regime of precipitations, annual temperature range and the precipitations of the coldest quarter (Table 1). The jackknife method is based on the comparison of models that were developed for each of the factors and models that were developed without the factors. If two factors are related, the exclusion of one of them would reduce the number of errors in the model.
The temperature and proportion of moisture characterize the climate of a particular territory, which is important for population preferences. Moisture is usually the primary factor influencing the distribution of the species near the southern borders of the range, while the northern borders are, to a higher degree, determined by temperature criteria.
The results we obtained suggest that the Z. tenebrioides species prefers humid conditions in a dry period and is negatively affected by large amounts of precipitations in winter. The increase in annual temperature fluctuations also reduces the distribution range of Z. tenebrioides (Figure 3).
According to the developed model, the most favourable territories for Z. tenebrioides, with the maximum favourability value (0.80-0.98), are in the east coast of Scotland and southeast England, specifically in arable lands with wheat and barley. The conditions in France, Belgium, the Netherlands and Denmark are highly favourable for Z. tenebrioides. This is related to the fields of grain crops that exist in these areas, among other factors. Southeast and northeast Germany is characterised as highly favourable for Z. tenebrioides and has fields with grain crops. Southeast and Eastern Europe include both more favourable and less favourable territories, and in most cases the favourability of a territory correspond to the presence or absence of arable fields. In Ukraine, the most suitable zone is the steppe, and as we progress northwards, (the forest-steppe zone and the Polissia) the favourability coefficient decreases. Favourable territories are found in Krasnodar Krai and Stavropol Krai (Figure 4). and has fields with grain crops. Southeast and Eastern Europe include both more favourable and less favourable territories, and in most cases the favourability of a territory correspond to the presence or absence of arable fields. In Ukraine, the most suitable zone is the steppe, and as we progress northwards, (the forest-steppe zone and the Polissia) the favourability coefficient decreases. Favourable territories are found in Krasnodar Krai and Stavropol Krai (Figure 4). The models of predicted bioclimatic ranges for 2050 and 2070 in the view of two climatic scenarios are presented in Figure 5. According to the estimated areas, both of the scenarios predict a decrease in the bioclimatic range of Z. tenebrioides by 2070 (Table 2).   The models of predicted bioclimatic ranges for 2050 and 2070 in the view of two climatic scenarios are presented in Figure 5. According to the estimated areas, both of the scenarios predict a decrease in the bioclimatic range of Z. tenebrioides by 2070 (Table 2).  The model, generated according to the RCP 2.6 scenario, indicates that the bioclimatic range of Zabrus tenebrioides will increase by 334,792 km2 by 2050, and reduce by 1,195,022 km2 by 2070. According to RCP 8.5 scenario, the range will decrease by 1,020,602 km2 by 2050, and slightly increase by the year 2070 compared with 2050. At the same time, the area with the most favourable conditions will decrease.
The modelling also generated maps of potential and predicted ranges of H. rufipes. It is especially interesting to analyse the results that indicate the contributions of the climatic parameters to the modelling of their potential range. This data allowed us to infer the significant ecological factors that limit the distribution of the species. The most significant factors influencing the distribution of this species are bio 1-the mean annual temperature, bio 2-the mean daily range of temperatures, bio 3-the isothermality, bio 6-the lowest temperature of the coldest month, bio 10-the mean temperature of the warmest quarter (Table  3).  The model, generated according to the RCP 2.6 scenario, indicates that the bioclimatic range of Zabrus tenebrioides will increase by 334,792 km 2 by 2050, and reduce by 1,195,022 km 2 by 2070. According to RCP 8.5 scenario, the range will decrease by 1,020,602 km 2 by 2050, and slightly increase by the year 2070 compared with 2050. At the same time, the area with the most favourable conditions will decrease.
The modelling also generated maps of potential and predicted ranges of H. rufipes. It is especially interesting to analyse the results that indicate the contributions of the climatic parameters to the modelling of their potential range. This data allowed us to infer the significant ecological factors that limit the distribution of the species. The most significant factors influencing the distribution of this species are bio 1-the mean annual temperature, bio 2-the mean daily range of temperatures, bio 3-the isothermality, bio 6-the lowest temperature of the coldest month, bio 10-the mean temperature of the warmest quarter (Table 3). Note: Bio 1-the mean annual temperature, bio 2-the mean daily range of temperatures, bio 3-the isothermality, bio 6-the lowest temperature of the coldest month, bio 10-the mean temperature of the warmest quarter. AUC-the parameter of the factor's significance.
The optimum mean annual temperature for the distribution of H. rufipes is within 5-20 • C the mean daily range of the temperatures equals 6 • C the isothermality is within 30-35 • C the lowest temperature of the coldest month is −5-0 • C the mean temperature of the warmest quarter is 15-25 • C ( Figure 6). Note: Bio 1-the mean annual temperature, bio 2-the mean daily range of temperatures, bio 3-the isothermality, bio 6-the lowest temperature of the coldest month, bio 10-the mean temperature of the warmest quarter. AUC-the parameter of the factor's significance.
The optimum mean annual temperature for the distribution of H. rufipes is within 5-20 °C; the mean daily range of the temperatures equals 6 °C; the isothermality is within 30-35 °C; the lowest temperature of the coldest month is -5-0 °C; the mean temperature of the warmest quarter is 15-25 °C (Figure 6).
The model demonstrates that, in the present-day climatic conditions, based on the existing set of the data, H. rufipes may occur throughout Europe, covering northern territories of Africa and spreading to China. The range of H. rufipes, with some exceptions, coincides to a high degree with the distribution of Z. tenebrioides. The most favourable territories for this species are in Eastern Scotland and Southeast England, Southern Sweden, France, Germany, the Baltic countries, Southeast and Eastern Europe. Suitable climatic conditions for the species also exist in countries within Central Asia (Figure 7).  The model demonstrates that, in the present-day climatic conditions, based on the existing set of the data, H. rufipes may occur throughout Europe, covering northern territories of Africa and spreading to China. The range of H. rufipes, with some exceptions, coincides to a high degree with the distribution of Z. tenebrioides. The most favourable territories for this species are in Eastern Scotland and Southeast England, Southern Sweden, France, Germany, the Baltic countries, Southeast and Eastern Europe. Suitable climatic conditions for the species also exist in countries within Central Asia (Figure 7). The predicted models indicate a decrease in the overall bioclimatic range and in the range of territories that are the most favourable for this species (Figure 8). This is confirmed by the estimation of the areas of present-day and predicted bioclimatic ranges ( Table 4). a b c d  The predicted models indicate a decrease in the overall bioclimatic range and in the range of territories that are the most favourable for this species (Figure 8). This is confirmed by the estimation of the areas of present-day and predicted bioclimatic ranges (Table 4).

Discussion
Our data, which suggests that winter precipitations, precipitations during the driest quarter of the year and the seasonality of precipitations have a great effect on Z. tenebrioides correlates with the data collected by other authors. Most researchers confirm the ideawhich was developed a long time ago-that the factors that limit the population level of Z. tenebrioides are precipitations, soil moisture and low temperatures in which the larvae overwinter [89,90].
The predictive maps, based on determined ecologic-climatic niches, which we obtained using MaxEnt algorithm, generally coincide with the distribution of the studied species across the continent. Almost the entire territory of Western and Central Europe has favourable conditions for the dispersal of Z. tenebrioides, which is coherent with the literature's data [89,91]. Some territories, Belarus for example, are less favourable, but nonetheless inhabitable territories. According to Trepashko and Bojko [92], Z. tenebrioides had not been discovered within Belarus Republic territory until 2015. In the spring of 2016, in the fields of winter triticale in Southern Belarus, large growth centers of Z. tenebrioides larvae were discovered for the first time [93]. For Z. tenebrioides, precipitations are crucial, as confirmed by the studies of Beleckij and his students [89,91,94,95]. They note that one of the limiting factors is drought. A multi-year reproduction analysis of the species in Ukraine suggests that the outbreaks occurred in the years of drought in 1863-1979, droughts in 1991-1993, and severe droughts in 2003-2007 (especially in 2003 and 2007) [90,93,94]. The temperature optimal for Z. tenebrioides ranges between +20 and +26 • C; whereas an increase in temperature of up to +30 • C decreased the activity of the beetles, and temperatures of +36 • C caused their death. During hot periods, the beetles hide in forests, fractures in soil and other shelters. The results of the study revealed that the population of Z. tenebrioides in agrocoenoses decreased during hot periods and increases three-fold in neighbouring forests.
Depending on the weather, a diapause of Z. tenebrioides lasts for an average of 20-30 days. After precipitations and temperature drops, the beetles emerge to the soil surface (second half of August-early September). Females lay eggs in August-September. A dissection of females revealed 6-14 eggs. Starting from mid-June, the number of generative females in the pit fall traps has increased, peaking in August-September. The larvae overwinter in fields of winter crops and have low cold resistibility, which limits their range. In spring, larvae feed on winter wheat, and then descend into soil. Pupation takes place in late April; the pupa phase lasts for 15-30 days.
Hilevskij [95], who studied Z. tenebrioides, notes that "in order to develop, the first instar larvae need the overall effective temperature to equal 34.5 • C, second instar larvae −37.5 • C, third instar −47.5 • C, pupa −18.0 • C. In September-October, when the mean daily temperature is 6-10 • C, the development of the first instar larvae completes in 30-45 days, and the feeding period lasts for 15-20 days. Then the larvae descend into the soil to a depth of 20-30 cm to mold which lasts 5-7 days. Molded larvae are low-active for another 5-8 days, and then begin to actively feed. The development period of the second and third instar larvae (the most harmful)-depending on the weather conditions, is 30-50 days, and the period of active feeding takes 15-20 days" [95].
Heavy precipitations following long droughts affect the early development of oocytes in the ovaries of Z. tenebrioides, leading to their resorption. Such weather conditions cause the unproductive expenditure of the reserve substances and a decrease in fertility of female beetles. Reductions in the soil moisture levels impede the exit of beetles to the surface and decreases their fertility and their periods of oviposition. The activity of beetles resumes following a decrease of precipitations and a temperature reduction, usually in the second half of August-early September [96]. Factors that limit their reproduction include a critical decrease in soil temperature in late spring and winter and drought during the emergence of larvae from the egg chamber, which is coherent with the data we obtained. Eggs do not develop in low humidity. The distribution of the ground beetles in areas of cultivated winter wheat is limited by the following climatic factors: a decrease in the mean soil temperature at a 20 cm depth during the coldest month to −3 • C and a decrease in the annual level of precipitations amounting to less than 400 mm.
While the significant factors for Z. tenebrioides are those related to precipitations (for the beetle inhabits open spaces), the significant factors influencing H. rufipes are temperature related. Moisture deficiency in the southern part of the range is compensated for by movement to forests or through their diet. Harpalus rufipes is a polytopic mesophile, which inhabits forest biotopes and anthropogenic landscapes. It is reported as an effective entomophage that kills the Colorado beetle, Sitona and owlet moths. According to our data, H. rufipes develops in one year in the southern part of the range, and in two years in the northern part. Imagoes occur from March-April to late October. Oviposition takes place from early May to July. Given its a two-year development cycle, the overwintering phases occur for both larva and imago.
Climate change causes mean annual temperatures to rise, and therefore leads to a change in other bioclimatic parameters, which in turn would lead to a change in the species range in southern inner continental territories, causing an insignificant northward shift.

Conclusions
In conclusion, global climate changes are becoming more notable every year, and may significantly shift arable crops lands northwards, and, naturally, decrease the yield in the areas with less fertile soils. The two most abundant species of ground beetlespests Z. tenebrioides and H. rufipes-will evidently shift their ranges northward. According to the climatic models we developed, the ranges of both species would reduce in area in their eastern, most continental, parts. Controlling the populations of Z. tenebrioides and H. rufipes and creating a detailed study of changes in their interactions with other elements of biocoenoses (plant species the insects consume, predators, parasites, disease pathogens) would help us to better understand the future changes in the system of "grain crops-phytophage ground beetles".

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.