Next Article in Journal
Diversity of the Piscicola Species (Hirudinea, Piscicolidae) in the Eastern Palaearctic with a Description of Three New Species and Notes on Their Biogeography
Next Article in Special Issue
The Red Coral Community in the Messina Strait: New Findings from the 1700s Lazzaro Spallanzani Collection
Previous Article in Journal
Macro-moth (Lepidoptera) Diversity of a Newly Shaped Ecological Corridor and the Surrounding Forest Area in the Western Italian Alps
Previous Article in Special Issue
Phylogeographic and Bioclimatic Determinants of the Dorsal Pattern Polymorphism in the Italian Wall Lizard, Podarcis siculus
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Diversity
Volume 15
Issue 1
10.3390/d15010096
diversity-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Big Five: Species Distribution Models from Citizen Science Data as Tool for Preserving the Largest Protected Saproxylic Beetles in Italy

by
Lara Redolfi De Zan
1,
Sarah Rossi de Gasperis
2,
Vincenzo Andriani
3,
Marco Bardiani
4,
Alessandro Campanaro
5,
Silvia Gisondi
5,
Sönke Hardersen
4,
Emanuela Maurizi
6,
Fabio Mosconi
6,
Gianluca Nardi
4,
Livia Zapponi
7,
Pasquale Rombolà
8 and
Federico Romiti
8,*
1
DREAm-Italia-Dimensione Ricerca Ecologia Ambiente Soc. Coop. Agr., Via Garibaldi 3, 52015 Pratovecchio Stia, AR, Italy
2
CNR Consiglio Nazionale Delle Ricerche, Istituto di Geoscienze e Georisorse, Via G. La Pira 4, 50121 Firenze, FI, Italy
3
Arma Dei Carabinieri-Comando Unità Forestali Ambientali e Agroalimentari (CUFA)-Reparto Carabinieri Biodiversità di Verona, Via C. Ederle 16/A, 37126 Verona, VR, Italy
4
Reparto Carabinieri Biodiversità di Verona, Centro Nazionale Carabinieri Biodiversità “Bosco Fontana”, Strada Mantova 29, 46045 Marmirolo, MN, Italy
5
CREA Consiglio Per la Ricerca in Agricoltura e l’Analisi Dell’Economia Agraria–Centro di Ricerca Difesa e Certificazione, Via di Lanciola 12/a, 50125 Firenze, FI, Italy
6
CREA Consiglio Per la Ricerca in Agricoltura e l’Analisi Dell’Economia Agraria–Centro di Ricerca Difesa e Certificazione, Via C.G. Bertero 22, 00156 Rome, RM, Italy
7
CNR Consiglio Nazionale Delle Ricerche, Istituto per la BioEconomia, Via F. Biasi 75, 38098 S. Michele all’Adige, TN, Italy
8
Istituto Zooprofilattico Sperimentale del Lazio e Della Toscana ‘M. Aleandri’, Via Appia Nuova 1411, 00178 Rome, RM, Italy
*
Author to whom correspondence should be addressed.
Diversity 2023, 15(1), 96; https://doi.org/10.3390/d15010096
Submission received: 3 December 2022 / Revised: 22 December 2022 / Accepted: 24 December 2022 / Published: 11 January 2023
(This article belongs to the Special Issue Biodiversity in Italy: Past and Future Perspectives)
Browse Figures
Versions Notes

Abstract

:
Background. Volunteers’ participation in scientific research has increased in recent decades. Citizen science (CS) data have been used in quantitative ecology to analyse species ranges by means of species distribution models. We investigated the Italian distribution of five large saproxylic beetles (big five), to describe their niche space, paramount areas for their conservation, and conservation gaps. Methods. CS data from two projects, climate and environmental variables were used to produce Habitat suitability (HS) maps for each species and averaged HS maps. The big five’s conservation status was assessed interpolating HS maps with the distribution of protected areas, concomitantly identifying conservation gaps. Results. The pre-alpine and Apennines arcs, north-eastern Sicily and eastern Sardinia, were identified as conservation’s hotspots. Ranking HS levels from minimum to optimal, the extent of conservation gaps decreases as environmental suitability for the big five increases. Conclusions. For the first time in Italy, CS data have been used to investigate niche space of the largest protected saproxylic beetles and analyse the distribution of their suitable habitat. The resulting HS raster maps and vector layers, reporting HS value in all Italian protected areas (n° 3771), were provided and discussed, reporting an application example for conservation purposes.

1. Introduction

The participation of non-scientist volunteers in scientific research, referred to as citizen science (CS), has consistently increased in recent decades [1]. Following a period of rapid growth and a process of professionalisation, the non-expert contribution to global research can now be considered fundamental in many disciplines [2,3,4,5]. Particularly, CS has proved to be a powerful tool for facing conservation challenges, from public awareness-raising regarding environmental issues to building scientific knowledge [6]. In recent years, CS data have been increasingly used in quantitative ecology to predict current species range by means of species distribution models (SDM) [5,7,8]. SDM, besides being useful in theoretical research on species ecology and distribution, represents a powerful tool in nature conservation and management [9]. Developing robust SDMs needs a large dataset, but collecting occurrence data over a large spatial scale poses a severe challenge to researchers; CS projects offer an effective alternative to tackle this issue [1,4]. The distribution of a species can change in response to a large number of interrelated natural and anthropogenic processes. Climate and land use changes, urbanisation, deforestation, increasing extent of lands earmarked for agriculture and invasive species introduction are the main factors driving the Anthropocene’s biodiversity loss [10,11,12,13]. Even though the processes contributing to the decline of many organisms are well known, many endangered species still have poorly studied geographic distributions [14]. Unfortunately, collecting occurrence data on invertebrates, particularly on endangered insect species, can be difficult, due to their short adult activity period, hidden life cycle, low detection probability, the great sampling effort needed and, in many cases, the engagement of expert entomologists for species identification. In this regard, citizens’ involvement in entomological studies has recently become a pillar of conservation-oriented research [15,16]. In Italy, several CS projects focusing on insects have been recently carried out: Farfalle in ToUr (2013-present), Life MIPP (2012–2017), InNat (2017–2022), Life ESC360 (2018–2022) and X-Polli:Nation (2021-present), just to cite a few. Some of these projects focused on endangered saproxylic species, mainly dealing with the conservation of saproxylic beetles (order Coleoptera). Saproxylic beetles depend upon dead or dying wood during some part of their life-cycle and represent one of the largest taxa contributing to forest saproxylic biodiversity [17,18]. Playing a remarkable role in the trophic web of the forest ecosystem [19,20], they are considered good indicators of forest health [21]. Indeed, saproxylic beetles have been the subject of CS and research projects, which demonstrated their relevance as flagship and umbrella species, whose conservation leads to the preservation of many other forest-dwelling species [22,23,24,25,26].
We selected data on saproxylic beetle occurrence from two CS projects carried out in Italy. Five saproxylic beetle taxonomic units (TU) were chosen: Morimus Brullé, 1832, Lucanus cervus (Linnaeus, 1758), Cerambyx cerdo Linnaeus, 1758, Rosalia alpina (Linnaeus, 1758) and Osmoderma Serville, 1828. These five TUs (hereafter the big five) are easily recognizable by the citizen due to their relatively large size (from 15 mm of the smallest R. alpina to ~90 mm of the largest L. cervus) [27] and thanks to in-depth sheets for each target TU provided in the project-related apps and websites.
For the TU Morimus (Cerambycidae), we considered two taxa occurring in Italy as subspecies of Morimus asper (Sulzer 1776) (hereafter M. asper asper/funereus), as asserted by many authors [28,29,30], omitting the formerly Morimus asper ganglbaueri Reitter, 1894. M. asper asper/funereus adults measure 15–40 mm in length [31], have an elongated oval body and can be chromatically and morphologically distinguished from similar longhorn beetles of the Italian fauna [32]. The European stag beetle L. cervus (Lucanidae) is one of the best-known coleopteran species and the largest stag beetle species in Europe, characterized by a pronounced sexual dimorphism in adults [33]. The species is widely distributed in northern and central regions, whereas only the congeneric Lucanus tetraodon Thunberg, 1806 occurs in the central and southern portion of the Italian peninsula and Sicily [34,35]. Differences in the morphology of the male mandibles make the species easily identifiable from the closely related L. tetraodon [33,36]. Adults of the great Capricorn beetle C. cerdo (Cerambycidae), can be distinguished from other similar (and co-occurring) congeneric species, such as C. welensii (Küster, 1845), by expert entomologists thanks to their heavy sculptured pronotum, blackish body and reddish elytral apex [31,37,38]. Adults of the Rosalia longicorn R. alpina (Cerambycidae) are easily identifiable thanks to their peculiar colour and the black spots on the elytra [31,39,40]. Due to the difficulties in morphological identification of the European Osmoderma spp. (Scarabaeidae), species assignment remains a task for specialists. Hence, for the TU Osmoderma, we considered the species of this genus occurring in Italy, which share the same micro-habitat, inhabiting mature and hallow broad-leaved trees: O. eremita (Scopoli, 1763), O. cristinae Sparacio 1994 and the putative Italian subspecies O. eremita italicum Sparacio, 2000 [41].
Although distributed across Europe in old-growth forests, populations of the big five have been, in the past centuries threatened by forest practices, such as the removal of dead wood, which dramatically reduced their larval habitat [32,42]. The conservation of these species is enforced by the European Habitats Directive 92/43/EC, which lists all TUs in the Annexes II (together with the subspecies M. asper funereus), and three TUs (C. cerdo, R. alpina and Osmoderma spp.) in Annexes IV, with R. alpina and O. eremita as priority species. Additionally, they are classified as Least Concern (C. cerdo, L. cervus, M. asper asper), Near Threatened (R. alpina), Vulnerable (M. asper funereus, O. eremita) and Endangered (O. cristinae) in the Italian IUCN Red List [43].
In view of their relevance as flagship and umbrella species, for the protection of the local saproxylic communities, and considering their conservation status, SDMs of the big five were developed to provide a useful tool to concentrate conservation efforts [22,43].
The present study aims to: (i) identify suitable areas and explore the ecological niches of the TUs; (ii) analyse TU habitat suitability maps in order to highlight hotspots areas for their conservation; (iii) match the identified suitable areas with the Italian protected areas to develop a useful tool for conservation purpose, providing a case of practical application.

2. Materials and Methods

SDMs for the big five were developed using a comprehensive modelling and simulation framework technique based on the R package ‘sdm’ [44]. The generated models’ ensemble was used to predict species potential distribution and highlight any possible gap in their conservation status by comparison with current extent and distribution of protected areas in Italy. Predicted big five distributions were obtained combining remote sensing data on climate and vegetation indices, with a high resolution (~1 km), with species occurrence data (presence and pseudo-absence). Details regarding the specific attributes of occurrence data, environmental variables and modelling are outlined henceforth.

2.1. Occurrence Data

Occurrence data came from two CS projects: a European project financed by the EU Life programme, Monitoring Insects with Public Participation (Life 11/NAT/IT 000252) and a national project financed by the Ministry of the Ecological Transition (InNat). The data came from non-expert participants that used the websites and mobile apps of the CS projects to send records of the target species, accompanied by photographs to allow the validation by experts. Pictures were uploaded on the project websites or through the mobile apps, identified and validated by expert entomologists. The raw dataset consisted of 3538 reports from 865 people, uploaded from 2014 to 2020 (accession data: January 2021), even though some records (n = 125) date back to the early 2000s and have been uploaded later. A data cleaning procedure was performed to increase the quality of occurrence data and the model’s performance [45,46]. Data preparation consisted in removing occurrence data lying outside Italian borders and in eliminating duplicate records. Furthermore, a grid with the same resolution (~1 km), extent and origin of environmental variables was created to select a single occurrence record in each occupied cell, avoiding multiple records which would lead to resampling the same environmental predictors. This procedure led to the following occurrence data per TU: 696 for M. asper asper/funereus, 894 for L. cervus, 124 for C. cerdo, 243 for R. alpina and 32 for Osmoderma spp. (Figure 1).

2.2. Environmental Variables Selection

Given the fundamental role of climate, soil properties, vegetation indices, forest structure and land use for the saproxylic beetle fauna [42,47,48,49], we compiled a dataset consisting of 47 biotic and abiotic variables (Table 1). The set of predictors for the SDMs consisted in 19 bioclimatic variables, 6 seasonal variables (6 × 4 = 24 in total), among which were vegetation and soil water indices, 3 geomorphological variables (altitude, slope and aspect) and the human modification of terrestrial systems index. Climate data were obtained from WorldClim website (https://www.worldclim.org/data/index.html) (accessed on 22 March 2021), recently updated to the 2.1 version (January 2020), which provides 19 bioclimatic variables at very high spatial resolution (~1 km at the equator). These bioclimatic variables represent historical climate data (1970–2000) on annual trends, seasonality and extreme or limiting environmental factors [50]. Elevation above sea level was obtained from WorldClim; slope and aspect were computed using the terrain function of the R package ‘raster’ [51]. Seasonal wind speed, water vapour pressure and solar radiation were calculated from the available monthly variables on WorldClim, dividing seasons according to meteorological breakdown into four three-month periods (i.e., winter starts in December to end in February) [52]. As soil properties represent fundamental factors for the larval stage of saproxylic beetles [53], the soil water index (SWI) from Copernicus Global Land Service database [54] was used as a predictor variable. Daily SWI data were downloaded for a six years period (2015–2020) and averaged to obtain seasonal SWI using the Geostatistical Analyst module of ESRI ArcGisTM. The SWI quantifies the moisture condition at various depths in the soil, providing a 1 km resolution product covering Europe (Version 1) [55], which has been demonstrated to yield high agreement with ground observations [56]. The normalised difference vegetation index (NDVI) and the enhanced vegetation index (EVI) were chosen as independent variables, being correlated to vegetation biomass and canopy biophysical parameters such as photosynthetic activity [57,58,59]. NDVI and EVI were downloaded through the online Data Pool at the NASA Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov) (accessed on 22 March 2021), provided by the MODIS 6 collection products. A monthly temporal aggregation over a ten years period (2000–2010) at a spatial resolution of 1 km provided the chosen features for the data acquisition. Seasonal averages of SWI, NDVI and EVI were then computed. To account for the human impact on the environment, data on global human modification of terrestrial systems (HMTS) were downloaded from the NASA data centre of Earth Observing System Data and Information System (EOSDIS). HMTS index quantifies human modification of lands, at a spatial resolution of 1 km, modelling 13 anthropogenic stressors within five major categories: human settlement, agriculture, transportation, mining and energy production, and electrical infrastructure [60]. The HMTS index output is a continuous 0–1 metric that reflects the proportion of modified landscape.

2.3. Species Distribution Model

As it is likely that the 19 bioclimatic variables from WorldClim are subject to multicollinearity issue [61], to check the overall correlation a variance inflation factor (VIF) stepwise procedure was performed using the vifstep function in the ‘usdm’ R package [62]. The spatial coordinates of occurrence data of each TU (separately) were used to extract values of the 19 bioclimatic variables, then the VIF analysis was performed. Vifstep, calculates VIF for all variables, excludes the one with highest VIF (threshold = 10), then repeats the procedure until no variable with VIF greater than threshold remains [62]. Only variables that had VIF > 10 were excluded from subsequent analyses [63,64]. Four among the most commonly used modelling methods were used in the SDM framework: generalised linear model (GLM), boosted regression tree (BRT), random forest (RF) and maximum entropy (MaxEnt) [65,66]. Pseudo-absence data were generated for each TU, obtaining more than twice the amount of occurrence data, in order to yield the most reliable distribution models [67]. For this purpose, a random method in geographic space was used, which removes points located in presence cells. Two data-splitting procedures, subsampling and bootstrapping, were used to produce partitions of the data. For the subsampling procedure, 30% of the initial data was used. Partitions were then used to train the models, whilst the remaining data were used for models’ evaluation. Both procedures were applied to each modelling method and repeated three times, leading to a total of 24 model outcomes for each TU, considering the whole SDM framework (4 modelling methods × 2 data-splitting procedures × 3 replications). Model outcomes were selected according to their Area Under Curve value (AUC), effective in summarizing the overall accuracy [68], choosing those with an AUC greater than 0.8, i.e., with an excellent discriminating ability [69]. The chosen model outcomes were then combined using a weighted averaging based on true skill statistic (TSS) and a threshold optimization criterion to maximise sensitivity and specificity [44,70]. The result of the best model outcomes’ combination was used: (i) to build five habitat suitability (HS) maps, (ii) to calculate the relative importance of environmental variables for each TU (based on Pearson correlation test) and (iii) to draw niche-space plots, characterising the big five’s ecological requirements. Finally, to provide useful tools for saproxylic beetle conservation, two HS maps were produced summarising the big five’s HS maps. The first HS map (HSmean), resulting from averaging the HS for each TU, highlighted those areas suitable for big five conservation as a whole. To generate the second HS map (HSmax), the maximum HS value of each 1 km2 cell was used, regardless of the TU. The HSmax map can be useful to pinpoint meaningless areas for big five conservation, together with areas where protection’s actions should target few Tus. The HS range (max-min value) of both maps has been divided into five equal ranks (hereafter HS levels), defined as follows: minimum, low, medium, high and optimal HS. Then, the raster layer unique value tool in QGIS [71] was applied to the projected HSmean and HSmax maps (EPSG: 6875, RDN2008/Italy zone) to calculate the extent (km2) of each HS level. Considering that L. cervus, Osmoderma spp. and R. alpina do not occur in Sardinia [28] and L. cervus does not occur in part of central and in southern Italy [33,72], their HS maps were clipped according to the known distributions, to avoid over/under estimation in calculating the HSmean and HSmax maps.

2.4. Gap Analysis

Boundaries of 3948 Italian reserves and national parks were downloaded, as three shapefile layers, from the world database on protected areas [73]. The protected areas predominantly or entirely marine were excluded, leading to a subset of 3771 protected areas. Portions of these protected areas on the sea were excluded by clipping the vector layers with the mainland boundary of Italy. The zonal statistic tool in QGIS was used to calculate basic statistics (mean, median and maximum value) for each protected area, based on HSmean and HSmax raster layers. Then, the resulting statistics were stored in the attribute table of the shapefiles containing the 3771 protected areas (provided in Supplementary Materials). Finally, HSmean and HSmax raster layers were clipped according to the maximum extent of protected areas, merging the three shapefiles.
The output raster layers were reprojected as above described to calculate HS levels’ extent (km2) within protected areas and, by subtraction, the conservation’s gaps. In order to show how our results can contribute to conservation action, we provided a concrete example at a smaller scale. We chose to use the SDM results to identify a conservation gap within a recognized hotspot area, by overlaying the HSmean map (raster file) with the borders of protected areas (shape file).

3. Results

3.1. Species Distribution Model

According to the VIF results, 7 out of 19 bioclimatic variables were excluded from all SDMs due to collinearity problems (VIF > 10). Among the 12 remaining bioclimatic variables, those showing a VIF < 10 for each TU were used as predictors together with the remaining 28 environmental variables (Table 1). Considering the whole SDMs, the RF modelling method reported the best performances, i.e., highest AUC, TSS values and lowest deviance (AUC = 0.90, TSS = 0.71, deviance = 0.67) (Table 2). Details on the 24 model outcomes, for each TU, are reported in Table S1. The outcomes with an excellent discriminating ability, chosen to be combined, were 17 for M. asper asper/funereus, 24 for L. cervus, 9 for C. cerdo, 24 for R. alpina and 13 for Osmoderma spp. The big five HS maps, resulting from the model ensemble, were provided as geoTIFF raster layers (Supplementary Materials and illustrated in Figure 1). Environmental variables’ importance, assessed by the model ensemble of each TU, are reported in Figure S1 (Supplementary Materials). Two niche space plots for each TU were reported to illustrate the role of the more relevant environmental variables (according to the Pearson correlation results) in predicting the occurrence of suitable habitat (Figure 2 and Figure S1). Morimus asper asper/funereus shows the widest distribution among the big five, with the highest HS values in northeast Italy and along the northern Apennines (Emilia-Romagna and Tuscany regions). Suitable habitat of this TU is distributed along the whole pre-alpine arc, within the northern and central Apennines, in southwestern Italy (Lucan and Calabrian Apennines) and along the mountain ranges of northeast Sicily (Figure 1a). Suitable habitat occurs at high autumn NDVI, low autumn SR and where annual mean temperature (BIO1) is in the range 7–13 °C (Figure 2a,b). The highest L. cervus HS values were predicted in the eastern portion of the pre-alpine arc and in the west/north-west Po plain (Piedmont and western Lombardian Prealps). Predicted HS slowly decreases along the northern Apennines till the southern border of the Umbria-Marche Apennines (central Apennines) (Figure 1b). Low values of spring and autumn SR and summer and winter WS positively correlate with high L. cervus HS (Figure 2c,d). The C. cerdo HS map shows an optimum in the central and western Po plain (Figure 1c). The distribution of suitable habitat for this species is positively correlated with HMTS, summer NDVI and summer mean temperature (BIO10) (Figure 2e,f). Although reporting the highest HS value (89%), R. alpina shows the narrowest extent of suitable habitat, limited to those localities in the altitude range 800–2000 m asl, characterised by moderate winter wind speed (4–9 m/s) and an annual mean temperature from 1 to 13 °C (Figure 2g,h). Hotspots for this species overlap the mountain ranges of northern and central Apennines, northeast Prealps, peaks of the southern Apennines (Samnite, Lucan and Calabrian Apennines) and mountain ranges of northeast Sicily (Figure 1d). The suitable habitat for Osmoderma spp., although distributed all over plain, hilly and mountainous areas of Italy, shows the lowest maximum HS value (73%) reported among the big five. The lowest HS values are predicted in the Po plain, along the Murge plateau (southeast Italy) and in the south, central and western part of Sicily (Figure 1e). Summer and autumn NDVI positively correlate with this TU’s presence (Figure 2i). Minimum temperature of the coldest month (BIO 6) within the range −7 °C and 2 °C positively correlates with Osmoderma spp. occurrence, as well as mean temperature of the wettest quarter (BIO 8) within the range −8 °C and 14 °C (Figure 2j). To better identify the conservation’s hotspots for each TU please refer to the raster layers provided in Supplementary Materials (HS_maps_raster_layers.zip). The HSmean map illustrates the distribution of the hotspots paramount for the big five’s conservation (Figure 3a). The Sulcis subregion and the eastern portion of Sardinia are added to the above-mentioned zones, due to the presence of suitable habitat for C. cerdo and M. asper asper/funereus. The HSmax map highlights that the central and northernmost portion of the Alps, the coastal area of Tuscany (the Maremma), the plain in northern Apulia region, the plains along the Italian shores of the Ionian sea, the south-central and western Sicily and the south-central Sardinia present the lowest HS for the big five (Figure 3b).

3.2. Gap Analysis

The overall result of the zonal statistic tool applied to the shapefiles of protected areas was reported as Supplementary Materials (Protected_areas_ITA_vector_layers.zip) and graphically represented showing the differences among protected areas in mean value of HSmean and HSmax maps (Figure 4a,b). The protected areas reporting high value of suitable habitat for the big five as a whole (HSmean) are located in the north-east, in the southern portion of the northern Apennine, in the central Apennines and in the southern Apennine (Calabria region) (Figure 4a). Considering the zonal statistic result for HSmean raster layer, the protected area reporting the highest mean, median and maximum value is the state nature reserve of Sasso Fratino (Emilia-Romagna region). Conversely, the lowest values were recorded in the regional nature reserve Monte Conca (Sicily region). As regards the HSmax map raster layer, the special area of conservation of Torrente Lerada (Friuli-Venezia Giulia region) reported the highest mean and median values, whereas the National Park of Foreste Casentinesi showed the highest maximum value (Emilia-Romagna and Tuscany regions) (Figure 4b). The amount of conservation gaps decreases at increasing HS level, from ~80% of the minimum to ~57% of the optimal HS level, considering the suitable habitat for the big five as a whole (HSmean) (Table 3). On the other hand, the extent of conservation gaps concerning HSmax does not show any decreasing trend (Table 3). A practical example of how protected area shapefiles can overlay HSmean map to highlight conservation gaps and possibly plan targeted actions is illustrated in Figure 5. A relevant area for the big five conservation (marked by red squares in Figure 5) is located in the northern Apennines (central Italy), between SACs (Special Areas of Conservation-Habitat Directive 92/43/CE). This area, consisting of approximately 100 km2, holds a wide territory characterized by high HSmean values (range 43–57%), deserving of being placed under environmental protection.

4. Discussion

This is the first study modelling the distribution of protected saproxylic beetles in Italy using CS data. Our findings confirmed that citizens’ participation in conservation projects, as long as their contributions are validated by experts, can play a crucial role in studying the distribution and habitat preferences of endangered insect species. Involvement of the citizen already proved to be a useful tool in monitoring saproxylic beetle species [74]. Besides, as previously suggested, CS projects can improve knowledge about saproxylic beetle distribution, accelerating the process of gathering occurrence data [16]. In particular, our results pinpointed that suitable habitat has a greater extent for some TUs (i.e., M. asper asper/funereus and L. cervus), while it might be restricted for some others (i.e., R. alpina and Osmoderma spp.). Indeed, M. asper asper/funereus and L. cervus are among the most widespread protected saproxylic beetle species in Italy [16]. A wide area, covering the Pre-alpine arc and the northern Apennines, reported the highest values of suitable habitat either for M. asper asper/funereus and L. cervus. This result is in line with the preferences of these species for mature and well-structured deciduous forests, especially for lowland and medium-altitude oak woodlands [33,43,75]. Conversely, the habitat of Osmoderma spp. and R. alpina can be considered so narrow to be properly described as a series of suitable woodland “islands” that remain in a “sea” of unsuitable land. According to previous findings [76], the suitable habitat for R. alpina occurs mainly in the mountainous altitudinal zone (Figure 2g and Figure 1d), following almost entirely the distribution of beech forest (Fagus sylvatica) [77]. Indeed, even if R. alpina can colonise plenty of deciduous tree species, from the coastline to about 2000 m asl, it is considered a montane species, associated with beech forests [23,78]. The preference of R. alpina for cold climates is supported by the niche space results which indicate that sites with an annual mean temperature within the range 1–13 °C and a low to medium spring solar radiation are suitable for this species. Contrastingly, other TUs, such as C. cerdo, prefer high summer temperatures and present high HS values in the lowland of the Po valley, characterised by medium-high values of HMTS owing to the intensive agricultural activity. These findings corroborate previous considerations about habitat openness and isolated trees as important parameters for predicting the presence of this species [79]. In Italy C. cerdo generally occurs in semi-open wood stands and viable populations can also be found in tree-lined avenues [37,43]. It has to be remarked that, our predicted HS maps are based on remote sensing data, which did not take into account some fine-scale habitat features, such as presence of hallow or dead trees, stands and logs, crucial for saproxylic beetles’ life-cycle. Nevertheless, the highlighted hotspots might facilitate the process of fine-scale assessment of suitable habitat, targeting field surveys where the HS index reported the highest values. In addition, it might be interesting to conduct surveys in those areas where high HS values have been predicted and no data on big five’s occurrence have been reported. This would help both to validate the proposed model and to add new data on the presence of the big five.
To summarise, as it is evident by comparing the HS maps of each TU (Figure 1), the big five showed some differences in the extent and geographical location of suitable habitat, a direct consequence of dissimilarities in their niche space. Nevertheless, the considerations made regarding HS values for each TU might also be applied to some congeneric species, as in the case of Osmoderma spp. For instance, Cerambyx welensii, which is not protected by Habitat Directive, shares the same ecological niche with C. cerdo.
As a consequence of the differences in niche space between TUs, the area of the highlighted hotspots (HSmean map in Figure 3a) has a reduced extent and is marked by lower HS values if compared to the HS maps of single TUs. However, we believe it is critical to focus conservation actions on the highlighted hotspots. Primarily because it might be a way to make the most of the frequently scarce resources earmarked for conservation; secondly, as an umbrella effect, the conservation status of many other species sharing the big five’s habitat, may also be improved. Furthermore, the provided HS maps of each TU (Supplementary Materials: HS_maps_raster_layers.zip) could be useful for conservation managers to highlight more accurately TU-specific sites, within the boundaries of the protected area under their own jurisdiction, crucial for conservation purposes.
Alongside their contribution to conservation management, the provided HS maps for the big five fill the gaps in their distribution and provide insights about their habitat preferences. Besides, HS maps can be considered a baseline for further research aiming at studying the range contraction of these endangered species (extinction risk) under future scenarios [80,81,82]. Monitoring changes in species distribution may indeed be considered an early warning signal for habitat and species loss.
As regards levels of suitable habitat within the network of protected areas, the results of the zonal statistic tool can be compared to highlight the difference between mean value of HSmean and HSmax within each reserve. This match highlights: those protected areas that represent hot spots for the conservation of the big five as a whole, those reserves holding relevant sites for just some target species and the protected areas less valuable for the big five protection (Figure 4a,b). In particular, some protected areas of the central Po plain and the central Apennines (Umbria-Marche Apennine), while being less pivotal for the big five (low HSmean value), own critical sites for the conservation of certain target species (high HSmax value). For instance, the Po plain contains protected areas that show marked improvements in the HS level–from low HSmean to high HSmax–such as the Special Protection Areas of Valli di Novellara and Valli di Gruppo (Emilia Romagna Region). This finding entails that, even an area characterised by a high level of human landscape modification (logging, farming and agriculture activities), might contain suitable sites for species such as C. cerdo or O. eremita, which benefit from individual old growth trees, albeit placed in an anthropized context. The same reasoning could be applied to some protected areas of central Italy, which have suitable sites for M. asper asper/funereus, L. cervus and R. alpina, turning from medium HSmean to high HSmax level (e.g., Gran Sasso e Monti della Laga National Park) (Figure 4a,b).
Finally, as regards the gap analysis results, conservation measures are focused on the big five’s hot spots (HSmean), since nearly the 50% of optimal HS level falls within the protected areas network. Considering HSmax, the absence of a decreasing trend in conservation gaps might be the result of the habitat preferences of some TUs for anthropized landscape, where protected areas are mainly small and far apart (e.g., the presence of suitable habitat for C. cerdo in the Po plain). In this regard, the creation of ecological corridors, to connect protected areas, could be a priority to prevent the risk of local extinction.
From a research point of view, CS data inherently present some biases, e.g., due to the different detection probability of the species in relation to their abundance and behaviours. Among the five TUs, Osmoderma spp. is probably the less suited for CS monitoring in Italy, due to its low detectability by non-expert entomologists. Furthermore, the difference in the number of visitors among protected areas (attractiveness, accessibility, etc.), should be taken into account to remove spatial autocorrelation among occurrence data [83]. By contrast, it should be mentioned that museum data-frequently used in SDM-could suffer from biases too, e.g., spatial and temporal problems due to old records, which require careful examination to avoid introducing error in distribution models [84,85]. On the other hand, we would like to remark that public participation represents a pivotal way to educate new generations and increase their awareness on nature conservation issues [86,87,88]. Indeed, even citizens, though to a lesser extent with respect to researchers, might notice changes in abundance and distribution of target species through a continuative participation in CS projects. Finally, the provided practical example in the northern Apennines allowed us to highlight the process to identify conservation gaps and to target managing actions. Hence, we demonstrated how CS data, once validated and analysed by researchers, could contribute to increase knowledge on species ecology and highlight those areas where conservation actions should be focused. Consequently, the final task of national and local authorities, together with park managers, should be to concentrate their efforts on expanding/establishing protected areas located in the highlighted hotspots.

Supplementary Materials

The following supporting information can be downloaded at: https://doi.org/10.5281/zenodo.7523948: HSmean, HSmax and HS maps of each TU (raster layers); shape files of protected areas with their mean, median and maximum values of HSmean and HSmax raster layers (vector layers); Table S1 and Figure S1.

Author Contributions

Conceptualization, L.R.D.Z., S.R.d.G. and F.R.; methodology, L.R.D.Z., S.R.d.G. and F.R.; software, F.R.; validation, F.R.; formal analysis, L.R.D.Z., S.R.d.G. and F.R.; investigation and resources, L.R.D.Z., S.R.d.G., F.R. and P.R.; data curation, M.B., A.C., S.G., S.H., E.M., F.M., G.N. and L.Z.; writing—original draft preparation, L.R.D.Z., S.R.d.G. and F.R.; writing—review and editing, L.R.D.Z., S.R.d.G., V.A., M.B., A.C., S.G., S.H., E.M., F.M., G.N., L.Z., P.R. and F.R.; visualization, F.R.; supervision, L.R.D.Z., S.R.d.G. and F.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Authors would like to thank the staff involved in the European Project MIPP (Life11/NAT/IT/000252) and in the Italian citizen science project InNat. Authors are particularly grateful to over 800 citizens without whose contribution this work would not have been possible. The authors would also like to thank three anonymous reviewers whose comments and suggestions contributed to the improvement of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Feldman, M.J.; Imbeau, L.; Marchand, P.; Mazerolle, M.J.; Darveau, M.; Fenton, N.J. Trends and gaps in the use of citizen science derived data as input for species distribution models: A quantitative review. PLoS ONE 2021, 16, e0234587. [Google Scholar] [CrossRef] [PubMed]
  2. Cooper, C.B.; Shirk, J.; Zuckerberg, B. The invisible prevalence of citizen science in global research: Migratory birds and climate change. PLoS ONE 2014, 9, e106508. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Danielsen, F.; Jensen, P.M.; Burgess, N.D.; Altamirano, R.; Alviola, P.A.; Andrianandrasana, H.; Brashares, J.S.; Burton, A.C.; Coronado, I.; Corpuz, N.; et al. Multicountry Assessment of Tropical Resource Monitoring by Local Communities. BioScience 2014, 64, 236–251. [Google Scholar] [CrossRef] [Green Version]
  4. Gardiner, M.; Allee, L.; Brown, P.M.; Losey, J.E.; Roy, H.E.; Smyth, R.R. Lessons from lady beetles: Accuracy of monitoring data from US and UK citizen-science programs. Front. Ecol. Environ. 2012, 10, 471–476. [Google Scholar] [CrossRef] [Green Version]
  5. Matutini, F.; Baudry, J.; Pain, G.; Sineau, M.; Pithon, J. How citizen science could improve species distribution models and their independent assessment. Ecol. Evol. 2021, 11, 3028–3039. [Google Scholar] [CrossRef]
  6. McKinley, D.C.; Miller-Rushing, A.J.; Ballard, H.L.; Bonney, R.; Brown, H.; Cook-Patton, S.C.; Soukup, M.A. Citizen science can improve conservation science, natural resource management, and environmental protection. Biol. Conserv. 2017, 208, 15–28. [Google Scholar] [CrossRef] [Green Version]
  7. Moreira, F.S.; Regos, A.; Gonçalves, J.F.; Rodrigues, T.M.; Verde, A.; Pérez, J.A.; Meunier, B.; Lepetit, J.-P.; Honrado, J.P.; Gonçalves, D. Combining Citizen Science Data and Satellite Descriptors of Ecosystem Functioning to Monitor the Abundance of a Migratory Bird during the Non-Breeding Season. Remote Sens. 2022, 14, 463. [Google Scholar] [CrossRef]
  8. Steen, V.A.; Elphick, C.S.; Tingley, M.W. An evaluation of stringent filtering to improve species distribution models from citizen science data. Divers. Distrib. 2019, 25, 1857–1869. [Google Scholar] [CrossRef]
  9. Hernandez, P.A.; Graham, C.H.; Master, L.L.; Albert, D.L. The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography 2006, 29, 773–785. [Google Scholar] [CrossRef]
  10. Clavero, M.; García-Berthou, E. Invasive species are a leading cause of animal extinctions. Trends Ecol. Evol. 2005, 20, 110. [Google Scholar] [CrossRef] [PubMed]
  11. Dirzo, R.; Young, H.S.; Galetti, M.; Ceballos, G.; Isaac, N.J.; Collen, B. Defaunation in the anthropocene. Science 2014, 345, 401–406. [Google Scholar] [CrossRef] [PubMed]
  12. Pimm, S.L.; Raven, P. Extinction by numbers. Nature 2000, 403, 843–845. [Google Scholar] [CrossRef] [PubMed]
  13. Steffen, W.; Grinevald, J.; Crutzen, P.; McNeill, J. The Anthropocene: Conceptual and historical perspectives. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2011, 369, 842–867. [Google Scholar] [CrossRef] [PubMed]
  14. Jetz, W.; McPherson, J.M.; Guralnick, R.P. Integrating biodiversity distribution knowledge: Toward a global map of life. Trends Ecol. Evol. 2012, 27, 151–159. [Google Scholar] [CrossRef]
  15. Silvertown, J. A new dawn for citizen science. Trends Ecol. Evol. 2009, 24, 467–471. [Google Scholar] [CrossRef] [PubMed]
  16. Zapponi, L.; Cini, A.; Bardiani, M.; Hardersen, S.; Maura, M.; Maurizi, E.; Redolfi De Zan, L.; Audisio, P.; Bologna, M.A.; Carpaneto, G.M.; et al. Citizen science data as an efficient tool for mapping protected saproxylic beetles. Biol. Conserv. 2017, 208, 139–145. [Google Scholar] [CrossRef]
  17. Speight, M.C.D. Saproxylic Invertebrates and Their Conservation; Council of Europe, Series 46 Publications and Documents Division: Strasbourg, France, 1989. [Google Scholar]
  18. Lassauce, A.; Paillet, Y.; Jactel, H.; Bouget, C. Deadwood as a surrogate for forest biodiversity: Meta-analysis of correlations between deadwood volume and species richness of saproxylic organisms. Ecol. Indic. 2011, 11, 1027–1039. [Google Scholar] [CrossRef]
  19. Pechacek, P.; Kristin, A. Comparative diets of adult and young threetoed woodpeckers in a European alpine forest community. J. Wildl. Manag. 2004, 68, 683–693. [Google Scholar] [CrossRef]
  20. Stokland, J.N. The Saproxylic Food Web. Biodiversity in Dead Wood; Stokland, J.N., Siitonen, J., Jonsson, B.G., Eds.; Cambridge University Press: New York, NY, USA, 2012; pp. 29–54. [Google Scholar]
  21. Dollin, P.E.; Majka, C.G.; Duinker, P.N. Saproxylic beetle (Coleoptera) communities and forest management and practices in coniferous stands in southwestern Nova Scotia, Canada. Zookeys 2008, 2, 291–336. [Google Scholar]
  22. Buse, J.; Ranius, T.; Assmann, T. An endangered longhorn beetle associated with old oaks and its possible role as an ecosystem engineer. Conserv. Biol. 2008, 22, 329–337. [Google Scholar] [CrossRef]
  23. Duelli, P.; Wermelinger, B. Rosalia alpina L.: Un Cerambicide raro ed emblematico. Sherwood 2005, 114, 19–23. [Google Scholar]
  24. Lachat, T.; Ecker, K.; Duelli, P.; Wermelinger, B. Population trends of Rosalia alpina (L.) in Switzerland: A lasting turnaround? J. Insect Conserv. 2013, 17, 653–662. [Google Scholar] [CrossRef]
  25. Rink, M.; Sinsch, U. Radio-telemetric monitoring of dispersing stag beetles: Implications for conservation. J. Zool. 2007, 272, 235–243. [Google Scholar] [CrossRef]
  26. Russo, D.; Cistrone, L.; Garonna, A.P. Habitat selection in the highly endangered beetle Rosalia alpina: A multiple spatial scale assessment. J. Insect Conserv. 2011, 15, 685–693. [Google Scholar] [CrossRef]
  27. Audisio, P.; Baviera, C.; Carpaneto, G.M.; Biscaccianti, A.B.; Battistoni, A.; Teofili, C.; Rondinini, C. Lista Rossa IUCN dei Coleotteri saproxilici Italiani; Comitato Italiano IUCN e Ministero dell’Ambiente e della Tutela del Territorio e del Mare: Rome, Italy, 2014. [Google Scholar]
  28. Stoch, F. CK2000. Checklist of the Species of the Italian Fauna. Online Version 2.1. 2003. Available online: https://www.faunaitalia.it/checklist/ (accessed on 15 March 2021).
  29. Danilevsky, M.L. A new species of the genus Morimus Brullé, 1832 (Coleoptera, Cerambycidae) from Central Europe. Humanit. Space Int. Alm. 2015, 4, 215–219. [Google Scholar]
  30. Danilevsky, M.L.; Gradinarov, D.; Sivilov, O. A new subspecies of Morimus verecundus (Faldermann, 1836) from Bulgaria and a new subspecies of Morimus asper (Sulzer, 1776) from Greece (Coleoptera, Cerambycidae). Humanit. Space Int. Alm. 2016, 5, 187–191. [Google Scholar]
  31. Bense, U. Longhorn Beetles: Illustrated key to the Cerambycidae and Vesperidae of Europe; Margraf Verlag: Weikersheim, Germany, 1995. [Google Scholar]
  32. Hardersen, S.; Cuccurullo, A.; Bardiani, M.; Bologna, M.A.; Mason, F.; Maura, M.; Maurizi, E.; Roversi, P.F.; Sabbatini Peverieri, G.; Chiari, S. Monitoring the saproxylic longhorn beetle Morimus asper-investigating wood characteristics, season, time of the day and odour traps. J. Insect Conserv. 2017, 21, 231–242. [Google Scholar] [CrossRef]
  33. Bardiani, M.; Chiari, S.; Maurizi, E.; Tini, M.; Toni, I.; Zauli, A.; Campanaro, A.; Carpaneto, G.M.; Audisio, P. Guidelines for the monitoring of Lucanus cervus. In Guidelines for the Monitoring of the Saproxylic Beetles protected in Europe; Carpaneto, G.M., Audisio, P., Bologna, M.A., Roversi, P.F., Mason, F., Eds.; Nature Conservation: Sofia, Bulgaria, 2017; Volume 20, pp. 37–78. [Google Scholar]
  34. Franciscolo, M.E. Coleoptera Lucanidae. Fauna d’Italia; Calderini Edizioni: Bologna, Italy, 1997; Volume XXXV. [Google Scholar]
  35. Lapiana, F.; Sparacio, I. I Coleotteri Lamellicorni delle Madonie (Sicilia) (Insecta Coleoptera Lucanoidea et Scarabaeoidea). Naturalista siciliano (S. IV) 2006, 30, 227–292. [Google Scholar]
  36. Romiti, F.; Redolfi De Zan, L.; Piras, P.; Carpaneto, G.M. Shape variation of mandible and head in Lucanus cervus (Coleoptera: Lucanidae): A comparison of morphometric approaches. Biol. J. Linn. Soc. 2017, 120, 836–851. [Google Scholar] [CrossRef]
  37. Redolfi De Zan, L.; Bardiani, M.; Antonini, G.; Campanaro, A.; Chiari, S.; Mancini, E.; Maura, M.; Sabatelli, S.; Solano, E.; Zauli, A.; et al. Guidelines for the monitoring of Cerambyx cerdo. In Guidelines for the Monitoring of the Saproxylic Beetles protected in Europe; Carpaneto, G.M., Audisio, P., Bologna, M.A., Roversi, P.F., Mason, F., Eds.; Nature Conservation: Sofia, Bulgaria, 2017; Volume 20, pp. 129–164. [Google Scholar]
  38. Villiers, A. Faune des Coléoptères de France. 1. Cerambycidae; Editions Lechavalier: Paris, France, 1978. [Google Scholar]
  39. Bense, U.; Klausnitzer, B.; Bussler, H.; Schmidl, J. Rosalia alpina (Linnaeus, 1758) Alpenbock. In Das europäische Schutzgebietssystem Natura 2000; Petersen, B., Ellwanger, G., Biewald, G., Hauke, U., Ludwig, G., Pretscher, P., Schröder, E., Ssymank, A., Eds.; Ökologie und Verbreitung von Arten der FFH-Richtlinie in Deutschland. Band 1: Pflanzen und Wirbellose. Schr. R. f. Landschaftspfl. u. Natursch; Bundesamt für Naturschutz: Bonn, Germany, 2003; Volume 69, pp. 426–432. [Google Scholar]
  40. Campanaro, A.; Redolfi De Zan, L.; Hardersen, S.; Antonini, G.; Chiari, S.; Cini, A.; Mancini, E.; Mosconi, F.; Rossi de Gasperis, S.; Solano, E.; et al. Guidelines for the monitoring of Rosalia alpina. In Guidelines for the Monitoring of the Saproxylic Beetles protected in Europe; Carpaneto, G.M., Audisio, P., Bologna, M.A., Roversi, P.F., Mason, F., Eds.; Nature Conservation: Sofia, Bulgaria, 2017; Volume 20, pp. 165–203. [Google Scholar]
  41. Audisio, P.; Brustel, H.; Carpaneto, G.M.; Coletti, G.; Mancini, E.; Trizzino, M.; Antonini, G.; De Biase, A. Data on molecular taxonomy and genetic diversification of the European Hermit beetles, a species complex of endangered insects (Coleoptera: Scarabaeidae, Cetoniinae, Osmoderma). J. Zool. Syst. Evol. Res. 2009, 47, 88–95. [Google Scholar] [CrossRef]
  42. Müller, J.; Bußler, H.; Kneib, T. Saproxylic beetle assemblages related to silvicultural management intensity and stand structures in a beech forest in Southern Germany. J. Insect Conserv. 2008, 12, 107–124. [Google Scholar] [CrossRef]
  43. Carpaneto, G.M.; Baviera, C.; Biscaccianti, A.B.; Brandmayr, P.; Mazzei, A.; Mason, F.; Battistoni, A.; Teofili, C.; Rondinini, C.; Fattorini, S.; et al. A Red List of Italian Saproxylic Beetles: Taxonomic overview, ecological features and conservation issues (Coleoptera). Fragm. Entomol. 2015, 47, 53–126. [Google Scholar] [CrossRef]
  44. Naimi, B.; Araujo, M.B. sdm: A reproducible and extensible R platform for species distribution modelling. Ecography 2016, 39, 368–375. [Google Scholar] [CrossRef] [Green Version]
  45. Graham, C.H.; Elith, J.; Hijmans, R.J.; Guisan, A.; Townsend Peterson, A.; Loiselle, B.A.; NCEAS Predicting Species Distributions Working Group. The influence of spatial errors in species occurrence data used in distribution models. J. Appl. Ecol. 2008, 45, 239–247. [Google Scholar] [CrossRef] [Green Version]
  46. Lobo, J.M. More complex distribution models or more representative data? Biodivers. Inform. 2008, 5, 14–19. [Google Scholar] [CrossRef]
  47. Della Rocca, F.; Bogliani, G.; Breiner, F.T.; Milanesi, P. Identifying hotspots for rare species under climate change scenarios: Improving saproxylic beetle conservation in Italy. Biodivers. Conserv. 2019, 28, 433–449. [Google Scholar] [CrossRef]
  48. Holuša, J.; Fiala, T.; Foit, J. Ambrosia beetles prefer closed canopies: A case study in oak forests in Central Europe. Forests 2021, 12, 1223. [Google Scholar] [CrossRef]
  49. Mazur, A.; Witkowski, R.; Kuźmiński, R.; Jaszczak, R.; Turski, M.; Kwaśna, H.; Łakomy, P.; Szmyt, J.; Adamowicz, K.; Łabędzki, A. The Structure of Saproxylic Beetle Assemblages in View of Coarse Woody Debris Resources in Pine Stands of Western Poland. Forests 2021, 12, 1558. [Google Scholar] [CrossRef]
  50. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1 km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 2017, 37, 4302–4315. [Google Scholar] [CrossRef]
  51. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling. R Package Version 3.4-5. 2020. Available online: https://CRAN.R-project.org/package=raster (accessed on 22 March 2021).
  52. Trenberth, K.E. What are the seasons? Bull. Am. Meteorol. Soc. 1983, 64, 1276–1282. [Google Scholar] [CrossRef]
  53. Janssen, P.; Cateau, E.; Fuhr, M.; Nusillard, B.; Brustel, H.; Bouget, C. Are biodiversity patterns of saproxylic beetles shaped by habitat limitation or dispersal limitation? A case study in unfragmented montane forests. Biodivers. Conserv. 2016, 25, 1167–1185. [Google Scholar] [CrossRef]
  54. European Union. Copernicus Land Monitoring Service. European Environment Agency (EEA). 2021. Available online: https://land.copernicus.eu/ (accessed on 23 March 2021).
  55. Bauer-Marschallinger, B.; Paulik, C.; Hochstöger, S.; Mistelbauer, T.; Modanesi, S.; Ciabatta, L.; Massari, C.; Brocca, L.; Wagner, W. Soil Moisture from Fusion of Scatterometer and SAR: Closing the Scale Gap with Temporal Filtering. Remote Sens. 2018, 10, 1030. [Google Scholar] [CrossRef] [Green Version]
  56. Paulik, C.; Dorigo, W.; Wagner, W.; Kidd, R. Validation of the ASCAT Soil Water Index using in situ data from the International Soil Moisture Network. Int. J. Appl. Earth Obs. Geoinf. 2014, 30, 1–8. [Google Scholar] [CrossRef]
  57. Campbell, J.B. Introduction to Remote Sensing, 2nd ed.; Taylor & Francis: London, UK, 1996. [Google Scholar]
  58. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.-M.; Tucker, C.J.; Stenseth, N.C. Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 2005, 20, 503–510. [Google Scholar] [CrossRef] [PubMed]
  59. Didan, K.; Munoz, A.B.; Solano, R.; Huete, A. MODIS Vegetation Index User’s Guide (MOD13 Series); University of Arizona, Vegetation Index and Phenology Lab: Tucson, AZ, USA, 2015. [Google Scholar]
  60. Kennedy, C.M.; Oakleaf, J.R.; Theobald, D.M.; Baruch-Mordo, S.; Kiesecker, J. Managing the middle: A shift in conservation priorities based on the global human modification gradient. Glob. Change Biol. 2019, 25, 811–826. [Google Scholar] [CrossRef] [PubMed]
  61. Graham, M.H. Confronting Multicollinearity in Ecological Multiple Regression. Ecology 2003, 84, 2809–2815. [Google Scholar] [CrossRef] [Green Version]
  62. Naimi, B.; Hamm, N.A.S.; Groen, T.A.; Skidmore, A.K.; Toxopeus, A.G. Where is positional uncertainty a problem for species distribution modelling? Ecography 2014, 37, 191–203. [Google Scholar] [CrossRef]
  63. Montgomery, D.C.; Peck, E.A. Introduction to Linear Regression Analysis; Wiley: New York, NY, USA, 1992. [Google Scholar]
  64. Craney, T.A.; Surles, J.G. Model-Dependent Variance Inflation Factor Cutoff Values. Qual. Eng. 2002, 14, 391–403. [Google Scholar] [CrossRef]
  65. Elith, J.; Leathwick, J.R.; Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 2008, 77, 802–813. [Google Scholar] [CrossRef]
  66. Elith, J.; Phillips, S.J.; Hastie, T.; Dudik, M.; Chee, Y.E.; Yates, C.J. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 2011, 17, 43–57. [Google Scholar] [CrossRef]
  67. Barbet-Massin, M.; Jiguet, F.; Albert, C.H.; Thuiller, W. Selecting pseudo-absences for species distribution models: How, where and how many? Methods Ecol. Evol. 2012, 3, 327–338. [Google Scholar] [CrossRef]
  68. Fielding, A.H.; Bell, J.F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 1997, 24, 38–49. [Google Scholar] [CrossRef]
  69. Hosmer, D.W., Jr.; Lemeshow, S.; Sturdivant, R.X. Applied Logistic Regression; John Wiley & Sons: Hoboken, NJ, USA, 2013. [Google Scholar]
  70. Allouche, O.; Tsoar, A.; Kadmon, R. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 2006, 43, 1223–1232. [Google Scholar] [CrossRef]
  71. QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. 2019. Available online: http://qgis.osgeo.org (accessed on 22 March 2021).
  72. Bartolozzi, L.; Norbiato, M.; Cianferoni, F. A review of geographical distribution of the stag beetles in Mediterranean countries (Coleoptera: Lucanidae). Fragm. Entomol. 2016, 48, 153–168. [Google Scholar] [CrossRef]
  73. UNEP-WCMC; IUCN. Protected Planet: The World Database on Protected Areas (WDPA) and World Database on Other Effective Area-based Conservation Measures (WD-OECM). 2022. Available online: www.protectedplanet.net (accessed on 22 March 2021).
  74. Thomaes, A.; Barbalat, S.; Bardiani, M.; Bower, L.; Campanaro, A.; Fanega Sleziak, N.; Gonçalo Soutinho, J.; Govaert, S.; Harvey, D.; Hawes, C.; et al. The European stag beetle (Lucanus cervus) monitoring network: International citizen science cooperation reveals regional differences in phenology and temperature response. Insects 2021, 12, 813. [Google Scholar] [CrossRef] [PubMed]
  75. Jurc, M.; Ogris, N.; Pavlin, R.; Borkovic, D. Forest as a habitat of saproxylic beetles on Natura 2000 sites in Slovenia. Rev. D’écologie (Terre Et Vie) 2008, 65, 53–66. [Google Scholar] [CrossRef]
  76. Bosso, L.; Rebelo, H.; Garonna, A.P.; Russo, D. Modelling geographic distribution and detecting conservation gaps in Italy for the threatened beetle Rosalia alpina. J. Nat. Conserv. 2013, 21, 72–80. [Google Scholar] [CrossRef]
  77. Wühlisch, G.V. European Beech (Fagus sylvatica). EUFORGEN Technical Guidelines for Genetic Conservation and Use; Bioversity International: Rome, Italy, 2008; 6 p, ISBN 978-92-9043-787-1. Available online: https://www.euforgen.org/publications (accessed on 23 March 2021).
  78. Di Santo, D.; Biscaccianti, A.B. Coleotteri saproxilici in Direttiva Habitat del Parco Nazionale del Gran Sasso e Monti della Laga. Boll. Della Soc. Entomol. Ital. 2014, 146, 99–110. [Google Scholar] [CrossRef]
  79. Buse, J.; Schröder, T.; Assmann, B. Modelling habitat and spatial distribution of an endangered longhorn beetle—A case study for saproxylic insect conservation. Biol. Conserv. 2007, 137, 372–381. [Google Scholar] [CrossRef]
  80. Keith, D.A.; Akçakaya, H.R.; Thuiller, W.; Midgley, G.F.; Pearson, R.G.; Phillips, S.J.; Regan, H.M.; Araújo, M.B.; Rebelo, T.G. Predicting extinction risks under climate change: Coupling stochastic population models with dynamic bioclimatic habitat models. Biol. Lett. 2008, 4, 560–563. [Google Scholar] [CrossRef]
  81. Pereira, H.M.; Leadley, P.W.; Proença, V.; Alkemade, R.; Scharlemann, J.P.W.; Fernandez-Manjarrés, J.F.; Araújopatricia, M.B.; Balvanera, P.; Biggs, R.; Cheung, W.W.L.; et al. Scenarios for global biodiversity in the 21st century. Science 2010, 330, 1496–1501. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  82. Dawson, T.P.; Jackson, S.T.; House, J.I.; Prentice, I.C.; Mace, G.M. Beyond predictions: Biodiversity conservation in a changing climate. Science 2011, 332, 53–58. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. Veloz, S.D. Spatially autocorrelated sampling falsely inflates measures of accuracy for presence-only niche models. J. Biogeogr. 2009, 36, 2290–2299. [Google Scholar] [CrossRef]
  84. Newbold, T. Applications and limitations of museum data for conservation and ecology, with particular attention to species distribution models. Prog. Phys. Geogr. 2010, 34, 3–22. [Google Scholar] [CrossRef]
  85. Soberon, J.M.; Llorente, J.B.; Onate, L. The use of specimen-label databases for conservation purposes: An example using Mexican Papilionid and Pierid butterflies. Biodivers. Conserv. 2000, 9, 1441–1466. [Google Scholar] [CrossRef]
  86. Battisti, C.; Poeta, G.; Romiti, F.; Picciolo, L. Small environmental actions need of problem-solving approach: Applying project management tools to beach litter clean-ups. Environments 2020, 7, 87. [Google Scholar] [CrossRef]
  87. Bonney, R.; Phillips, T.B.; Ballard, H.L.; Enck, J.W. Can citizen science enhance public understanding of science? Public Underst. Sci. 2016, 25, 2–16. [Google Scholar] [CrossRef]
  88. Locritani, M.; Merlino, S.; Abbate, M. Assessing the citizen science approach as tool to increase awareness on the marine litter problem. Mar. Pollut. Bull. 2019, 140, 320–329. [Google Scholar] [CrossRef] [PubMed]
Diversity 15 00096 g001 550
Figure 1. Maps of the suitable habitat for each TU, where the colour scale represents the habitat suitability (HS) from minimum (red) to optimal HS (blue). Black points represent the distribution of the big five occurrence data used in the SMD framework. Reports from the Life MIPP and InNat citizen science projects were selected excluding duplicate occurrences and doubletons within the same ~1 km2 cell. In (a), the HS map of Morimus asper asper/funereus; in (b), Lucanus cervus; in (c), Cerambyx cerdo; in (d), Rosalia alpina; in (e), Osmoderma spp. The L. cervus HS map was clipped in central Italy according to its known occurrence data and considering that only L. tetraodon is present in southern Italy.
Figure 1. Maps of the suitable habitat for each TU, where the colour scale represents the habitat suitability (HS) from minimum (red) to optimal HS (blue). Black points represent the distribution of the big five occurrence data used in the SMD framework. Reports from the Life MIPP and InNat citizen science projects were selected excluding duplicate occurrences and doubletons within the same ~1 km2 cell. In (a), the HS map of Morimus asper asper/funereus; in (b), Lucanus cervus; in (c), Cerambyx cerdo; in (d), Rosalia alpina; in (e), Osmoderma spp. The L. cervus HS map was clipped in central Italy according to its known occurrence data and considering that only L. tetraodon is present in southern Italy.
Diversity 15 00096 g001
Diversity 15 00096 g002 550
Figure 2. Graphical representation of the niche space for the big five: (a,b) Morimus asper asper/funereus, (c,d) Lucanus cervus, (e,f) Cerambyx cerdo, (g,h) Rosalia alpina and (i,j) Osmoderma spp. The colour scale indicates habitat suitability (HS) according to a two-dimensional environmental space, based on the selected predictors, from low HS (dark red) to high HS (dark blue). Environmental variables are reported in the following units: solar radiation (SR) (kJ/(m2 × day)); wind speed (WS) [m/s]; BIO1, BIO6, BIO8 and BIO10 (°C × 10), altitude (m asl). Human modification of terrestrial systems (HMTS) is a continuous 0–1 metric. The valid range and scale factor for NDVI and EVI are, respectively, −2000–10,000 and 0.0001.
Figure 2. Graphical representation of the niche space for the big five: (a,b) Morimus asper asper/funereus, (c,d) Lucanus cervus, (e,f) Cerambyx cerdo, (g,h) Rosalia alpina and (i,j) Osmoderma spp. The colour scale indicates habitat suitability (HS) according to a two-dimensional environmental space, based on the selected predictors, from low HS (dark red) to high HS (dark blue). Environmental variables are reported in the following units: solar radiation (SR) (kJ/(m2 × day)); wind speed (WS) [m/s]; BIO1, BIO6, BIO8 and BIO10 (°C × 10), altitude (m asl). Human modification of terrestrial systems (HMTS) is a continuous 0–1 metric. The valid range and scale factor for NDVI and EVI are, respectively, −2000–10,000 and 0.0001.
Diversity 15 00096 g002
Diversity 15 00096 g003 550
Figure 3. Maps of the suitable habitat for the big five, where the colour scale represents the habitat suitability (HS) from minimum (red) to optimal HS (green). In (a), the HSmean map indicates the hot spots for the big five conservation as a whole, derived from the raster layer produced by averaging the HS of each TU. In (b), the HSmax map, produced by selecting the maximum value of each cell among the TUs’ HS maps, shows those relevant areas for the conservation of at least one TU (green) and meaningless areas for big five conservation (red). For the HS maps of each TU and the raster layers of HSmean and HSmax please refer to the Supplementary Materials (HS_maps_raster_layers.zip).
Figure 3. Maps of the suitable habitat for the big five, where the colour scale represents the habitat suitability (HS) from minimum (red) to optimal HS (green). In (a), the HSmean map indicates the hot spots for the big five conservation as a whole, derived from the raster layer produced by averaging the HS of each TU. In (b), the HSmax map, produced by selecting the maximum value of each cell among the TUs’ HS maps, shows those relevant areas for the conservation of at least one TU (green) and meaningless areas for big five conservation (red). For the HS maps of each TU and the raster layers of HSmean and HSmax please refer to the Supplementary Materials (HS_maps_raster_layers.zip).
Diversity 15 00096 g003
Diversity 15 00096 g004 550
Figure 4. Maps showing the distribution of the Italian protected areas and their relevance for the conservation of big five’s suitable habitat. In (a), each protected area was coloured according to the average HSmean value, dividing the range (min–max) in five HS levels: red = minimum, orange = low, yellow = medium, light green = high and dark green = optimal HS. In (b), the shape of each protected area was coloured according to the average HSmax value, following the same ranking procedure. Average, together with median and maximum value, of both HS maps by protected areas were stored in Supplementary Materials (Protected_areas_ITA_vector_layers.zip).
Figure 4. Maps showing the distribution of the Italian protected areas and their relevance for the conservation of big five’s suitable habitat. In (a), each protected area was coloured according to the average HSmean value, dividing the range (min–max) in five HS levels: red = minimum, orange = low, yellow = medium, light green = high and dark green = optimal HS. In (b), the shape of each protected area was coloured according to the average HSmax value, following the same ranking procedure. Average, together with median and maximum value, of both HS maps by protected areas were stored in Supplementary Materials (Protected_areas_ITA_vector_layers.zip).
Diversity 15 00096 g004
Diversity 15 00096 g005 550
Figure 5. Map showing the habitat suitability hotspot in northern Apennines of the Tuscany region. The area marked by red squares indicates the conservation gap highlighted by overlaying the national protected area boundaries to the HSmean map of the big five.
Figure 5. Map showing the habitat suitability hotspot in northern Apennines of the Tuscany region. The area marked by red squares indicates the conservation gap highlighted by overlaying the national protected area boundaries to the HSmean map of the big five.
Diversity 15 00096 g005
Table
Table 1. Variables used to model habitat suitability of the big five saproxylic beetle species. The results of the variance inflation factor (VIF) stepwise procedure indicate the BIO variables used in the SDM framework for each TU (VIF < 10).
Table 1. Variables used to model habitat suitability of the big five saproxylic beetle species. The results of the variance inflation factor (VIF) stepwise procedure indicate the BIO variables used in the SDM framework for each TU (VIF < 10).
Source (Time Interval)CodeDescriptionVIF Results
WorldClim (1970–2000)BIO 1Annual mean temperatureMaf, Lc, Ra
BIO 2Mean diurnal range: mean of monthly (max temp.-min temp.)Cc
BIO 3Isothermality (BIO 2/BIO 7) (×100)Maf, Lc, Cc, Ra, Os
BIO 4Temperature seasonality (standard deviation × 100)Maf, Lc, Ra, Os
BIO 5Max temperature of warmest month
BIO 6Min temperature of coldest monthOs
BIO 7Temperature annual range (BIO 5–BIO 6)
BIO 8Mean temperature of wettest quarterMaf, Lc, Cc, Ra, Os
BIO 9Mean temperature of driest quarterMaf, Lc, Cc, Ra, Os
BIO 10Mean temperature of warmest quarterCc
BIO 11Mean temperature of coldest quarter
BIO 12Annual precipitation
BIO 13Precipitation of wettest monthMaf, Lc, Cc, Ra, Os
BIO 14Precipitation of driest month
BIO 15Precipitation seasonality (coefficient of variation)Maf, Lc, Cc, Ra, Os
BIO 16Precipitation of wettest quarter
BIO 17Precipitation of driest quarter
BIO 18Precipitation of warmest quarterMaf
BIO 19Precipitation of coldest quarterMaf, Lc, Cc, Ra, Os
AltitudeDigital elevation model (DEM)
SRseasonalSolar radiation
WSseasonalWind speed
WVPseasonalWater vapour pressure
Computed with R softwareSlopeSlope of the cell from WorldClim altitude
AspectOrientation of the cell from WorldClim altitude
Copernicus Land Service (2015–2020)SWIseasonalSoil water index
MODIS Vegetation index products (2000–2010)NDVIseasonalNormalized difference vegetation index
EVIseasonalEnhanced vegetation index
NASA-EOSDIS (2016)HMTSHuman modification of terrestrial systems index
Maf: Morimus asper asper/funereus; Lc: Lucanus cervus; Cc: Cerambyx cerdo; Ra: Rosalia alpina; Os: Osmoderma spp.
Table
Table 2. Summary of the SDM framework results, reporting the actual number of occurrence data, the generated pseudo-absence for each taxonomic unit (TU), together with the means of each calculated metric per modelling method. Metrics’ means were calculated from the model outcomes (n° 24) for each TU (reported in Supplementary Materials Table S1).
Table 2. Summary of the SDM framework results, reporting the actual number of occurrence data, the generated pseudo-absence for each taxonomic unit (TU), together with the means of each calculated metric per modelling method. Metrics’ means were calculated from the model outcomes (n° 24) for each TU (reported in Supplementary Materials Table S1).
Taxonomic Unit (n° Occurrence/Pseudo-Absence)MethodAUCCORTSSDeviance
Morimus asper/funereus (696/1400)GLM0.790.470.471.04
BRT0.800.500.481.10
RF0.880.660.620.80
MaxEnt0.820.530.521.02
Lucanus cervus (894/1800)GLM0.910.700.690.71
BRT0.910.700.690.88
RF0.950.800.780.54
MaxEnt0.930.720.720.69
Cerambyx cerdo (124/300)GLM0.770.440.461.15
BRT0.780.500.481.06
RF0.840.610.570.86
MaxEnt0.770.440.501.11
Rosalia alpina (243/600)GLM0.960.820.850.48
BRT0.960.800.820.68
RF0.980.860.890.38
MaxEnt0.970.840.870.47
Osmoderma spp. (32/120)GLM0.640.240.3816.5
BRT0.840.570.680.82
RF0.850.580.670.76
MaxEnt0.840.540.670.86
AUC: area under the ROC function; COR: correlation; TSS: true skill statistic; GLM: generalized linear model; BRT: boosted regression tree; RF: random forest; MaxEnt: maximum entropy.
Table
Table 3. Results of the raster layer unique value tool applied to the habitat suitability (HS) maps. The extent of the five HS levels was reported for the Italian mainland and for the protected areas (considering their maximum extent). The conservation gaps were reported as the percentage of non-protected areas for each HS level over the Italian mainland.
Table 3. Results of the raster layer unique value tool applied to the habitat suitability (HS) maps. The extent of the five HS levels was reported for the Italian mainland and for the protected areas (considering their maximum extent). The conservation gaps were reported as the percentage of non-protected areas for each HS level over the Italian mainland.
HS LevelHSmean MapHSmax Map
HS Range (%)Italian Mainland (km2)Protected Areas (km2)Gaps (%)HS Range (%)Italian Mainland (km2)Protected Areas (km2)Gaps (%)
Minimum5.13–16.73107,885.4421,899.1579.705.92–22.5658,542.7911,938.4279.61
Low16.74–28.33108,664.6022,049.6079.7122.57–39.1986,013.5819,559.9277.26
Medium28.34–39.9465,765.2514,611.8277.7839.20–55.8391,029.8718,971.2179.16
High39.95–51.5416,074.495244.2467.3855.84–72.4652,795.3710,128.7480.82
Optimal51.55–63.141531.39665.8656.5272.47–89.0911,450.403181.2272.22
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Redolfi De Zan, L.; Rossi de Gasperis, S.; Andriani, V.; Bardiani, M.; Campanaro, A.; Gisondi, S.; Hardersen, S.; Maurizi, E.; Mosconi, F.; Nardi, G.; et al. The Big Five: Species Distribution Models from Citizen Science Data as Tool for Preserving the Largest Protected Saproxylic Beetles in Italy. Diversity 2023, 15, 96. https://doi.org/10.3390/d15010096

AMA Style

Redolfi De Zan L, Rossi de Gasperis S, Andriani V, Bardiani M, Campanaro A, Gisondi S, Hardersen S, Maurizi E, Mosconi F, Nardi G, et al. The Big Five: Species Distribution Models from Citizen Science Data as Tool for Preserving the Largest Protected Saproxylic Beetles in Italy. Diversity. 2023; 15(1):96. https://doi.org/10.3390/d15010096

Chicago/Turabian Style

Redolfi De Zan, Lara, Sarah Rossi de Gasperis, Vincenzo Andriani, Marco Bardiani, Alessandro Campanaro, Silvia Gisondi, Sönke Hardersen, Emanuela Maurizi, Fabio Mosconi, Gianluca Nardi, and et al. 2023. "The Big Five: Species Distribution Models from Citizen Science Data as Tool for Preserving the Largest Protected Saproxylic Beetles in Italy" Diversity 15, no. 1: 96. https://doi.org/10.3390/d15010096

APA Style

Redolfi De Zan, L., Rossi de Gasperis, S., Andriani, V., Bardiani, M., Campanaro, A., Gisondi, S., Hardersen, S., Maurizi, E., Mosconi, F., Nardi, G., Zapponi, L., Rombolà, P., & Romiti, F. (2023). The Big Five: Species Distribution Models from Citizen Science Data as Tool for Preserving the Largest Protected Saproxylic Beetles in Italy. Diversity, 15(1), 96. https://doi.org/10.3390/d15010096

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop