Analysis of geographic centrality and genetic diversity in the declining grasshopper species Bryodemella tuberculata (Orthoptera: Oedipodinae)

Human-induced ecological and climatic changes have led to the decline and even local extinction of many formerly widely distributed temperate and cold-adapted species. Determining the exact causes of this decline remains difficult. Bryodemella tuberculata was a widely distributed orthopteran species before the mid-19th century. Since then, many European populations have suffered drastic declines and are now considered extinct or critically endangered. We used ecological niche modelling based on a large dataset of extant and extinct occurrence data to investigate whether poor climatic suitability in the periphery of its global range was a possible cause of the local extinction of the European populations of B. tuberculata. We also used population genetics based on the COI marker to estimate and compare the genetic diversity of extant populations. We found that Europe still provides highly suitable habitats close to the climatic optimum, contradicting the assumption of climate change as major driver of this decline. Instead, changes in land-cover and other anthropogenic modifications of the habitats at the local scale seem to be the major reasons for local extinctions. Genetic analysis suggests Central Asia as center of diversity with a stable population size, whereas the effective sizes of the remaining European populations are decreasing. We found European genetic lineages nested within Central Asian lineages, suggesting a Central Asian source distribution area. Our results suggest that the declining European populations represent relics of a formerly wider distribution, which was fragmented by changes in land-use. These relics are now threatened by limited connectivity and small effective population sizes. Specific conservation actions, such as the restoration of former or potential new habitats, and translocation of individuals from extant populations to these restored sites may help slow, stall, or even revert the extinction process.


Introduction
Biodiversity decline, specifically insect decline, has been recognized as a major threat to global ecosystems (Dirzo et al. 2014;Conrad et al. 2006;Hallmann et al. 2017). Insects are key elements of most terrestrial trophic interactions, provide essential ecosystem services of global relevance, and generate significant economic and aesthetic benefits with cultural value to human society (Wagner 2019). At the same time, insect species and populations are threatened by a variety of factors, the most crucial being habitat loss, fragmentation and deterioration, followed by climate change (Franco et al. 2006;Potts et al. 2010; Vanbergen and the Insect Pollinators Initiative 2013).
Species are affected differently by environmental change depending on their ecological tolerance and other species-specific life history traits, determining the likelihood of decline and local extinction (Cahill et al. 2013). In most cases, the extinction of local populations depends on the combination of the geographic location of the population within the distribution range, the ecological tolerance of the species and shifts in habitat conditions (Lawton 1994;Cahill et al. 2013). Hutchinson's concept of the ecological niche (Hutchinson 1957) divides the multidimensional parameter space of habitat variation roughly into two regions: (1) the region in which the number of local births exceeds that of local deaths (source populations), close to the optimum niche of the species, and a source of emigration; and (2) the region in which the number of deaths exceeds that of local births (sink populations), typically far from the environmental optimum of the species, where populations are maintained by immigration from source populations (Holt 1996). Consequently, the proximity to the optimal environmental space (centroid) is expected to correspond to higher genetic diversity, growth rates and population stability (Pulliam 1988;Vanderwal et al. 2009). Known as center-periphery hypothesis, this concept has provided a baseline for many studies of population dynamics and genetic variability at species distribution limits.
We know of only few studies that have explored the center-periphery hypothesis regarding genetic variability in relation to environmental suitability or niche centrality-with contrasting results. Diniz-Filho et al. (2009), for example, related genetic diversity to average ecological niche model (ENM) suitability scores derived from multiple correlative niche modeling algorithms. Lira-Noriega and Manthey (2014) focused on testing linear regressions of genetic diversity measures (e.g., allelic richness, nucleotide diversity) in comparison with Euclidean distances between the environmental centroid of the suitable area estimated by ENM, and assumed to be the fundamental niche of the species. The correlations were tested on 40 species, including plants, birds, mammals, worms and insects. No clear relationship was recovered in any study, with the majority of cases showing negative correlations and a few cases displaying positive correlations (Lira-Noriega and Manthey 2014). The mixed results might have been a product of the lack of a direct assessment of the species environmental centroid. Hence, a possible alternative to the approaches previously used is to characterize the ecological niche from a Grinnellian perspective-a set of environmental conditions that allow a species to maintain populations for long periods of time without immigration events . This concept has been successfully applied to study the relationship between geographical distances of populations to the niche centroid and their abundances (Yañez-Arenas et al. 2012;Urena-Aranda et al. 2015). However, no study has explored thus far the relationship of Grinnellian niche centroids with the geographic centroid and/or genetic diversity to explain patterns of population decline.
Due to well documented population declines within a huge part of its natural range, the speckled buzzing grasshopper Bryodemella tuberculata is well suited to investigate this relationship. This species occupies a wide global range across the northern Palearctic from Central Europe to eastern Asia (Bagachanova et al. 2011;Budrys et al. 2008;Budrys and Pakalniškis 2007;Fartmann et al. 2008;Reich 1991;Sergeev 1992;Srinivasan and Prabakar 2013). Up to the mid-20th century, B. tuberculata was still common in the European parts of its range, as documented by many specimens deposited in museum collections. It was frequently recorded from heath areas and river banks (Zacher 1919;Krauss 1883;Graber 1872). Over the course of the 20th century, many populations became extinct because their habitats (heath areas and unregulated river banks) were destroyed due to human demographic expansion, specifically a change in land-use due to intensified agriculture and shipping industry, and infrastructural developments (e.g. regulation of rivers; Laussmann et al. 2010;Bischoff 1997;Breunig and Thielmann 1992;Quinger 2014). Today, B. tuberculata has become a flagship species for the conservation of its threatened habitat types, mostly due to its large size for a grasshopper, its bright reddish hindwings and its characteristic buzzing flight noise. The most recent evaluation by the IUCN in 2016 (Zuna-Kratky et al. 2016) indicated that the populations in Denmark, Latvia, Poland and Switzerland are extinct (Berg et al. 2005;Budrys et al. 2008;Budrys and Pakalniškis 2007;Liana 2004;Maas et al. 2002;Monnerat et al. 2007), and the same probably applies to those in Italy (Massa et al. 2012). Furthermore, the relic populations of most other European countries have strongly declined despite frequent monitoring and legal protection (Binot-Hafke et al. 2011;Maas et al. 2011). In Central Europe, B. tuberculata is now restricted to small areas in Germany and Austria along the rivers Isar and Lech. In Northern Europe it is restricted to the island Øland in Sweden ( Fig. 1; Bieringer and Weißmair 2017;Pfeuffer 2004;Voith et al. 2016). At the same time, the species remains common in North and Central Asia, occupying almost all suitable areas in the region (pers. observed by the authors).
This recent and steep decline of B. tuberculata populations offers an optimal setting to test the center-periphery hypothesis, and in connection with analysis of the land-use change potentially providing a framework to understand the factors leading to its population decline. Hence, we calculated ellipsoid envelope models based on a large dataset of present and past occurrences of the species to characterize the climatic niche of B. tuberculata and test the center-periphery hypothesis by comparing relationships between genetic diversity, climatic niche, and geographic centrality between extinct and extant populations. We further evaluated landscape change in the European area of the distribution and evaluate several extinct and extant locations for their habitat structure and suitability. We hypothesize that (1) changing climatic conditions over the last 100 years are not the reason for the decline of the species, but rather (2) changes in land-use and land-cover are responsible for the decline and local extinction. Assuming that genetic diversity and proximity to the niche centroid are related to fitness and survivability (Reed and Frankham 2003;Leimu et al. 2006), we specifically seek to understand whether the relationship between those factors could provide any insight into the susceptibility of populations (and possibly species) to decline and extinction. Based on this, we hypothesize that (3) the genetic diversity in the Asian core populations is higher than in the European relict populations.
Due to the low dispersal capacity of B. tuberculata (Reinhardt et al. 2005) every single record was considered as a population. If several specimens from the same location were recorded, duplicates were removed. We considered any European population as extinct if it was listed as extinct in the relevant literature (Budrys and Pakalniškis 2007;Voith et al. 2016;Ruffo 2003;De Carlini 1889;Schmidt and Lilge 1997;Baur and Museum 2006;Bakker et al. 2015;Zuna-Kratky et al. 2017;Glowacinski and Nowacki 2006) or if no findings were documented after 1999 at the sampling location or within a 10 km buffer zone around the occurrence. All other populations from Europe and all known populations from Asia were considered extant. After data cleaning and trimming, 280 extant and 100 extinct locations were used to calculate the environmental distance to the centroid. All records used for modelling are supplied in Supplementary Information SI 3 and visualized in Fig. 1.

Environmental data
Environmental layers were obtained from Worldclim v. 1.4 (www. world clim. org; Hijmans et al. 2005) in 2.5 arc-min resolution. WorldClim is based on interpolations of weather station data (i.e., monthly precipitation and minimum and maximum temperatures) over the period . From the 19 variables available, we excluded four (mean temperature of wettest quarter, mean temperature of driest quarter, precipitation of warmest quarter, precipitation of coldest quarter) a priori due to known spatial artifacts between adjacent grid cells (Escobar et al. 2014;Campbell et al. 2015). Seeking to avoid bias regarding the combination of variables used to characterize the species niche centrality, we tested three distinct environmental sets: 'set 1' included all 15 variables; 'set 2' included only temperature variables (i.e., annual mean temperature; mean diurnal range; isothermality; temperature seasonality; max. temperature of warmest month; min. temperature of coldest month; temperature annual range; mean temperature of warmest quarter and mean temperature of coldest quarter); and 'set 3' included only precipitation variables (i.e., annual precipitation; precipitation of wettest month; precipitation of driest month; precipitation seasonality; precipitation of wettest quarter; precipitation of driest quarter). To avoid overfitting, overly dimensional environmental space and collinearity among variables, we performed principal component analysis (PCA) using the function kuenm_rpca in the package 'kuenm' (Cobos et al. 2019b) in R 3.6.3 (R Core team 2020). For each set, we retained as many components as necessary to explain > 95% of the total variance in the dataset for model calibration. As a result, set 1 was reduced to five principal components (PCs) summarizing the 15 bioclimatic variables, while set 2 and 3 were reduced to three PCs, summarizing temperature and precipitation variables (see SI 6).

Ellipsoid models
To characterize the environmental niche of B. tuberculata we created an ellipsoid envelope model representing the niche shape assumed when multiple dimensions are considered (Jiménez et al. 2019) using the 'ellipsenm' package (Cobos et al. 2019a). Based on these 1 3 models, we calculated the niche centroid for all collected extant records of the species and extracted the Mahalanobis distances between grid cells representing local environments of study populations and environmental conditions of the optimum (ellipsoid centroid). Models were calibrated using the 95% pairwise confidence region for the ellipsoid and were evaluated as candidate models using the function 'ellipsoid_calibration' of the 'ellipsenm' R package (Cobos et al. 2019a). Two distinct methods were employed to construct ellipsoid models: (1) 'covmat', which creates ellipsoids based on the centroid and a matrix of co-variances of the variables and (2) 'mve1', which generates an ellipsoid that reduces the volume contained in it without losing the data contained (i.e., minimum volume ellipsoid; Van Aelst and Rousseeuw 2009). Best model selection was based on statistical significance (partial ROC; Peterson et al. 2008); the proportion of testing data known to be in suitable areas and prediction of unsuitable areas was based on omission rates (E = 5%; Anderson et al. 2003) and prevalence. To calculate the partial ROC metric, we used 500 bootstrap iterations with 50% of testing data to be used in each bootstrapped process with 5% of testing data error in the data due to uncertainty. The calibration area (i.e., region accessible to the species; Barve et al. 2011) included Eurasia, except Southeast Asia, where the species is not found.
Final parameters were selected based on best evaluated models and used to create the final models using 10 replicates with bootstrapped subsamples of 75% of the data using the function 'ellipsoid_model' available in the 'ellipsenm' R package (Cobos et al. 2019a). The replicates were produced by excluding one occurrence record at a time. The confidence level of a pairwise confidence region for the ellipsoid was 99%, with 1% of the occurrence data considered as potential environmental outliers, not included in the ellipsoid envelope model for the species' ecological niche.

Environmental and geographic distances
The environmental distance of each record (= population) to the ellipsoid centroid (= ecological optimum) was calculated using the Mahalanobis distance metric, which is considered wel-suited for non-spherically symmetric distributions and deals robustly with variables with different scales, such as the bioclimatic variables we used for characterizing niche space (De Maesschalck et al. 2000;Hijmans 2020). We generated histograms showing the Mahalanobis distance (D 2 ) of all populations to the centroid.
To test whether the distribution of Mahalanobis distances of extinct populations to the centroid are significantly distinct from the distances recovered for extant populations, we performed a non-parametric Mood's Median test, which tests if two or more independent samples originate from populations with the same median. We used this to compare the shape and scale of distributions for all groups, based on the hypothesis that the Mahalanobis distances of extant populations are lower compared to those of extinct populations (p-value > 0.05). Analyses were performed in R. Distribution maps and Mahalanobis plots were created in QGIS v. 3.10.3 (QGIS Development Team 2020).

Land-cover change
Furthermore, we used the HILDA dataset v. 2 (HIstoric Land Dynamics Assessment 2gross land-cover changes; Fuchs et al. 2013Fuchs et al. , 2014Fuchs et al. , 2015 http:// www. geo-infor matie. nl/ fuchs 003/#) to display landscape change through time from 1900 to 2000, the main timeframe of local extinctions in Europe. Based on this, we calculated the percent land-cover change for forest, grassland, settlement and cropland area from 1900 to 2000 using the QGis plugin LecoS-Landscape Ecology Statistics v. 3.0.0 (Jung 2016).

Molecular analyses
For this study, we tested several genetic markers. Most gene fragments were highly conserved and did not show any variability; genomic data is currently not available to us and remains challenging because of the extremely large genome sizes in Orthoptera and specifically Oedipodinae (Husemann et al. 2021). In turn, DNA barcode data is known to provide meaningful estimates of genetic diversity, and additional genetic data of this marker was available from other studies in BOLD and GenBank. Hence, we decided to use a fragment of the Cytochrome Oxidase 1 for the genetic comparison of the European and Asian populations. We compiled COI sequences of 34 specimens of B. tuberculata from populations across its global range (Table 1). DNA was extracted from muscle tissue of a hind femur of each specimen with a high salt extraction protocol (Paxton et al. 1996). We amplified the barcoding fragment of the mitochondrial Cytochrome C Oxidase subunit I (COI) gene using the primers provided by Pan et al. (2006). For PCR reactions we used DreamTaq DNA polymerase (Thermo Fischer Scientific, Schwerte, Germany) following the manufacturer's protocol. Cycling conditions are shown in SI 4. PCR products were visualized using 1 µl of EZ-Vision Two DNA Dye (Amresco, Darmstadt, Germany); samples were run on 1% agarose gels. We purified successfully amplified samples using the Thermo Fisher protocol for PCR product clean-up prior sequencing (Werle et al. 1994). The samples were sequenced at Macrogen Europe (Macrogen, Amsterdam, Netherlands) using the EZ-Seq service. Further, sequences from earlier projects of MH and OH (Husemann et al. 2012;Hawlitschek et al. 2017) were obtained from BOLD (see Table 1).

Analysis of genetic diversity
The MUSCLE algorithm (Edgar 2004) as implemented in Geneious v. 10.0.9 (http:// www. genei ous. com, Kearse et al. 2012) was used to edit, trim and align the sequences, resulting in three datasets for molecular analyses. The first dataset included all 34 sequences (total alignment length 609 bp). The second dataset only comprised samples from Europe (Germany N = 15; Sweden N = 1; Lithuania N = 1) with a total alignment length of 609 bp. A third dataset was generated comprising the samples from Central Asia (Russia N = 9; Mongolia N = 5; China N = 3; Korea N = 1) with a total alignment length of 609 bp. Geneious was used to check for internal stop codons (none were present). All sequences were blasted in NCBI GenBank and BOLD systems v. 3 (Ratnasingham and Hebert 2007) to exclude  (Excoffier et al. 2005). We used PopArt (Leigh and Bryant 2015) to construct a TCS haplotype network (Clement et al. 2000) of the full dataset without outgroup and with seven geographical traits (Germany, Sweden, Lithuania, Russia, Mongolia, China. Korea).
To test for demographic stability we conducted Bayesian Skyline analyses using BEAST v. 2.5.1, based on the reduced datasets of Europe and Central Asia. HKY was used as substitution model. A relaxed log-normal clock using the substitution rate provided by Papadopoulou et al. (2010) was implemented as clock model to calculate a Coalescent Bayesian Skyline plot; the analysis was run for 100 million generations sampling every 1,000 generations. We ran BEAST v. 2.5.1 with automatic thread pool size and BEAGLE library (Ayres et al. 2012). Convergence was assessed with TRACER v. 1.7.1 to avoid insufficient sampling. TRACER was also used to calculate the Bayesian Skyline reconstruction (Bouckaert et al. 2014) discarding 10% of samples as burn-in.

Species distribution
We assembled 280 extant distribution records from Austria (28) (2), Latvia (5), Poland (31) and Switzerland (5). Location files including coordinate data are provided as supplementary information SI 3. Figure 1 shows all records we used for modelling. The northernmost location was found in Siberia (Russia) at 67°N; the southernmost location was in Tibet (China) at 28°N. Most of the occurrences were recorded between 57°N and 45°N. The longitudinal extension stretches across most of Eurasia from 152°E in Russia to 12°E in Germany over a range of around 7500 km.

Ellipsoid model
The best fitting method to construct the climatic ellipsoids was 'covmat'; mean AUC, p-value of partial ROC and omission rates were significantly better than random expectations (p-value < 0.05; Table 2). The complete report of ellipsoid characteristics (e.g., centroid, covariance matrix, semi-axes length, etc.) can be found in SI 5.

Environmental and geographic distances
The geographic prediction of Mahalanobis distance (D 2 ) in the calibration area is shown in Fig. 2 and SI 6 and SI 7. The comparison of D 2 on all environmental sets shows all extinct populations, except for one in northern Austria (Ried im Innkreis: 48.222461°N 13.484835°E) to be as close to the climatic centroid as extant populations (Table 3). The Table 2 Calibration and evaluation of ellipsoid models used to characterize the niche of Bryodemella tuberculata Bold rows highlight the method selected to create final models; we provide the evaluation metrics (mean AUC, p-value partial ROC, omission rate), valid iterations and mean prevalence calculated in environmental ('Prevalence on E-space') and geographical space ('Prevalence on G-space')   Fig. 3). Moreover, the average distances to the centroid were found to be higher for the extinct and extant European locations (distance: 38.00-66.28, Ø = 59.53 (± 6.08)) than for Asian populations, which were geographically closer to the centroid (distance: 3.86-80.42, Ø = 30.39 (± 18.27); overall p-value < 0.001; Fig. 3).

Land-cover change
Using the HILDA dataset we were able to track changes in land-cover in Central Europe from 1900 to 2000 (SI 8): a massive change from cropland (1900: 33.9% to 2000: 29.3%)

Analysis of genetic diversity
We included 34 COI sequences of B. tuberculata in the analyses. These represented a total of 28 haplotypes with a haplotype diversity of 0.978 (± 0.015) and a nucleotide diversity of 0.031 (± 0.015). The pairwise distance between central (Asia) and peripheral (Europe) populations was moderate, but significant (Φ ST of 0.104; p-value < 0.001). The haplotype network (Fig. 4) shows a close relationship of all haplotypes with few evolutionary steps between them. Two main groups are recognizable, albeit with only few separating mutations. These two groups do not reflect a split between European and Asian samples, as European samples are grouping with sequences from Central Asia. The samples from China are more closely related to those from Korea and to most of the Mongolian samples representing the first group; the other group is composed of the Russian, Mongolian and European samples (Germany, Sweden, Lithuania). Furthermore, the Russian, German, Lithuanian and one Mongolian sample share the same haplotype.
Bayesian Skyline plots tracing female effective population size through time were calculated for European and Central Asian specimens separately (Fig. 5). The analysis showed a clear trend of population decline during the last 1000 years for the European populations. The same analysis performed for the Central Asian samples yielded a 10 times higher estimate of female effective population size and demographic stability through time. Colors represent countries, mutational steps between haplotypes are shown as marks across connection lines

Discussion
In this study, we investigated potential reasons for the decline of the European populations of the formerly widely distributed band-winged grasshopper B. tuberculata using ecological niche modelling in combination with a basic population genetic assessment. Specifically, we aimed to test if (1) changing climatic conditions over the last 100 years, or (2) changes in land-use and land-cover are responsible for the decline and local extinction. For this we assessed the habitat requirements and the decline of the species and modelled the environmental niche. Furthermore, we tested if (3) the genetic diversity in the Asian core populations is higher than in the European relic populations, which may support the centerperiphery hypothesis. We discuss the results in regards of conservation management of the species.

Climate models reveal suitable conditions in extinct locations
We used an ellipsoid envelope approach to estimate the climatic niche of B. tuberculata and measure the proximity of its populations to their optimum. All sets recovered low distances of the extinct populations to the climatic niche centroid, and none recovered any difference between extinct and extant populations (Figs. 2, 3; SI 6-7). Set 1 offered higher variation in distances to the centroid (range: 0.417-32.611), possibly due to the higher variation and multi-dimensionality of the principal components included (i.e., 15 variables).
While most extinct populations were located in areas with lower Mahalanobis distance (D 2 ) to the centroid (i.e., high suitability), two extinct locations with higher values were recorded (D 2 = 22.306 and 17.142; Fig. 2, Table 3). These locations represent a population close to the German-Austrian border in the Berchtesgaden national park (D 2 = 17.142) and a population in Austria close to Ried im Innkreis (D 2 = 22.306). The higher D 2 at these two locations could have several reasons, e.g. imprecise locality data, errors in georeferencing, individuals dispersed from neighboring locations (e.g. migrating individuals; Barton and Hewitt 1982), artefacts of the bioclimatic layers used (e.g. Marchi et al. 2019;Poggio et al. 2018), or even just a high tolerance regarding abiotic conditions (e.g. Preston and Johnson 2020). The last reason is supported by data from Asia with similar D 2 range. Since our models predict areas situated in the distribution gap between the European and Asian populations as climatically suitable, the destruction of natural habitats is a likely reason for the loss of stepping stones and for the extinction of any connecting populations in the area.

Habitat requirements and historical decline of B. tuberculata
Bryodemella tuberculata is a species of open landscapes, preferring habitats with sparse vegetation and open ground. According to our observations, all Russian collection sites in the Altay-Sayan Mountains, the Kulunda Plain (the south-eastern part of West Siberian Plain between the Irtysh and Ob Rivers) and within the south-eastern Ural Mountains (pers. obs. MGS; for further details, see SI 9) were mainly composed of dry meadows and steppes with sparse vegetation (~ 10-30 cm height) and spots of open ground (generally, the vegetation cover was less than 80%). The ground is mostly sandy to rocky, with smaller pebbles. In most cases occurrence sites were close to water systems. The species enters transformed ecosystems, such as the steppe variants with moderate livestock grazing and the dry hayfields. This type of habitat was found at almost all visited collection sites in Central Asia and Europe (for habitat images see Supplementary Information SI 10). In Mongolia, the species is found in low abundance, but in almost all habitats with open spaces and sparse vegetation, mainly covered by shorter plant species up to 30 cm and at the edges to forests (pers. observ., LSD 2015-2019; Dey et al. in press; Supplementary Information SI 10). Pastures for grazing represent the major agricultural use in the country, hence, the landscape is characterized by bare steppe habitats with low vegetation providing good conditions for B. tuberculata.
All extant populations of Central Europe (Southern Germany, Austria) inhabit alpine river valleys with unmodified or little modified flood regimes; specifically zones of alluvial gravel beds with sparse or without any vegetation. Occasionally, individuals may be observed in adjacent structures similar to those described for the Asian populations (Landmann 2017). Since the middle of the 19th century, populations of B. tuberculata in Europe have been declining (e.g. Reich 1991;Zuna-Kratky et al. 2017). Despite the availability of a large database, we are unable to make any statements about the exact timing of extinction of B. tuberculata at any particular location, but we provide the last recorded findings for several regions ( Fig. 1; Supplementary Information SI 11). Altogether, our data suggest that most Central European extinction events took place between the 1920's and 1960's. We therefore suggest that each local extinction event in Central Europe may be closely connected with local landscape changes. Historically, at least in the mountainous regions, the naturally occurring annual floods following snowmelt removed most vegetation growing in the riparian habitats. Anthropogenic reductions of the flood volume, mostly due to the construction of hydroelectric dams, led to the establishment of shrub and tree vegetation in these habitats and a decline of the areas suitable for B. tuberculata (Reich et al. 2008;Juszczyk et al. 2020). A variety of restoration measures have been put in place over the last decades, including a reduced deduction of water for hydroelectric energy generation, mechanical removal of vegetation, and education of the public; the mid-and long-term efficiency of these measures remains to be demonstrated (Juszczyk et al. 2020). A different type of change to the habitat structure was recorded in Denmark. There, most heathlands in the vicinity of Abild (an extinct location) have been converted to agricultural areas or were overgrown by shrubs (pers. observ. LSD; SI 10). Although the agricultural land-use decreased from 74% in 1915 to 61% in 2015, afforestation and urban expansion still led to a decline of heath areas previously typical for many parts of Denmark (Pedersen and Møllenberg 2017). The species has been listed as extinct in Denmark since the 1950s. Similarly, in Northern Germany, restoration of some old heath landscapes has been performed. Based on literature and museum surveys, we were able to find several records from the Lüneburger Heide heathland areas, for which the gradual conversion to arable land and forests at the end of the 19th century is historically documented (Koopmann 2000). This change in land-use was mainly triggered by a decrease in the local production of wool and honey due to external competitors, making the maintenance of large pastures uneconomic (Naturpark-Lueneburger-Heide.de; Koopmann 2000). During the early 20th century, the first areas were restored to the original heath systems. Nowadays, the Lüneburger Heide is one of the largest heath areas in Northern and Central Germany. Although most heath areas in Central Europe declined through time, some of these cultural landscapes are still intact. Many of them are now under military use and kept open by military activity, which practically equals to legal protection of the landscape. As a result, B. tuberculata was rediscovered on 03 August 2008 in the Pabradė military training area in Lithuania, where it had been assessed as possibly extinct until 2008 (Budrys et al. 2008;Budrys and Pakalniškis 2007). This was the first record since the first half of the 20th century (Budrys and Pakalniškis 2007); yet, the status of the population remains unknown and it will have to be monitored to check its population establishment. In Northern Europe the population on Øland (Sweden) shows similar habitat preferences for barren land. The implementation of specific conservation measures appears to keep the population relatively stable.
Bryodemella tuberculata needs larger areas with open habitat, often with frequent natural disturbance. Anthropogenic changes of the European landscape in the last centuries (Plieninger et al. 2016;Hersperger and Bürgi 2009;Antrop 2004) likely led to the current patchy distribution of the species and will probably cause the extinction of further populations in the European range. The time series of land-cover plots (Supplementary Information SI 8) supports land-use change as a driver of the decline of B. tuberculata as it displays a rapid change from grassland and agricultural areas to settlement and reforestation in the middle of the 20th century and still ongoing. This change from open areas to more closed habitat types, and an increase of human pressure due to high-intensity land-use, may have led to the extinction of B. tuberculata at many Central European locations. Even though some habitats of extinct populations have now been restored, they are highly fragmented, and the dispersal capabilities of B. tuberculata appear insufficient to allow for quick colonization of these habitats without support. Similar scenarios of decline have been described for many other European species, e.g., carabid beetles with different ecological preferences in Belgium, Denmark and the Netherlands (Kotze and O'Hara 2003); common and widespread butterflies in the Netherlands (Van Dyck et al. 2009); or the ground nesting Black Grouse populations in Lower Saxony, Germany (Ludwig et al. 2009). In turn, the availability of large stretches of natural habitats suitable for B. tuberculata currently remains much higher in the Asian range. However, vast areas of habitat are also being destroyed in Russia, especially in the European part of the country, due to urbanization, agriculture and mining (Smelansky and Tishkov 2012). This process is much slower due to lower human activity in these areas and simply a larger area, but can also be expected to fragment populations and impact genetic diversity in the near future.

Distribution of genetic diversity and historical demography
While our genetic data is limited, also due to the rareness of the species, it suggests that the majority of genetic diversity is found in Central Asia. The European populations show lower diversity and are nested within the Asian populations. Hence, based on our COI data, European populations only represent a subset of the species' genetic diversity, as expected for relic populations, which may have gone through a population bottleneck (e.g. Chen et al. 2016;Gaublomme et al. 2013;Hájková et al. 2007;Lucchini et al. 2004). We did not find any evidence supporting the status of these populations as distinct genetic lineages. This may suggest a previously more continuous, potentially panmictic distribution across much of Europe and Asia and a relatively recent fragmentation and decline in Central Europe, as also suggested by our historical distribution data.
The genetic patterns do not support the center-periphery hypothesis, i.e., low genetic variability in areas of low environmental suitability, as low genetic variability was also detected in areas of high environmental suitability. However, due to strong contemporary restrictions on migration, gene flow is likely completely interrupted and the European populations are likely under a drift regime, rendering them prone to extinction. On the other hand, Bayesian Skyline plots provide some evidence that supports the center-periphery hypothesis driven by geographic distance. Our results show that female effective population size is larger in Central Asia and also has been more stable in the last 10 ka (Fig. 5), while the European populations have smaller estimated population sizes and are in decline. The results we obtained match other taxa with similar distributions, e.g., the leaf-beetle Cheilotoma musciformis, which has some relic populations in Poland. These populations are declining, probably because habitat destruction caused a disjunction from the main distribution, resulting in an observed genetic bottleneck (Kajtoch et al. 2016). Such patterns of loss of genetic diversity associated with habitat destruction, which have led to the disjunction of populations, can also be seen in the decline of the Danish population of the otter Lutra lutra (Pertoldi et al. 2001), among others (e.g. Cremene et al. 2005;Kotze and O'Hara 2003).
Insufficient genetic samples and lack of more fine-scaled genome-wide markers inhibited the direct assessment of a relationship between genetic diversity and niche centrality. Nevertheless, our results support the findings of Lira-Noriega and Manthey (2014), who equally failed to detect any clear relationship between genetic diversity, and geographic and environmental distance to the geographic centroid. However, they did find a strong relationship between genetic diversity and climatic niche centrality. In some taxa, Lira-Noriega and Manthey (2014) describe a negative tendency for the relationship of genetic diversity and geographical distance to the source population, reflecting environmental impact on the population dynamics, rather than a fundamental ecological relationship. In case of B. tuberculata, we suggest a similar effect of habitat destruction in the periphery of the distribution supported by decreasing female effective population size in the European populations (geographically distant from the main distribution; source) in contrast to the stable demography in the Asian populations (geographically close to the source). A more detailed study of other habitat characteristics, especially land-use patterns, may shed light on the reasons for the local decline and could help to manage the conservation activities to save the remaining European populations from decline, which may also facilitate the recolonization of former locations.

Aspects of conservation management
Local conservation efforts for species with isolated and declining populations may often be seen as unsustainable, because these populations may be threatened by the effects of climate change, which cannot be counteracted by localized efforts (e.g. Sinervo et al. 2010, Malcom et al. 2006. Our results show that extinct or threatened European populations of B. tuberculata are in areas of optimal climatic conditions for this species, whereas some Asian populations appear to be thriving also in areas of poor climatic suitability. This suggests that B. tuberculata is most likely resilient to the purely climatic effects of global warming in the Central European part of its range. If our results hold true, this means that threatened populations may survive the next decades of changing climate as long as their habitats are locally protected, and that these local efforts can be considered sustainable. This includes not only threatened extant populations, but also means that reintroduction efforts into the localities of extinct populations may be a promising option as long as the habitat structures necessary for B. tuberculata have been restored and long-term management measures are in place. We suggest that monitoring of extant and potentially translocated populations should be accompanied by genetic studies using deeper populationlevel sequencing methods, such as microsatellites or RAD sequencing. Once a high-quality whole genome becomes available, genome resequencing may become a viable option.
The speckled buzzing grasshopper is just one species of an assembly of taxa sharing a habitat that is highly threatened by anthropogenic land-cover change. We believe that our study reinforces that status of B. tuberculata as a flagship and umbrella species for this organism community, and we hope that our results help making a case for increased efforts to their protection. Data availability All genetic data will be available in GenBank. All distribution data will be available in the