The Spread of the Invasive Locust Digitate Leafminer Parectopa robiniella Clemens, 1863 (Lepidoptera: Gracillariidae): The Ukrainian Context

: The spread of phytophagous pests are often associated with global warming. These species may be of interest in terms of biological indications of climate change. We considered the locust digitate leafminer P. robiniella . In Ukraine it was ﬁrst recorded in 2003. In 2020–2021, we found areas of massive R. pseudoacacia leaf damage caused by the pest on Trukhaniv island (Kyiv) and some places in the region. Using 592 georeferenced records of P. robiniella across Europe and a Bayesian additive regression trees algorithm, we modeled the distribution of the moth. As predictors, a current climate (WorldClim v.2) and a black locust habitat suitability map were employed. Western and south-western regions of Ukraine and Transcarpathia are considered the most favorable for the pest. Amongst the factors determining its niche, summer moisture and warm conditions are the most important for facilitating the spread and naturalization of the moth. Under progressing climate change, the species is expected to move eastward.


Introduction
Earth's temperature has risen by 0.08 • C per decade since 1880, and the rate of warming over the past 40 years is more than twice that: 0.18 • C per decade since 1981 [1]. Numerous studies provide evidence for biological responses to recent climate change [2] and, in particular, focus on changes in the distribution and range shifts of species [3]. For pest insects, it is expected that climate warming will affect both the incidence and geographical extent and intensity of population outbreaks, with potentially severe economic and ecological consequences [4]. In this context, we considered the locust digitate leafminer Parectopa robiniella Clemens, 1863 (Lepidoptera: Gracillariidae) as a model species indicating climate change. It is native to North America, but was accidentally introduced to Italy, where it was first found in Milano in 1970 [5]. The moth has now been recorded from Italy, France, Germany, Slovenia, Croatia, Austria, Hungary, Slovakia, Romania, Lithuania, Ukraine and Russia [6][7][8][9][10][11][12][13][14]. On average, the spread of the pest occurs at a speed of about 100 km per year [15]. The moth is associated with the black locust (Robinia pseudoacacia L.), a tree introduced to Europe from North America at the beginning of the XVII century for decorative reasons, and to Ukraine it was brought at the end of the XVIII century [14]. In some countries, such as Romania and Hungary, it has become an essential landscape element [16]. In Ukraine too, particularly in the south, the tree has been widely cultivated for a variety of purposes. R. pseudoacacia is considered an excellent plant for growing in highly disturbed areas as an erosion control plant [17], but also not in the least because the black locust is a major honey-producing plant [18].
In Ukraine the species was recorded for the first time in Kyiv and Chernivtsi (a regional capital in the south-west of the country) in 2003 [19]. In 2020, we discovered an extensive outbreak of the pest on Trukhaniv Island, a floodplain island located within the city boundaries of Kyiv. In 2021 a second such foci was found in Obukhiv district, some 50 km south of the capital city. P. robiniella larvae develop in chambers ("mines") inside the leaves, causing "finger-like" excavations all around the margins of the central blotch of the mine [20].
Although P. robiniella is now widely prevalent in Europe, there are heterogeneities in the pest's distribution that are likely linked to environmental parameters. In this respect, species distribution models (SDMs) have proven to be useful tools for predicting pest distribution and elucidating the importance of a wide range of environmental covariates considered to affect pest species occurrence [21]. Using SDMs, our objective was to identify areas of habitat suitability of the moth, identify conditions that constrain the geographic distribution of P. robiniella in Europe, particularly in Ukraine, against the background of climate change.

Materials and Methods
As input, SDMs require georeferenced biodiversity observations. Localities for P. robiniella were gathered from Global biodiversity information facility (GBIF) [22] online resource, records found in the literature, and our personal field surveys. Because many uncertainties are associated with SDM projections, particularly when it comes to building a SDM for a species expanding its home range in a new area [23], we used for the analysis only records of European localities. In total there were 592 such records.
SDMs commonly utilize associations between environmental variables and known species occurrence records to identify environmental conditions within which populations can be maintained. SDMs extrapolate in situ habitats both in space and time to obtain spatially explicit and continuous surfaces indicating the probability of species occurrence [24]. SDMs are primarily climate-driven, meaning that the variables used to develop them typically portray climatic factors [21]. This makes sense because climate is a chief driver of environmental suitability [21,25]. Information on the bioclimatic parameters was collected as raster layers from the WorldClim website (http://www.worldclim.com/version2) (accessed on 21 January 2022). These 19 bioclimatic variables indicate a general trend of precipitation and temperature, extremity and seasonality of temperature [26]. Although environmental conditions may predict a species' potential geographic distribution broadly, climatic factors typically used to represent those conditions are seldom adequate surrogates for factors such as host availability [21]. Therefore, we included to the set of environmental predictors a raster reflecting the habitat suitability of the black locust in Europe. For this purpose, occurrences were extracted from the GBIF database [27], yielding 137,294 records. For technical reasons, this number of records was reduced to 2353 using the "point thinning" module in System for Automated Geoscientific Analyses geoGraphic Information System (SAGA GIS) software [28].
SDMs were generated by employing Bayesian additive regression trees (BART), a powerful machine learning approach. Running SDMs with BARTs has recently been substantially facilitated by the development of an R package, "embarcadero" [29], including an automated variable selection procedure highly effective at identifying informative subsets of predictors. Additionally, the package includes methods for generating and plotting partial dependence curves, illustrating the effect of selected variables on habitat suitability. The algorithm computes habitat suitability values ranging from 0, for a fully unsuitable habitat to 1, for a fully suitable habitat. Model performance was assessed using measures of accuracy. Firstly, the threshold-independent receiver operating characteristic (ROC) approach, where the calculated area under the ROC curve (AUC, Area Under the Curve) is considered as a measure of prediction success. The ROC curve is a graphical method that represents the relationship between the false-positive fraction (one minus the specificity) and the sensitivity for a range of thresholds. It has a range of 0-1, with a value greater than 0.5, indicating a better-than-random performance event. A rough classification guide is the traditional academic point system: poor (0.5-0.6), fair (0.6-0.7), good (0.7-0.8), very good (0.8-0.9) and excellent (0.9-1.0). Secondly, the true skills statistic (TSS). TSS values of −1 indicate predictive abilities if not better than a random model; zero indicates an indiscriminate model and +1 a perfect model.
Against the background of the appearance of P. robiniella in Kyiv, positive dynamics of the average annual temperature in the city were hypothesized to act as one of the possible factors for promoting the invasion of this species; data from the E-OBS (European Observations) website [30] were used for this purpose.
Maps of habitat suitability in the GeoTIFF format were processed and visualized in SAGA GIS and time series data were processed in R [31].

Results
Because ecological niche models are sensitive to sample bias and spatial autocorrelation, which would produce models of lower rather than higher quality, initial sets or records, both for the moth and black locust, were filtered out using the corresponding module in the "embarcadero" package. In the end, the total number of records was reduced to 204 and 2334 for the moth and tree species, respectively.
Despite the notion that climatic factors typically used to represent habitat conditions are rarely adequate surrogates for factors such as host availability [23], both measures of accuracy showed each of the two SDMs built for the moth (one using bioclimatic predictors alone, while the other in addition incorporated the habitat suitability raster for the black locust) performed very well: AUC = 0.940 and TSS = 0.740, and AUC = 0.940 and TSS = 0.729, respectively. In addition, both models are highly correlated (R 2 = 94.20%, p < 0.05).
The automated variable selection procedure identified an informative subset of seven predictors ranked in order of importance for P. robiniella: (1) Bio 18-Precipitation of Warmest Quarter, an index providing the total precipitation (in mm) during the warmest three months of the year (in short, these are summer months in the study area), (2) Bio 8-Mean Temperature of Wettest Quarter, a quarterly index that approximates mean temperatures ( • C) that prevail during the wettest season (once again, these are summer months). Third in the rank, and not surprisingly, is the distribution in the study area of bioclimatic habitat suitability for the black locust, although the model obtained for the tree species itself performed relatively poorly, it was nonetheless in an acceptable range: AUC = 0.740, TSS = 0.368, and showed a strong dependence on temperature variables (Bio 7-Annual Temperature Range, featuring the seasonal amplitude in ambient temperature, and Bio 6-Minimum Temperature of Coldest Month).
Plots of partial dependence curves of identified informative predictors manifesting the model for P. robiniella are presented in Figures 1 and 2. Both curves clearly emphasize the dependence of summer precipitation and temperature on the well-being of the moth, where habitat suitability sharply increases after precipitation reaches the mark of around 250 mm, and the mean ambient temperature exceeds 15 • C. Interestingly, winter conditions turned out to be of less importance: Bio 6-Minimum Temperature of Coldest Month and Bio 14-Precipitation of Driest Month occupy the 5th and 7th ranks, respectively. Most likely this is because the moths successfully pass the winter at the pupal stage under the protection of covers of leaf litter [14,19]. Regarding the host plant, the partial dependence curve ( Figure 1A) shows a rapid increase in suitability for the moth in areas where suitability for the tree exceeds a predicted 50% (in other words, "better for the host-better for the pest").
For linking the appearance of P. robiniella in Kyiv with climatic dynamics, the average annual temperature in the city for 1950-2017 was downloaded from the E-OBS database. Average annual temperature has a strong uphill (positive) linear relationship with precipi-tation in the warmest quarter (Bio 18) (Pearson's r = 0.706, p < 0.05) and therefore can be considered a proxy of the former.
According to E-OBS, the average annual temperature in Milan in the 1970s, where P. robiniella was first detected, was 12.99 • C, which assumes that the local conditions were favorable for its naturalization, which, in fact, they were. At the same time, the average annual temperature in Kyiv was 8.73 • C, making the difference with Milan at that time 4.26 • C. Since 1988, there has been a significant (p < 0.05, 1000 bootstrap) trend of increasing average annual temperature in Kyiv ( Figure 1B) and now the difference with Milan in the 1970s is one and a half less: precisely 2.73 • C. For linking the appearance of P. robiniella in Kyiv with climatic dynamics, the average annual temperature in the city for 1950-2017 was downloaded from the E-OBS database. Average annual temperature has a strong uphill (positive) linear relationship with precipitation in the warmest quarter (Bio 18) (Pearson's r = 0.706, p < 0.05) and therefore can be considered a proxy of the former.
According to E-OBS, the average annual temperature in Milan in the 1970s, where P. robiniella was first detected, was 12.99 °C, which assumes that the local conditions were favorable for its naturalization, which, in fact, they were. At the same time, the average annual temperature in Kyiv was 8.73 °C, making the difference with Milan at that time 4.26 °C. Since 1988, there has been a significant (p < 0.05, 1000 bootstrap) trend of increasing average annual temperature in Kyiv ( Figure 1B) and now the difference with Milan in the 1970s is one and a half less: precisely 2.73 °C.   For linking the appearance of P. robiniella in Kyiv with climatic dynamics, the average annual temperature in the city for 1950-2017 was downloaded from the E-OBS database. Average annual temperature has a strong uphill (positive) linear relationship with precipitation in the warmest quarter (Bio 18) (Pearson's r = 0.706, p < 0.05) and therefore can be considered a proxy of the former.
According to E-OBS, the average annual temperature in Milan in the 1970s, where P. robiniella was first detected, was 12.99 °C, which assumes that the local conditions were favorable for its naturalization, which, in fact, they were. At the same time, the average annual temperature in Kyiv was 8.73 °C, making the difference with Milan at that time 4.26 °C. Since 1988, there has been a significant (p < 0.05, 1000 bootstrap) trend of increasing average annual temperature in Kyiv ( Figure 1B) and now the difference with Milan in the 1970s is one and a half less: precisely 2.73 °C.   In the 2000s, when P. robiniella was already detected in Kyiv, the average annual temperature was increasingly favorable for the species: 9.70 • C, and in the past decade it has risen to 10.26 • C. In light of our modelling exercises, we consider that the rising average annual temperatures (or more precisely, rising average summer temperatures), by large, are involved in setting the scene for the spread and naturalization of the pest, although in a recent study [32], a negative correlation was found between temperature and spread.

Discussion
Using the SDM and the 10th percentile threshold, locations can be found in Europe that are assumed under current climatic conditions to support the spread and naturalization of the moth. In Figure 2, these areas are shaded in dark green and cover a vast range in Western and Central Europe, much of Italy and the Balkans, reaching as far as the eastern Black Sea region. In Ukraine, particularly favorable for the pest are areas in the west and south-west of the country, and Transcarpathia bordering the Pannonian Plain. Under progressing climate change, the moth can be expected to move further towards the east.

Conclusions
SDMs seeking to identify features that characterize a species' known distribution are likely to provide basic quantitative information about a species' apparent habitat preferences and potential distribution. In our study, the extraction of such information has been accomplished using the SDM response curves for a number of environmental covariates identified as an informative subset of predictors. An emphasis is made on summer moisture and warm conditions facilitating the spread and naturalization of the considered pest.