Priority conservation areas for protected saproxylic beetles in Romania under current and future climate scenarios



Climate change threatens species and ecosystems globally, including forest ecosystems that support rich invertebrate diversity. Saproxylic beetles, that depend on old-growth trees and deadwood, are facing increasing pressure. Consequently, conserving these beetles has become a priority for EU Member States. We developed ensemble species distribution models for five saproxylic beetles for current and three future time horizons under two Shared Socioeconomic Pathways and two Global Circulation Models. We used a systematic conservation planning approach to assess the effectiveness and resilience to climate change of the Romanian Natura 2000 network for saproxylic beetles while identifying areas for prospective protected area expansion to meet EU conservation targets. Our study revealed that under all scenarios and time horizons, the saproxylic beetles may lose over 80% of their suitable habitat and restrict their distribution to higher elevations. According to the conservation prioritization analysis, we found that, when considering 30% of the landscape as being protected, an average of 85% of species distribution is retained within priority areas overlapping the Carpathian Mountains, while for the current protected area coverage (18% of Romania’s terrestrial area), the existing Natura 2000 network does not perform satisfactorily, with only ~ 30% of the saproxylic species distributions falling within the network. Our results corroborate previous findings on saproxylic beetle range shifts and contractions due to climate change. Furthermore, our findings question the effectiveness of the current Natura 2000 network, as it is currently inadequate for protecting these species. To achieve the goals of the EU Biodiversity Strategy 2030 of protecting at least 30% of the EU’s territory, we advocate the expansion of the Natura 2000 sites to future suitable saproxylic beetle habitats.


Human-induced climate change acts synergistically with other threats to globally imperil habitats and species (Mantyka-Pringle et al. 2012; Harvey et al. 2023). The most unmistakable consequence of climate change is the increase in global temperatures when compared to pre-industrial levels, with Europe being the fastest-warming continent (Kjellström et al. 2018). Under these conditions, habitats important for saproxylic beetles, such as old-growth forests, are expected to experience severe pressure and threats (La Porta et al. 2008; European Commission 2021). Because of the importance of forests for human well-being and viable biodiversity, the European Union (EU) not only set an overall goal of protecting a minimum of 30% of the EU’s land area (of which one-third should be strictly protected), but also considered old-growth forests as a priority for stricter protection (European Commission 2021).

Saproxylic beetles are deadwood specialists and keystone species in forest ecosystems. The diversity of saproxylic beetles is an indicator of a healthy forest ecosystems (Jansson et al. 2009; Mazzei et al. 2018), and thus, the conservation of beetles is a priority for EU Member States. Yet, insects, in general, are declining at an alarming rate in both altered and intact landscapes (Seibold et al. 2015; Wagner 2020). In Europe, species associated with deadwood habitats are among the most threatened taxa due to intensive forest management practices and habitat fragmentation (Nieto and Alexander 2010), alongside broader-acting threats such as climate change. Because of the affinity of saproxylic beetles for deadwood, their diversity is positively correlated with the amount and diversity of deadwood (Lassauce et al. 2011; Lachat et al. 2012; Seibold et al. 2018); therefore, a decline in old-growth forest habitats has led to range contractions for several species (Seibold et al. 2015). In contrast, research has indicated that saproxylic beetles benefit from increased temperatures (Müller et al. 2015) and that beetles differ in their habitat requirements for deadwood based on life history stages and latitude (Winiger et al. 2023). In Europe, 21 saproxylic beetle species are listed in the EU Habitats Directive, including Rosalia alpina, Cerambyx cerdo, and Lucanus cervus (Cálix et al. 2018). The EU has developed strategies for the recovery of saproxylic beetles, such as the EU Biodiversity Strategy for 2030 (European Commission 2020), yet range shifts and forest structure changes have not yet been accounted for in conservation planning.

Among European countries, Romania embraces the most contiguous forested areas, which include the Carpathian Mountains, a biodiversity hotspot for saproxylic beetles (Munteanu et al. 2022; Stanciu et al. 2023). Carpathian forests shelter viable populations of several iconic saproxylic beetles listed in Annex II of Habitats Directive (Directive/92/43/EEC 1992): Rosalia alpina (Coleoptera: Cerambycidae), Lucanus cervus (Coleoptera: Lucanidae), Cerambyx cerdo (Coleoptera: Cerambycidae), Osmoderma eremita (Coleoptera: Scarabaeidae) and Morimus funereus (Coleoptera: Cerambycidae). Past forestry management practices in the Carpathians, that include selective logging and removal of deadwood and old-growth trees, degraded the forest structure and led to a local decrease in species abundance (Prunar et al. 2013; Olenici and Fodor 2021). Complicating matters, Romanian forest managers are faced with conflicting mandates regarding these species, i.e., either considering them as pest species and applying eradication methods to lower their populations or protecting them as threatened and endangered species per the EU environmental mandates (Brodie et al. 2019). As a result, forest managers have neither the incentives nor technical information on the saproxylic beetle community to promote concrete conservation actions. Thus, the most important tool in protecting these species is the inclusion of their habitats in the Natura 2000 network, one of the most extensive networks of conservation areas in the world, which has been created to operationalize EU Birds (Directive 2009/147/EC 2009) and Habitats Directives (Directive/92/43/EEC 2013). However, the designation of Natura 2000 sites often lacks clear, quantifiable conservation objectives or planning, as highlighted in other studies (Iojă et al. 2010; Kukkala et al. 2016a; Miu et al. 2020; Cazzolla Gatti et al. 2023). Frequently, Natura 2000 designations result from the pursuit of area-based targets established by the European Union for country-level protection, with a specific goal of protecting 30% of each EU country by 2030 (European Commission 2020).

Systematic conservation planning is an effective approach for identifying a network of protected areas that align with EU targets (Margules and Pressey 2000; Wintle et al. 2019). This framework optimizes conservation benefits while minimizing adverse influences on other resources. Spatial conservation prioritization, as part of systematic conservation planning, typically employs algorithms that account for complementarity and representativeness of species and communities to identify areas that complement each other to prevent redundant conservation efforts (Mikusiński et al. 2007; Kukkala et al. 2016b; Kujala et al. 2018). In Europe, the overall conclusion of studies that assessed the distribution of saproxylic insects and the effectiveness and representativity of Natura 2000 has shown that planning is suboptimal (D’Amen et al. 2013; Zehetmair et al. 2015; Bosso et al. 2018). For example, Bosso et al. (2013), indicates that the Natura 2000 network protects less than 56% of Rosalia alpina’s suitable habitat in Italy, while Lachat et al. (2013), found that only 11% of its suitable habitat is protected in Switzerland. In their study, Bosso et al. (2018), found that only 25% of Rosalia alpina’s potential distribution is covered by Natura 2000 sites in France, and only 35% of its suitable areas are protected in Austria. For Romania, studies assessing the coverage and effectiveness of Natura 2000 in protecting species listed in Annex II of the Habitats Directive are scarce (Iojă et al. 2010; Popescu et al. 2013; Miu et al. 2018, 2020), mostly due to species occurrence data limitations. Specifically, for saproxylic insects, studies have focused on evaluating the diversity of insects in smaller areas (Stan and Nitzu 2013; Bărbuceanu et al. 2015; Manu et al. 2016, 2017, 2019; Stan et al. 2016; Maican et al. 2019; Brodie et al. 2019). Thus, there is a pressing demand for national-level systematic conservation approaches to identify the representation of these taxa in the current protected area network and to highlight gaps and additional areas for protection to achieve national and European protection targets. Such approaches would also benefit forest biodiversity conservation because saproxylic beetles are considered habitat engineers (Buse et al. 2008) and can function as surrogates for other taxa (Ranius 2002; Holland 2007; Foit et al. 2016).

The aim of our study was to identify conservation priorities for five listed saproxylic beetle species (Rosalia alpina, Lucans cervus, Cerambyx cerdo, Osmoderma eremita, and Morimus funereus) in Romania under current and future climate change scenarios. First, we developed species distribution models for the current and three future horizons (2021–2040, 2041–2060, and 2061–2080) and two Shared Socioeconomic Pathways. We then used a systematic conservation planning approach, in combination with forest cover information, to identify the representation of these species in the current protected area system and to identify additional priorities to meet EU conservation targets. Specifically, our objectives are as follows: (1) to map the distribution of five target species using an ensemble modeling approach; (2) to assess the changes in the distribution of saproxylic beetles due to climate change; (3) to assess the resilience to climate change of Romanian Natura 2000 network covering saproxylic beetles and (4) to identify future areas for protected area expansion. Overall, this research aims to aid the current efforts of Romanian and EU decision-makers to identify optimal areas for protected area expansion and to develop management plans that include forest conservation strategies for saproxylic beetles.

Materials and methods

Study area, species, and occurrence data

Romania is a biodiversity hotspot in Europe, overlapping five European biogeographical regions, i.e., Alpine, Continental, Pannonian, Steppic, and Black Sea (Rozylowicz et al. 2019). Over 27% of its territory is covered by forests (Munteanu et al. 2016), with more than 3.5% comprising old-growth forests (Knorn et al. 2013). Forest habitats are dominated by broadleaved species (70% of the forest habitats, mostly Quercus ssp. and Fagus sylvatica as dominant species) and deciduous species (30% of the forest habitats, mainly Picea abies and Abies alba as dominant species) (Veen et al. 2010). Due to the large extent of old-growth and less intensively managed forests, Romania is also considered a hotspot of saproxylic beetle biodiversity (Nieto and Alexander 2010). In Romania, over 20 saproxylic beetles are protected by the Habitats Directive (Gîdei and Popescu 2012, 2014; Directive/92/43/EEC 2013; Fusu et al. 2015), although the saproxylic beetle diversity is more than one order of magnitude higher even for specific areas (Ilic and Curčić 2013; Brodie et al. 2019). Natura 2000 protected areas created for Habitats Directive species and habitats include 425 Sites of Community Interest, covering 40,500 km2, i.e., 17% of Romania (European Environment Agency 2021).

The five study species, the alpine longicorn Rosalia alpina Linnaeus, 1758, the Morimus longicorn Morimus funereus Mulsant, 1863, the great Capricorn beetle Cerambyx cerdo Linnaeus, 1758 (Cerambycidae), the hermit beetle Osmoderma barnabita Motschulsky 1845, part of the Osmoderma eremita complex (Scarabaeidae), and the stag beetle Lucanus cervus (Linnaeus 1758) (Lucanidae) are listed in Annex II of the EU Habitat Directive (Directive 92/43/EEC, 1992). These species have a relatively extensive range in Romania compared with other saproxylic species listed in Annex II of the EU Habitat Directive (e.g., Cucujus cinnaberinus, Buprestis splendens, Pseudogaurotina excellens, Rhysodes sulcatus) and they may act as umbrella species for other taxa that are data-deficient and have lesser distributions.

Saproxylic beetle species occurrences were retrieved from (i) public biodiversity databases (Global Biodiversity Information Facility ( 2023), (ii) peer-reviewed articles and technical reports (supplementary file S1), and (iii) citizen science data from social media (Facebook entomology groups, e.g., Insects of Romania and Europe) for the period 2005–2022. Following Marcer et al. (2022), we discarded occurrence data with more than 5 km uncertainty in the GBIF database. To further minimize sampling and spatial bias, we applied a thinning process to limit the data to one record per grid cell, followed by spatial thinning using the spThin R package to remove clustered occurrence records within a 2 km radius (Boria et al. 2014). The spThin package uses a randomization approach and returns a dataset with the maximum number of records for a given thinning distance (Aiello-Lammens et al. 2015). The final saproxylic beetle occurrence database contained 530 occurrence records: 60 for Cerambyx cerdo, 190 for Lucanus cervus, 118 for Morimus funereus, 38 for Osmoderma eremita, and 124 for Rosalia alpina. The grid used in the thinning process was created at the same resolution as the WorldClim database, i.e., ~ 1 km2, and was used to resample all datasets used in the study, i.e., occurrence records, environmental data, and protected areas dataset.

Environmental data for species distribution modeling

To model the current distribution of the five saproxylic species in Romania, we used 1970–2000 WorldClim 2.1 database at 30 s spatial resolution, i.e., ~ 1 km2 (Hijmans et al. 2005; Fick and Hijmans 2017). We extracted 19 bioclimatic variables, and, to avoid overfitting, we selected for analysis only the variables with a Pearson pairwise correlation <|0.75| (Dormann et al. 2013). The final set of variables used as environmental predictors includes six variables: Isothermality (BIO3) in percent, Temperature Annual Range (BIO7) in °C, Mean Temperature of Driest Quarter (BIO9) in °C, and Precipitation Seasonality (BIO15), Precipitation of Warmest Quarter (BIO18), and Precipitation of Coldest Quarter (BIO19), each measured in mm.

To infer potential changes in species distribution due to climate changes, we used two IPCC emission scenarios (Shared Socioeconomic Pathways, SPP): SSP1-2.6 sustainability pathway (i.e., best-case scenario) and SSP5-8.5 fossil-fueled development (i.e., worst-case scenario) from two general circulation models (HadGEM3-GC31-LL and MIROC6) (Hideo et al. 2019; Ridley et al. 2019). HadGEM3 and MIROC6 were selected because they offer diverse climate projections and help cover a wide range of probable future climates that incorporate uncertainty in predictions (Thuiller et al. 2019). For each general circulation model under an emission scenario, we selected three time horizons (2021–2040, 2041–2060, and 2061–2080), resulting in 12 modeling cases (Table 1). For example, “Had worst-case 2021–2024” indicates a climate model prediction from HadGEM3-GC31-LL for the time frame 2021–2024, based on the high-emission scenario SSP5-8.5 (i.e., worst-case scenario)

Saproxylic beetles distribution modeling

The current potential distribution of the five saproxylic beetles and predicted future changes were modeled using the BIOMOD2 R package for ensemble modeling, which consists of running a group of algorithms simultaneously (Thuiller et al. 2014). For each species, we fitted an ensemble SDM based on five modeling techniques: two regression methods, generalized linear model (GLM) and generalized additive model (GAM), two machine learning methods Random forests (RF) and Generalized boosted model (GBM), and maximum entropy (MAXENT). The use of five algorithms allows for a comprehensive and robust approach to predicting species distributions (Thuiller et al. 2014). Because we lacked true absence data for each species, we created pseudo-absences datasets. For each species, we generated 10 sets of random pseudo-absence records outside a buffer of 20 km from the presence points. The ratio between pseudo-absence and presence records was 3:1. This ratio was chosen to ensure robust modeling, as a balanced or slightly biased dataset toward absences can improve model performance by reducing overfitting and increasing the model’s ability to discriminate between presence and absence locations (Barbet-Massin et al. 2012).

Species distribution models were calibrated using 80% random samples from occurrence data, whereas the model performance was assessed using the remaining 20% of data (Gholamy et al. 2018). We evaluated the results of the SDM ensemble models using the area under the curve (AUC) of the receiver operating characteristic (ROC) and the true skill statistic (TSS) (Allouche et al. 2006). We also evaluated the contribution of dependent variables in predicting the species range for individual and ensemble models. We transformed the probability of occurrence for the models into binary present/absent format using the TSS cut-off value calculated by BIOMOD2 (Hao et al. 2019).

For each species, we created 13 ensemble species distribution models, one for the current climate and 12 predictions for the 2021–2040, 2041–2060, and 2061–2080 horizons (general circulation model × emission scenario × time horizon, see Table 1). To develop the ensemble models, we implemented the mean ensemble modeling algorithm (Hao et al. 2019), incorporating only individual models with a TSS > 0.4.

Priority areas for conservation of saproxylic beetles

Current and 2041–2060 future species distribution data were further used to select priority areas for conservation of saproxylic beetles. Spatial conservation prioritization was performed using the Zonation 5 app, a decision-support tool for spatial conservation planning (Moilanen 2022; Moilanen et al. 2022). Zonation produces a priority ranking by iteratively removing grid cells with the lowest total marginal loss of conservation value while accounting for the total and remaining distributions of protected saproxylic beetles. It produces a uniform distribution hierarchical ranking of the landscape from the highest (1) to the lowest (0) conservation value (Kujala et al. 2013; Moilanen 2022).

Spatial conservation prioritization was implemented using the mean ensemble occurrence probability for each species. We accounted for uncertainty in ensemble model predictions by subtracting the standard deviation of the ensemble model predictions from the mean ensemble values (Moilanen 2022). We included species distribution under the current climate scenario and four species distribution predictions for 2041–2060 timeline (two GCMs × two emission scenarios, see Table 1). We used only one time horizon because the earlier (2021–2040) time horizon incorporates our baseline data, and the later (2061–2080) time horizon is associated with higher levels of uncertainty and unpredictability in climate forecasts (Lee et al. 2023). For each scenario, we ran two prioritization analyses: (1) considering the Natura 2000 network (Sites of Community Interest and Special Areas for Conservation) as areas with the highest conservation value (i.e., constrained prioritization using protected areas as a hierarchical mask) and (2) selecting a priority area for conservation irrespective of the Natura 2000 network (i.e., unconstrained prioritization). The former identifies conservation priorities that complement the current protected area network (i.e., best areas to expand the network to accommodate the highest value area for beetle conservation). The latter identifies the areas with highest conservation value irrespective of their protection status; this prioritization can then be used to identify the representation of beetle species in the current protected area network and highlight gaps in protection.

As a prioritization algorithm, we used the Core Area Zonation type marginal loss rule, CAZ2 algorithm (Moilanen et al. 2022). CAZ2 maintains a relatively high average coverage of features while not significantly compromising the performance of the worst-performing features, making it a more suitable approach for achieving the actual conservation goals (Moilanen 2022; Moilanen et al. 2022). We also used condition with renormalization, an additional analysis option that represents information about local habitat requirements and their influence on biodiversity features. As a raster layer for habitat requirements, we used Corine Land Cover 2018 data (European Environment Agency 2019) to extract two main types of habitats: (1) all types of forests, and (2) all non-forested areas, such as places where agriculture and forestry mix, natural grasslands, and shrubs. The reasoning behind using forest as a condition layer is twofold: (1) our focal species are forest specialists and (2) the amount of forest cover is not predicted to decrease over the next four decades (Kucsicsa et al. 2020), despite potential changes in forest composition. We calculated the proportion of forest using a 1000 m moving window via function focal in the package raster for the program R (Core Team 2024; Hijmans 2024). We used the proportion of forest within 1000 m for each map cell as the condition value, with higher values denoting good habitat condition and values of 0 (i.e., no forest) denoting non-habitat. We used a 1000 m moving window based on estimates of the maximum distance traveled by our focal species (Drag et al. 2011; Dodelin et al. 2017; Drag and Cizek 2018) and to fully include heterogeneous habitats, for example, traditional wood-pastures in Central Romania (Hartel et al. 2014; Plieninger et al. 2015).

We considered as top spatial conservation priorities all grid cells falling in the top 30% of the predicted priority ranks (rank values = 0.7–1), which maximizes saproxylic beetles’ representation at the national level and corresponds to the European Union Biodiversity Strategy to achieve 30% protected lands by 2030 (European Commission 2020).

Natura 2000 gap analysis

We evaluated whether the most suitable habitats, as predicted from binary ensemble species distribution, overlapped with the Natura 2000 protected areas (hereafter PAs) sites in Romania. We did not consider other Special Protection Areas Natura 2000 sites in our analysis (i.e., SPAs), as they are specifically aimed at protecting bird species under the Birds Directive. We extracted the potential numbers of Natura 2000 PAs in which each species has a potential presence based on our ensemble projections. We further calculated the percentage of area coverage inside the Natura 2000 PA from the binary distributions. The evaluation was implemented by quantifying the spatial extent of each species distribution area within the current Natura 2000 network and calculating the percentage of the total projected distribution area for each species from the total Natura 2000 Sites of Community Importance area. The area of terrestrial Natura 2000 Sites of Community Importance is approximately 40,500 km2 or ~ 17% of Romania (~ 18% when converting to 1 km2 raster).


Beetle species distributions

The distribution models were evaluated using TSS (True Skill Statistic) and AUC (Area Under the Curve) metrics: TSS ranges from − 1 to + 1, where a TSS of + 1 indicates perfect discrimination, 0 implies random chance performance, and − 1 indicates complete disagreement between the model’s predictions and actual outcomes, while AUC measures the model’s ability to discriminate between presence and absence (scoring from 0 to 1). The species distribution models performed well according to the TSS and AUC evaluation values for all species under the ensemble models (Table 2). The most accurate model was for Cerambyx cerdo, followed by Osmoderma eremita, Lucanus cervus, Morimus funereus, and Rosalia alpina.

Table 2 Evaluation of TSS and AUC for ensemble models by species for current and future scenarios

The relationships between climatic variables and the prediction values for species showed distinct patterns. (supplementary file S2). All evaluated species showed an increase in the probability of occurrence with higher Mean Temperature of the Driest Quarter (BIO9). Most species showed a decrease in the probability of occurrence with higher precipitation levels during the coldest quarter (BIO19), except for Osmoderma eremita, whose distribution is positively influenced by higher precipitation quantities in the same period. All species exhibited a reasonably stable probability of occurrence for all temperature ranges in relation to the Temperature Annual Range (BIO7). The probability of occurrence, however, significantly declined as the temperature increased, with one exception: Rosalia alpina, whose probability of occurrence slightly increased when the temperatures increased.

Regarding the species-specific climatic variable importance included in SDMs, Cerambyx cerdo and Morimus funereus distributions were best explained by the mean temperature of the driest quarter (BIO9). Lucanus cervus distribution was most influenced by the precipitation seasonality (BIO15) and the precipitation of the warmest quarter (BIO18), while Rosalia alpina was influenced by the temperature annual range (BIO7) and by the precipitation of the warmest quarter (BIO18). In contrast, the distribution model of the species Osmoderma eremita was influenced by the precipitation of the coldest quarter (BIO19) (Fig. 1).

Fig. 1

The relative importance of environmental variables used to model saproxylic beetles’ distribution. Box plot = IQR of relative importance of individual models, vertical line = median importance of individual models, whiskers = outliers; black diamond = median importance of ensemble model. Variables used as environmental predictors: BIO9 – Mean Temperature of Driest Quarter in °C, BIO7 – Temperature Annual Range in °C, BIO3 – Isothermality in percent, BIO19 – Precipitation of Coldest Quarter in mm, BIO18 – Precipitation of Warmest Quarter in mm, BIO15 – Precipitation Seasonality in mm

Under the current climate, our models identified suitable habitats for Cerambyx cerdo in the western, southern, and eastern parts of Romania, where most of the oak forests are present (Fig. 2). For Lucanus cervus, substantially uninterrupted suitable habitats were identified in the western part of the country. Most of the suitable habitats for Morimus funereus were found in the west, with some hotspots in the southeastern part of Romania. A more uniform pattern of suitable habitats for Osmoderma eremita was found in the western half of the country, with emphasis on the southwestern and western Carpathian. For Rosalia alpina, the suitable habitats were restricted to the mountainous area in the Carpathians as well as in the central part of Romania.

Fig. 2

Species distribution models for the baseline in Romania (A – Cerambyx cerdo; B –Lucanus cervus; C – Morimus funereus; D – Osmoderma eremita; E – Rosalia alpina). Occurrence data were collected between 2005 and 2022 from public biodiversity databases, peer-reviewed articles and technical reports, and citizen science data from social media

Presence/absence maps also showed larger suitable areas for the saproxylic beetles than suggests published occurrences. In this mapping approach, Osmoderma eremita had the largest distribution area of all saproxylic beetles, > 62,000 km2, followed by Rosalia alpina with a distribution of ~ 57,000 km2, Lucanus cervus with ~ 45,000 km2, Cerambyx cerdo with ~ 42,000 km2, and Morimus funereus with the smallest area of ~ 37,000 km2.

Changes in species distributions

We found the most significant decline in suitable habitats for Cerambyx cerdo, showing an average area loss of 75% for all projections, with a maximum loss of 100% in Had best-case 2061–2080 scenario, and a minimum loss of 35% in the Had best-case 2021–2040 scenario. Lucanus cervus is expected to experience a 99% decrease in suitable habitat according to the Had worst-case 2061–2080 scenario and an increase of 3.5% in suitable habitat under the Had best-case 2021–2040 scenario. In comparison, Osmoderma eremita is likely to lose a maximum of 84% under the Had worst-case 2061–2080 scenario and a minimum loss of 18% under the Had best-case 2021–2040 scenario. Rosalia alpina would lose a maximum of 88% under the Had worst-case 2061–2080 scenario and a gain of 15% in suitable habitat under the Had best-case 2021–2040 scenario. Furthermore, Morimus funereus is likely to suffer the smallest reduction in range, with an average decrease of 25% in the suitable habitat (Table 3).

Table 3 Range change (%) under future climate conditions (two general circulation models × two emission scenarios × three-time horizons)

The HadGEM3-GC31-LL and MIROC6 climate models differed significantly. The Had best-case 2021–2040 showed habitat increases for some species; however, other species are likely to lose over 50% of their habitats, with some losing up to 100%. MIROC6 also forecasted substantial habitat loss, with potential losses between 50% and 90% for several species.

Spatially, the future distribution pattern of Cerambyx cerdo and Lucanus cervus was characterized by a relatively widespread distribution across Romania. In contrast, Osmoderma eremita and Rosalia alpina showed a more localized presence, predominantly in mountainous regions. Finally, Morimus funereus presented a unique distribution pattern, blending aspects of the other four species but still focusing on the Carpathians and northern regions. The HadGEM3-GC31-LL scenario predicted a more drastic change, characterized by a more dispersed distribution of species and a larger loss of suitable area. MIROC6 indicated a more concentrated distribution of species, primarily within and around the Carpathians Mountains, where their habitats are currently located (Fig. 3).

Fig. 3

Projection distribution patterns for five saproxylic beetles for the time horizon 2041–2060 under two emission scenarios: (SSP1-2.6 – best-case scenario & SSP5-8.5- worst-case scenario) and two Global Circulation Models

Conservation priorities for protected saproxylic beetles

Based on the results from the zonation analysis, when integrating the species distribution models for unconstrained prioritization (i.e., without considering Natura 2000 sites as high priority), we found that the top 30% spatial conservation priorities for saproxylic beetles shifted between the current and future climate change scenarios. In the current unconstrained prioritization, the top priority areas for saproxylic beetles (85% of the species distribution retained) overlap the southern, southeastern, and western Carpathians, the Carpathian Foothills, the central part of Romania (Transylvania), and several hotspots in the eastern and southeastern parts of the country (Fig. 4a).

Fig. 4

Priority rank maps for the current climate, with unconstrained prioritization (A) and constrained prioritization (B). Areas have been graded according to their priority rank, with the highest priorities in red (top 30%, ranks 0.7–1). Natura 2000 protected areas (SCIs – Sites of Community Interest) are represented in gray

When comparing the current unconstrained top priority areas with the top priority areas for the Had best-case 2041–2060 unconstrained scenario, the analysis revealed a shift of high conservation values to the eastern Carpathian Mountains, whereas the Had worst-case 2041–2060 unconstrained scenario showed a relatively aggregated pattern of the top priority areas in the Carpathians and a drastic reduction of the priority areas in the Carpathian Foothills and the southeastern region. For the Miroc best-case 2041–2060 scenarios, the top priority areas were distributed in the Carpathians and the eastern part of the country. For the Miroc worst-case 2041–2060 scenario, the analysis revealed more aggregated areas in the Carpathians. For both the best-case and worst-case scenarios, an average of 90% of species distribution is retained within the top 30% of Romania’s land area. Regarding the species, the top priority areas maintained a higher representation (almost 95%) for Osmoderma eremita under Had worst-case 2041–2060 and Miroc worst-case 2041–2060 (Fig. 5).

Fig. 5

Priority rank maps with unconstrained prioritization (A) and constrained prioritization (B) for Had best-case 2041–2060, Had worst-case 2041–2060, Miroc best-case 2041–2060 and Miroc worst-case 2041–2060 scenarios. Areas have been graded according to their priority rank, with the highest priorities in red (top 30%, ranks 0.7–1). The Natura 2000 protected areas (SCIs – Sites of Community Interest) are represented in gray

In the current constrained prioritization, the Natura 2000 network (Sites of Community Interest) encompasses almost 18% of the highest conservation value, retaining an average species representation of approximately 35% for all the species (supplementary file S3). By retaining Natura 2000 sites as areas with the highest priority ranks, the next 12% of the priority conservation areas identified as optimal locations for network expansion (to 30%) revealed potential areas for the Natura 2000 sites in the western Carpathians, southern and southeastern sub-Carpathians, and central part of the country (Fig. 4b), with an average of 85% of the saproxylic beetles’ species representation retained (supplementary file S3).

When comparing the current constrained top priority areas with the top priority areas for the Had best-case 2041–2060 and the Had worst-case 2041–2060 constrained scenarios, the 12% priority areas for conservation revealed a shift of the priority habitats to the higher elevation areas of the Carpathians (Fig. 5). The 12% top priority areas for the Miroc best-case 2041–2060 scenario revealed slightly more areas extended to the eastern Carpathians, whereas for the Miroc worst-case 2041–2060 scenario, the top priority areas revealed a more aggregated pattern to the eastern Carpathians and a reduction of areas in the southwestern Carpathians (Fig. 5).

The constrained prioritization analysis for current and future scenarios revealed that the species with the best representation is Osmoderma eremita, which, when considering expanding to 30% of the landscape as protected, will retain 90% of species representation (supplementary file S3).

Natura 2000 gap evaluation

The Natura 2000 network covers only a small part of the current saproxylic beetle distribution. Specifically, only 13.5% of the Cerambyx cerdo current distribution is protected by the Natura 2000 network. Similarly, the habitat of Lucanus cervus is only 20% covered, while in the case of Morimus funereus, 25% of its habitat falls within the network. The habitats of both Osmoderma eremita and Rosalia alpina are better represented in the Natura 2000 network, with 40% within the network.

The trend is also present in the number of PAs that currently include these species, compared with the number of possible PAs in which the species occur. For example, Cerambyx cerdo is currently listed in only 52 Natura 2000 Standard Data forms, but the current projection overlaps more than 196 PAs from the current Natura 2000 network. In comparison, Lucanus cervus is currently included in 86 Natura 2000 PAs, whereas the current projection overlaps more than 190 PAs. Similarly, Morimus funereus is included in 40 PAs, with the potential to be in more than 130 sites. Out of all the species, the representation of the species Osmoderma eremita is noticeably insufficient within PAs, while it is included in only 16 sites, and the species has a potential distribution in over 180. Lastly, Rosalia alpina is included in only 42 sites, with a potential presence in 194 sites for the current projection.

Future projections also show a drastic reduction in species distribution coverage within the Natura 2000 network. Cerambyx cerdo shows the biggest loss of suitable habitats, with reductions in both the number of protected areas available for its conservation (supplementary file S4) and the extent of which the current Natura 2000 network can cover the species distribution in the future. In most cases, the existing Natura 2000 network inadequately represents this species, covering only 10% of its distribution. In contrast, Rosalia alpina has the best-represented distribution by the Natura 2000 network in future projections. Most projections indicate that, in all future scenarios, the current protected area network can cover approximately 40% of the species distribution (Fig. 6).

Fig. 6

Proportion of saproxylic beetles’ distributions within the protected area network (Natura 2000 SCIs) under current and future projections


Our study revealed that forecasted climate change projections would induce significant shifts in the range of saproxylic beetles, with up to 80% loss of suitable habitats under most future projections. As a general pattern, the saproxylic beetles will restrict their distribution to the Carpathians Mountains due to favorable climatic and environmental conditions and reduce their distribution in the lowland areas. We used a systematic conservation planning approach to identify the conservation priority areas for saproxylic beetles, and we found that for both unconstrained (i.e., current PAs network not considered) and constrained prioritization (i.e., current PAs network not retained as high priority areas) scenarios, the top 30% spatial conservation priorities will shift across the current and future climate change scenarios, with top priority areas overlapping the Carpathian Mountains with an average of 85% of species distribution retained. Under current conditions, the existing Natura 2000 PAs do not perform well for saproxylic beetles, as only a small percentage of their distribution is covered by the current network. Contrary to Bosso et al. (2018), for certain saproxylic beetles (e.g., Rosalia alpina), we found that for future projections, a large area of suitable habitats will be lost, and Cerambyx cerdo will be the species with the greatest potential area lost, followed by Lucanus cervus.

Future distribution changes

Future projections for the analyzed saproxylic beetles’ distributions showed that a large area of suitable habitats will be lost by 2041–2060, but that the analyzed species responded differently to global change. While most species are predicted to lose 80–90% of their current range by 2080, the most drastic changes occurred under the HadGEM3-GC31-LL model predictions. Under the MIROC6 model, the habitat is reduced for most species, which tend to concentrate in the mountainous regions of the Carpathians (Fig. 3). These responses support previous findings in Europe, indicating that under future climate scenarios, there will be a reduction in the distribution of saproxylic beetles, and the species habitat will be restricted to higher altitudinal areas (Poloni et al. 2022).

We found that Cerambyx cerdo is likely to lose almost 75% of its suitable areas due to the species requirements for old-growth oak forested habitats (Redolfi De Zan et al. 2017; Manu et al. 2017; Sabatelli et al. 2023), which can be found at lower elevations and are an ecosystem most affected by future climate change and fragmentation (Parisi et al. 2018), as well as the decline in the number of old growth trees found in wood-pasture and semi-open habitats (Hartel et al. 2013; Redolfi De Zan et al. 2017; Torres-Vila 2017). For Morimus funereus and Rosalia alpina, the averages for all losses or gains in suitable areas are ~ 26% and 30%, respectively. Osmoderma eremita and Lucanus cervus are likely to lose over 50% of their habitats across all projections, a pattern also highlighted by other studies in Europe (Della Rocca and Milanesi 2020) (Table 3).

While some species may lose up to 90% of their suitable areas in the future, several projections show a slight increase (Had best-case 2021–2040 – Lucanus cervus, Morimus funereus, Rosalia alpina), which corroborates findings by Della Rocca and Milanesi (2020), who found that climate change can have a positive effect on beetle distribution. Most projections show a decrease in the distribution of suitable areas across Romania. The least impacted species under both GCMs and SSPs are Morimus funereus and Rosalia alpina. While these species are losing suitable areas, they show a more stable pattern compared with other species with dramatic losses, such as Lucanus cervus and Cerambyx cerdo. Our results highlight the importance of the Carpathians Mountains as suitable habitats for saproxylic beetles under future climate conditions. Saproxylic beetle distributions are predicted to concentrate in the Carpathian region due to favorable climatic conditions, while disappearing from isolated lowland forested patches because of extreme climatic conditions and land use change (Mikolāš et al. 2023). These findings corroborate other studies in Europe, which highlighted that future emission scenarios show a general reduction in suitable habitats for saproxylic beetles and a shift toward higher altitudes (Della Rocca et al. 2019).

Priority areas for saproxylic beetles

The unconstrained spatial prioritization analysis under current conditions showed that the Carpathians Mountains and the western part of the country had consistently high conservation values for saproxylic beetles. These areas represent a refuge for these species under both current and future climate conditions, as they encompass most of the remaining primary and secondary old-growth forests in Romania (Veen et al. 2010; Knorn et al. 2013) (Figs. 4 and 5). Furthermore, based on the two HadGEM3-GC31-LL future scenarios, climate change is anticipated to lead to a decline in priority areas for saproxylic beetles in the lowland areas, especially in the central and eastern regions. In comparison, for the two MIROC6 future scenarios, climate change is predicted to be less impactful, preserving priority areas in the Carpathians and some priority hotspots in the eastern part of the country and western Transylvania (Fig. 5). Some of the factors leading to the contraction of suitable habitat toward higher elevation areas are increased droughts and higher temperatures, as well as a possible reduction in the lowland forested habitats threatened more by agricultural intensification and development (Seibold et al. 2015). Therefore, our findings that higher altitudinal sites will offer favorable climate conditions for saproxylic beetles corroborate those of Bosso et al. (2018) and Della Rocca et al. (2019). The current unconstrained scenario showed that the top 30% of the landscape included ~ 80% of the species distribution and did not overlap completely with Natura 2000 sites. Similarly, Miu et al. (2020) demonstrated that a high proportion of invertebrate species are not covered by Natura 2000 sites. The highest conservation value was observed outside the protected areas in the western and southwestern parts of the country, as well as in the southern and eastern Carpathian Foothills, with small hotspots in the central and eastern parts of the country. The most vulnerable areas were in southern and eastern Romania, with a large reduction in suitable habitats for species in future scenarios. Although Natura 2000 PAs encompass only 18% of Romania’s territory (Miu et al. 2020), ~ 30% of the saproxylic species distributions fall inside protected areas, with a high number of potentially suitable PAs for saproxylic beetles uncovered by their distribution (supplementary file S3).

For the current constrained scenario, the prioritization analysis retained an average species representation of approximately 40% for all saproxylic species for current and future scenarios, less than half of their distribution in Romania (supplementary file S3). The performance of the Natura 2000 coverage in protecting saproxylic beetles was also questioned by several other studies in Europe, which highlighted that a similar percentage or less of the distribution is protected (Bense and Bussler 2003; Viñolas and Vives 2012; Bosso et al. 2013, 2018; Lachat et al. 2013). Rosalia alpina, one of the most charismatic beetle species in Europe (Campanaro et al. 2017), had only 35% of the distribution retained in the Natura 2000 sites, which is close to values stated in the previous studies (supplementary file S3). The species with the least retained distribution in Natura 2000 (32%) is Lucanus cervus, a species inhabiting mature deciduous forests from lowlands and oak forests with a high volume of deadwood at ground level (Bardiani et al. 2017; Méndez and Thomaes 2021). This situation mirrors other locations in Europe. Thomaes et al. (2008) highlighted that 11% of the Lucanus cervus distribution is covered by the Natura 2000 network in Belgium. Osmoderma eremita was the only species whose coverage in PAs increased slightly from 30% under the current scenario to 45% under the Had worst-case 2041–2060 and Miroc worst-case 2041–2060 scenarios.

Finally, under the constrained scenario, when considering an additional 12% of the landscape for the optimal expansion of Natura 2000, an average of 80% of the saproxylic beetle species representation was retained, with the best representation for Osmoderma eremita (90% of distribution protected). Under the constrained scenario, most of the priority areas will be restricted to the Carpathian region, which harbors continuously forested habitats suitable for saproxylic insect development. Therefore, the Carpathians will represent a refuge for the saproxylic beetles in the future.

Effectiveness of the Romanian Natura 2000 PA network for protecting saproxylic beetles

The current Natura 2000 PAs network does not achieve the EU conservation targets for the five listed saproxylic beetle species in Romania. However, it is anticipated that under some future emission scenarios, the saproxylic beetles will experience increased levels of representation. The decrease in species distributions because of climate changes leads to a change in the distribution of species, pushing the distribution of species to higher elevation areas, such as the Alpine biogeographic region, which includes a high number of Natura 2000 PAs (Popescu et al. 2013). As such, we predict a significant increase in the representation of saproxylic insect species by the Natura 2000 network under Had best-case 2021–2040, Had worst-case 2021–2040, and all Miroc best-case and worst-case scenarios.

The current Natura 2000 network covers only a small part of the saproxylic beetle distributions, with gap analysis revealing insufficient representation. The gap analysis results revealed that the current representation of the saproxylic beetles is insufficient within PAs, with most of the species reported in fewer protected Natura 2000 sites than their current potential distribution (supplementary file S4). For example, Osmoderma eremita is reported from only 16 sites, whereas the species has a potential distribution in over 180 sites. There is a knowledge gap regarding this elusive species in Romania, due to its cryptic life history and restrictive habitat requirements. This species requires hollow-bearing old trees with deadwood-specific microclimates for larval development, which can take 3 years or longer (Chiari et al. 2013; Maurizi et al. 2017; Parisi et al. 2019).

The protected area representation by biogeographic region is uneven. A recent study by Cazzolla Gatti et al. (2023) on strictly protected areas in Europe found that in Romania, the Steppic and the Continental biogeographic regions offer very limited protection to biodiversity and rare species; these findings corroborate studies of Miu et al. (2020) and Popescu et al. (2013), which found an urgent need for expanding the protected areas in these regions, as some of the protected saproxylic beetles are dependent on old trees in open or semi-open landscapes (Stan and Nitzu 2013; Manu et al. 2017, 2019). Our predictions are likely to overestimate distributions because microhabitats were not included in the analysis (e.g., old trees, deadwood). These findings underscore the limitations of relying solely on occurrence records for conservation gap assessments. Our results emphasize the significance of incorporating species distribution modeling and spatial planning techniques to estimate the probability of the presence of saproxylic species and the coverage of protected areas, enabling effective management and planning. Our study also provides valuable avenues for future investigations, including the identification of connectivity corridors used by the species and the prediction of suitable habitats in response to climate change. To best achieve the goals of the EU Biodiversity Strategy 2030 of protecting at least 30% of the EU’s land, we urge the expansion of Natura 2000 sites or establishing new protected areas covering the suitable habitats of protected saproxylic beetles.

Data availability

No datasets were generated or analysed during the current study.



IVM is thankful to the Bureau of Educational and Cultural Affairs, United States Department of State with the cooperation of the Romanian – U.S. Fulbright Commission and the Council for International Exchange of Scholars for the Fulbright Visiting Scholar Award.


MDM, LR, SC were supported by the European Commission, LIFE19 NAT/RO/000023 Conservation of saproxylic beetles in the Carpathians LIFE ROsalia (2020 – 2025). IVM was supported by a grant from the Romanian Ministry of Research, Innovation, and Digitisation, CNCS – UEFISCDI, project number PN-III-P1-1.1-PD2021-0378 – “Spatial planning of alien species management activities in Romania”

Author information

Authors and Affiliations
Center for Environmental Research, University of Bucharest, 1 N. Balcescu, Bucharest, 010041, Romania

Marian D. Mirea, Iulia V. Miu, Viorel D. Popescu & Laurentiu Rozylowicz

Ecology, Evolution and Environmental Biology, Columbia University, New York, NY, 10027, USA

Viorel D. Popescu & Bekka S. Brodie

Vrancea Environmental Protection Agency, 2 Dinicu Golescu, Focsani, Romania

Silviu Chiriac

Doctoral School in Geography Simion Mehedinti – Nature and Sustainable Development, University of Bucharest, Bucharest, Romania

Marian D. Mirea


LR and MDM conceived the study; LR, MDM and IVM designed the methodology; MDM and IVM collected the data; MDM, IVM and VDP analyzed the data; LR, MDM and IVM led the writing of the manuscript; VDP, BSB, and SC contributed to the writing of the manuscript. MDM and IVM contributed equally to this work. All authors contributed to the drafts and approved the final version for publication.

Source Priority conservation areas for protected saproxylic beetles in Romania under current and future climate scenarios

Abstract Climate change threatens species and ecosystems globally, including forest ecosystems that support rich invertebrate diversity. Saproxylic beetles, that depend

Directorul general al Regiei Naționale a Pădurilor – Romsilva, Marius – Dan Sîiulescu, a asistat la o aplicație practică în

Au fost desemnați câștigătorii ediției 2024 a concursului interșcolar, intitulat : „Supereroii din lumea insectelor saproxilice” : premiul I –