GIS-based modelling reveals the fate of antlion habitats in the Deliblato Sands

The Deliblato Sands Special Nature Reserve (DSSNR; Vojvodina, Serbia) is facing a fast successional process. Open sand steppe habitats, considered as regional biodiversity hotspots, have drastically decreased over the last 25 years. This study combines multi-temporal and –spectral remotely sensed data, in-situ sampling techniques and geospatial modelling procedures to estimate and predict the potential development of open habitats and their biota from the perspective of antlions (Neuroptera, Myrmeleontidae). It was confirmed that vegetation density increased in all parts of the study area between 1992 and 2017. Climate change, manifested in the mean annual precipitation amount, significantly contributes to the speed of succession that could be completed within a 50-year period. Open grassland habitats could reach an alarming fragmentation rate by 2075 (covering 50 times less area than today), according to selected global climate models and emission scenarios (RCP4.5 and RCP8.5). However, M. trigrammus could probably survive in the DSSNR until the first half of the century, but its subsequent fate is very uncertain. The information provided in this study can serve for effective management of sand steppes, and antlions should be considered important indicators for conservation monitoring and planning.

and bushes in Middle Europe 25 and South Africa 26,27 , and some antlions are documented as indicators of highly endangered habitats in South Africa 28,29 . In South Africa, a research programme has been undertaken with the goal of monitoring the status of antlions and related lacewings to ensure the conservation and survival of these particularly rich and unique south-African habitats 29 . The results are being used to develop protective measures, such as nature reserve areas with protected species, definition of areas that are inhabited by endangered species, or location and identification of areas that need protection because they are particularly rich in antlions 28,29 . Based on these facts and examples of good practice, we aimed to develop a predictive geospatial model to evaluate the prospects for the DSSNR and its antlions (Myrmecaelurus trigrammus (Pallas, 1771) in particular) if conservation actions are not implemented. Geospatial data (environmental and socioeconomic) and complex predictive models designed in the geographic information system environment are gaining acceptance among ecological studies 30 because they provide applicable results with high value for environmental planners and decision makers in the field of conservation biology 31,32 . It is thus crucial to bridge the gap between geospatial science and ecology to apply innovative spatial approaches to support for effective conservation planning strategies in all vulnerable ecosystems. In that regard, the following research questions were addressed: (1) How rapid is the succession process in the DSSNR? (2) How many open habitats have been lost over the last 25 years, and how many will survive by the end of the century? (3) How will the succession process, triggered by land abandonment and climate change, affect the existing antlion population in future years? And (4), can antlions be considered as an indicator species in the sand steppes for conservation monitoring and planning?
Material and Methods insect sampling. This study used antlion species, characteristic of sand steppe and grassland habitats.
Although most adult antlions are night active (e.g. Myrmeleon formicarius (Linnaeus, 1767); Euroleon nostras (Geoffroy in Fourcroy, 1785); Distoleon tetragrammicus (Fabricius, 1798)), they can also be easily collected during the day with an insect net 33 . Collecting larval antlions in the DSSNR study area would not be successful, since only the larvae of Myrmeleon inconspicuous (Rambur, 1842), M. formicarius and E. nostras construct pits, while the rest of the species listed in this study are non-pit-builders. We sought to use a robust sampling technique that could be carried out by employees of similar protected areas around the World in order to ensure the replicability, transferability, and above all, applicability of our results. Thus, random sweep-netting turned out to be the most convenient sampling method. Moreover, the focus of this research is not oriented towards species composition change on selected sites but to the currently most abundant diurnal antlion species Myrmecaelurus trigrammus (Pallas, 1771) (hereafter M. trigrammus), a potential indicator for regionally rare open sandy grassland habitats. Such an approach minimizes the sampling error and thus excuses the robust methodological procedure employed.
Adult antlions were collected in the study area within the DSSNR in July 2016. Open grassland habitats (N = 7) of varied size were selected and GPS recorded with a GARMIN GPS device (Fig. 1b). The abundance of antlions was determined by counting captured individuals on these sites, collected by the sweep-netting method for one hour by three persons in a plot of 50 m × 100 m. All open grassland habitats were visited at the same hour (12:00) of the day. After identification to a species level, specimens were returned to their natural habitat. For antlion identification, we used the keys in Badano and Pantaleoni 34 and Aspöck et al. 35 and some collected larvae excavated by spoon. However, the most abundant antlion species (occurring in all sampled open grassland habitats) -M. trigrammus -was then chosen for assessment of its current and potential future abundance in the whole DSSNR study area.
Time series file development. To detect the successional dynamics across the whole study area of the DSSNR, freely available cloud-free LANDSAT multispectral imagery was obtained from the EarthExplorer web . Source: Own data and ArcGIS Ocean and Imagery (WGS84) base maps URLs: https://www. arcgis.com/home/item.html?id=52bdc7ab7fb044d98add148764eaa30a, https://www.arcgis.com/home/item. html?id=5ae9e138a17842688b0b79283a4353f6. platform (https://earthexplorer.usgs.gov/). Data from 1992 to 2017 were downloaded and converted to TerrSet's 36 time series (TSF) format. However, for some time windows (1993, 1995-1999, 2002, 2004, 2008 and 2012), satellite images were not available, owing to cloud cover or sensor error in the target months (July and August). Thus, LANDSAT 4, 5, 7 and 8 products were applied to develop the first time series database (NDVI) containing 17 time frames.
Environmental and/or anthropogenic factors can accelerate or decelerate a succession process. Modern climate change, triggered by human activity, can be considered as a key determinant from both perspectives. We tested whether climate variables contribute significantly to the succession trend in the DSSNR study area. The CHELSA geospatial climate database (http://chelsa-climate.org/) was used to produce annual mean air temperature (T) and precipitation variables (RR). These were gathered for almost the same period as the satellite data (1993-2016, [2017 not yet available on the CHELSA database]) in order to compare climate variability and biomass change dynamics in the DSSNR. For that purpose, a second time series database (current climate) was developed.
A pixel level regression. Initially, the visible red and near-infrared bands (30 m horizontal resolution) were used to calculate the Normalized Difference Vegetation Index (NDVI) for each time window (17) by applying the VEGINDEX tool in Terrset software 36 . NDVI is a vegetation index commonly used in ecological studies because it carries information about vegetation density or biomass and minimizes the topographic effect, while producing a linear measurement scale 38,39 . The NDVI value for a given pixel always returns a number ranging from minus one (−1) to plus one (+1); however, the absence of green leaves gives a value close to zero, whereas zero means no vegetation, and values close to +1 (0.8-0.9) indicate the highest possible density of green leaves 38 . Prior to the pixel level regression analysis, mean NDVI values, mean T and mean RR values for 16 time windows (from 1993 to 2016) were calculated for the whole DSSNR study area by applying the zonal statistics tool in ArcGIS 40 . Thus, the NDVI-time, NDVI-RR and NDVI-T relationships were obtained. Since time and RR had a significant linear impact on the NDVI dynamics in the DSSNR study area (Fig. 2a,b), both were used as NDVI predictors on a pixel level (30 m 2 ) basis. Before that, all the missing geospatial data were statistically interpolated by using the linear temporal interpolation algorithm within the Earth Trends Modeler (ETM) module in TerrSet 36 . In order to identify the speed of the overgrowing process, the linear NDVI trends with corresponding slope (β coefficient), R 2 and p-values were calculated with the predictors time ( Fig. 2c) and RR (Fig. 2d) for each pixel (30 m 2 ) in the study area by applying the Series Trend Analysis tool within ETM 36 . Both NDVI slope variables (time-and RR-derived) were then combined with the weighted sum algorithm by considering corresponding weights according to the residual standard error values resulting from the current NDVI-time and NDVI-RR relationship ( Fig. 2a,b). Thus, the final combined NDVI slope coefficient enabled not only the objective estimation of the possible disappearance of open habitats in the study area under several climate change scenarios, but also the estimation of potential future abundance of the target species M. trigrammus, which is linked to those open areas within the DSSNR. Finally, future time-derived and RR-derived NDVI variables were calculated and then summed with corresponding weights for scenarios RCP4. 5  Both EG variables were tested for possible significant correlation across the study area with the Band Collection statistics tool in ArcGIS to avoid possible variable redundancy 41 . Finally, a Poisson regression model with those two EG predictors was developed to estimate the abundance of M. trigrammus in the DSSNR study area and to predict its future abundance by considering the speed of the succession process triggered by land abandonment and climate change, manifested in the mean annual precipitation variability. The Poisson model (2020) 10:5299 | https://doi.org/10.1038/s41598-020-62305-3 www.nature.com/scientificreports www.nature.com/scientificreports/ was fitted in R statistical software 42 by applying the glm function. In addition, the residual deviance parameter was used to perform a goodness of fit test for the overall model. The residual deviance is the difference between the deviance of the current model and the maximum deviance of the ideal model, where the predicted values are identical to those observed. Therefore, if the residual difference is small enough, the goodness of fit test will not be significant, indicating that the model fits the data. After assuring a properly specified regression model, future EG predictors (2050 and 2075) were calculated by implementing the slope coefficients of the combined time-and RR-derived NDVI trends for each RCP scenario. Consequently, the geographic position and total area of open habitat across the study area in selected future time windows and climate change scenarios were estimated (Fig. 3)

Results
Antlion assemblages in the Deliblato sands. Ten (10) antlion species have been recorded in the DSSNR; among them, three species (Palpares libelluloides (Linnaeus, 1764), Acanthaclisis occitanica (Villers, 1789) and Nohoveus punctulatus (Steven in Fischer v. Waldheim, 1822)) are suspected of being locally extinct in Serbia (Table 1). Among recently detected antlions, three species frequently occurred during the survey in July 2016: Myrmecaelurus trigrammus, Creoleon plumbeus (Olivier, 1811) and Myrmeleon inconspicuus. One of these, M. trigrammus, was considered as the most appropriate for geospatial modelling, given that it had the highest abundance in each sampled open grassland habitat. However, all three frequently occurring antlions are indicators of sand steppes in the Pannonian Basin. Based on these results and existing literature, we estimated the conservation status of antlion species in Serbia ( Table 1).

pace of succession. Historical maps (The First and the Second military survey of the Habsburg Empire)
revealed that our study area was practically a sandy desert-like environment by the end of the 18 th century and was sparsely covered with wooded plants by the second half of the 19 th century.
However, over the last 25 years, the forest has constantly progressed and is slowly fragmenting the existing open grassland habitats in the DSSNR (Fig. 1b). The pixel level regression analysis revealed the succession dynamics in detail. Time and mean annual precipitation (RR) are key natural determinants controlling the speed of the succession process in the study area (Fig. 2). Thus, the NDVI-time and NDVI-RR relationship, and the developed geospatial time series datasets (NDVI, current climate and future climate), enabled the objective estimation of future open grassland habitat distribution in the DSSNR.
In 2017, open grassland habitats covered 13815 hectares (23.5%). By considering the RCP4.5 climate scenario, open grassland habitats in the DSSNR study area could occupy nine times less area (1580 Ha or 2.69%) by 2050 or only 274 hectares (0.47%) by 2075 (Fig. 3a). The RCP8.5 scenario, the most pessimistic from the greenhouse gas emission perspective, offered a slightly better outcome for open areas in the DSSNR. Nevertheless, by 2050 these could shrink by 11.7% compared to the situation in 2017, covering 1614 hectares by then. The fragmentation process of open grassland habitats in the study area could potentially reach an alarming state by 2075, if RCP8.5 or a similar emission scenario is realized, covering an area of only 284 hectares (0.48%) (Fig. 3b). Since  www.nature.com/scientificreports www.nature.com/scientificreports/ slowly overgrown with time but are still within the open habitat threshold value. In contrast, the change in the mD2F EG variable is much more intensive. Because open habitats are losing area and are increasingly invaded by forest, the distance variable is taking over lower values in each successive future time window.

Ecogeographical (EG) variables and the Poisson regression model.
However, the weak correlation (r < 0.3) between the two predictor variables allowed the proper fitting of a Poisson regression model (Table 2), which was later applied in the geographical space to predict the potential future abundance of M. trigrammus across the study area. The insignificant goodness-of-fit chi-squared test (p > α, α = 0.05), based on residual deviance, additionally proved that the model was properly specified. Both predictor variables had a significant influence on the dependent variable (M. trigrammus abundance). As expected, estimates (log value) were negative, since higher abundance was clearly linked to larger open areas that were less densely vegetated.

Discussion
Massive afforestation of the DSSNR during the last two centuries had strong deleterious effects on the ecosystem composition, where more than three-quarters of the former natural habitat has been lost. While significant erosion control is a feasible positive effect, habitat loss is likely to be a negative consequence of afforestation in grasslands and sand dune ecosystems 44 . Habitat loss, as a consequence of afforestation and natural succession, is a crucial threat to the DSSNR. Maintenance of a steady state within an ecosystem requires continuous monitoring of biodiversity and ecosystem function 45 . In steppe grasslands, this can be realized by monitoring indicator species of antlions, whose population levels are easily measured. Antlions are sand-and grassland dwellers with behavioural, morphological and physiological adaptations for their environment 28 . They could have a substantial value for conservation needs 25,28,29 . It is widely accepted that a nationally rare or threatened species present in a surveyed area has a higher conservation value than a ubiquitous species because it contributes more to regional biodiversity than a common species 46,47 . Dry grasslands and sand steppes harbour a high proportion of these insects, whose larvae require a sandy microhabitat for completion of their larval development.
Based on our results, M. trigrammus should be regarded as an indicator species, characteristic of the sand steppe in the southern rim of the Pannonian Plain. Outside the Mediterranean, another two antlion species (C. plumbeus and M. inconspicuus) commonly occur in the sand steppe of the Pannonian Plain, in sand dunes or the remnants of formerly widespread sand dunes. Both have considerable potential to be considered as other ecological indicator species in these environments. However, they are rare in adjacent countries 25,[48][49][50][51][52][53][54][55][56] . Two antlion species, A. occitanica and N. punctulatus, which have not been found in the DSSNR since the beginning of the 20 th century, are declining throughout their range and are near extinction in some European countries 54,57,58 . One rare species, P. libelluloides, was found in the DSSNR only once; it is an insufficiently known antlion, probably extinct in the reserve 59 .
Our model predicts that natural succession in the form of afforestation could be completed within a 50-year period. This prediction leads to the conclusion that all grassland-linked taxa, including the antlions typical of the sand steppe, could become locally extinct. However, this process could be decelerated by permanent temperature increase and precipitation reduction, trends detected in the last half of the century 60 . We tested this assumption by integrating past, recent and future geospatial climatic data with the remotely sensed biomass change dynamics in the DSSNR study area. Some direct climate change impact on the target species M. trigrammus and others, can be expected in the coming years, but these will be negligible compared to those climate drivers influencing the www.nature.com/scientificreports www.nature.com/scientificreports/ succession process. The key factor in this case is mean annual precipitation, which is positively correlated with the biomass gain rate in the study area. Although the selected global climate models and emission scenarios predict less precipitation for the DSSNR in the second half of the 21 st century, open grassland habitats could reach an alarming fragmentation rate by 2075. Nature conservation actions are urgently needed if we want to preserve the DSSNR's current biodiversity.
From this perspective, a partial restoration of the DSSNR would be the only way to prevent the disappearance of natural habitats along with their fauna and flora. These restored open habitats must then, of course, be properly maintained. In this regard, fires have a positive effect in maintaining grassland habitats. Natural fire is considered as a natural disturbance that has many benefits for grassland ecosystems 61 . It can be a significant evolutionary force, which has driven the evolution of species 62 . Many plant species require natural fire to germinate and reproduce. Thus, a carefully chosen and controlled fire regime has a role in maintaining biodiversity in grassland or savanna ecosystems 63 . Moreover, controlled or prescribed burning is a tool that is receiving additional attention in ecosystem studies and protected area management. Dormant-season prescribed burning was applied to grasslands in East-Hungary in order to test its effects from the nature conservation perspective 64 . This experiment resulted in changed soil characteristics along with changes in plant biomass and composition. On the burnt sites, soil organic matter, pH, potassium and phosphorous did not change significantly, but soluble salt content increased drastically 64 . Other positive effects were also detected: an increased number of flowering shoots and plant diversity, along with stable arthropod abundance and diversity 64 .
In contrast to grazing, which is accompanied by nitrification of the soil, consequently influencing grassland species composition 65 , forest fires have no effect on that process. Indeed, forest fires and cattle grazing did occur in the past in the DSSNR 23,24 . Finally, we can conclude that controlled burning could be the most appropriate grassland management tool in the case of the DSSNR, because it supports plant diversity and does not threaten the invertebrate fauna; such conservation actions should be implemented there as soon as possible, in order to avert high local extinction risk and the worst-case scenario.

Data availability
All the data and results in the form of tables and maps are available on request.