The study was conducted in an area of about 1,025 km² within the Varese Province in the Lombardy region of northern Italy. Elevation in the study area ranges from a few metres above sea level (a.s.l.) near the Ticino River to about 600 m a.s.l at the foothills of Campo dei Fiori, above the Varese lake. The climate is temperate sub-continental (Ferré et al. 2005), generally humid with a mean annual temperature of about 13°C and an annual precipitation of between 700mm and 1,200mm (Prigioni 1995). This area encompasses twelve Natura 2000 sites amongst which were included nine Sites of Community Importance (SCI) and three Special Protection Areas (SPA) (Fig. 1; Suppl. material 1: Table S.1). Due to its strategic geographical location, this area connects the most extensive biodiversity source of the Po Plain, represented by the Valle del Ticino, to the massive mountains of Campo dei Fiori (Casale et al. 2015) and it is an important ecological corridor between the Alps and the Apennines (Casale and Brambilla 2008), the two major “reservoirs of biodiversity” present in the Italian peninsula and included amongst the worldwide Priority Ecoregions for Global Conservation (Olson and Dinerstein 1998). The role of the study area as an ecological corridor is confirmed by the fact that two very ecologically exigent forest species such as the pine marten (Martes martes) and the black woodpecker (Dryocopus martius), which only lived in mountainous and forest areas until a few years ago, recently began colonising lowland areas of the Po Valley (Casale et al. 2015).

Figure 1.

Study area and stag beetle sampling locations (black dots with white circles). Filled black dots indicate investigated sites where the species was not recorded. Black lines indicate the borders of the study area.

Despite its important ecological role and the high biodiversity value of the Ticino Valley Regional Park, which has been acknowledged as the MAB Biosphere Reserve “Valle del Ticino” (UNESCO 2005), the study area crosses the most urbanised area of the country, including the Milano-Malpensa International Airport. The forests of this area are very fragile, fragmented and mainly concentrated along the riverbanks. Wherever broad roads and clearings have been opened up, the dominant tree species, including the pedunculated oak (Quercus robur), black alder (Alnus glutinosa) and hornbeam (Carpinus betulus), have been rapidly replaced by a scrub of exotic species consisting of a dense undergrowth of black cherry (Prunus serotina), black locust (Robinia pseudoacacia), tree of heaven (Ailanthus altissima) and red oak (Quercus rubra) (Bruno et al. 2008).

The target species, the stag beetle, is considered as threatened in several countries within its geographical range (Harvey et al. 2011), listed in the IUCN Red List of Threatened Species as “near threatened” and in the EU Habitats Directive as a priority species of community interest (Appendix II) (Nieto and Alexander 2010). It is considered a focal species for the conservation of suitable habitats for saproxylic beetles (Thomaes et al. 2008) and as an umbrella species, as it represents a diverse and highly vulnerable fauna associated with a wide range of broad-leaved trees (Percy et al. 2000) and coarse woody material (Buse et al. 2008, Ducasse and Brustel 2008). Therefore, by directing management efforts toward the requirements of this species, the requirements of many cohabitant species that use the same habitat are addressed (Roberge and Angelstam 2004). The dataset consists of locations of stag beetles collected in the context of the TRANS INSUBRIA BIONET (TIB) LIFE project (LIFE10/NAT/IT/241) and the AMBROSIANO Regional project. Both of these projects aimed to identify, improve and protect ecological corridors for endangered species in highly fragmented and human-dominated areas. Species occurrences were investigated from 2012 to 2015 in 22 sites, from 2013 to 2015 in 6 sites and from 2014 to 2015 in 6 more sites for a total of 34 sites within the study area. Each site was surveyed every year, for two months (June and July), once a week. To minimise detection bias due to changes in stag beetle activity throughout the day, each site was surveyed in three different time periods: once a month at sunset and three times per month alternatively in the morning or in the afternoon. Each survey consisted of one transect, 500 m length and 10 m width, walking along roads or paths within wooded areas at a constant speed of 17m/min, for a total duration of 30 minutes in dry and low-windy days (Campanaro et al. 2016). The exact locations (X, Y coordinates) and the gender of all specimens observed flying or walking on the ground within transects were recorded.

A set of 18 predictor variables were derived that were contiguously available for the entire study area (Table 1). Specifically, the percentage of oak forests, native broad-leaved forests, invasive broad-leaved forests, mixed forests, beech forests, chestnut woodlands, black cherry stands, black locust stands and other woods were derived from the map of forest types in Lombardy 2016 (vector data; http://www.geoportale.regione.lombardia.it/). The percentage of shrub-lands, croplands, grasslands and dense and sparse human settlements were derived from the map of destination of agricultural usage and forest soils of Lombardy 2012 (vector data; http://www.geoportale.regione.lombardia.it/). The Shannon diversity index of habitats was also calculated considering all predictors mentioned above. The distance to watercourses and to human settlements was calculated from the map of destination of agricultural usage and forest soils of Lombardy 2012, while the distance to roads was estimated from the Open Street Map (http://www.openstreetmap.org). All of these predictor variables were re-sampled at 100 × 100m scale resolution to match the average home-range size of the target species (Sprecher-Uebersax 2003; Thomaes et al. 2008).

Table 1.Download as CSV 

Variables used in the development of stag beetle (Lucanus cervus) Ensemble of Small Models. Variables with Variance Inflation Factor (VIF) > 3 have to be removed due to multi-collinearity with other variables. Average values ± standard deviations at sampled and presence sites are also shown.

Variable	VIF	Sampled sites	Presence sites
Native broad-leaved forests (%)	2.854	26.35 (±38.21)	28.76 (±40.71)
Mixed woods (broad-leaved and coniferous) (%)	2.507	19.06 (±38.25)	19.62 (±39.54)
Oak forests (%)	1.081	9.87 (±23.74)	14.66 (±29.31)
Beech woods (%)	1.001	6.34 (±26.93)	8.85 (±34.07)
Chestnuts woods (%)	1.629	3.76 (±17.21)	6.11 (±21.75)
Distance to watercourses (m)	1.416	804.24 (±576.81)	859.34 (±626.99)
Shrub-lands (%)	1.121	7.85 (±6.11)	9.27 (±7.01)
Grasslands (%)	1.843	14.12 (±22.02)	10.48 (±19.12)
Croplands (%)	2.155	10.01 (±21.31)	10.67 (±23.13)
Invasive broad-leaved forests (%)	1.002	3.14 (±2.25)	3.09 (±2.22)
Black cherry woods (%)	1.006	3.12 (±2.41)	3.21 (±2.51)
Black locust woods (%)	1.979	9.41 (±22.71)	8.95 (±23.32)
Other woods (%)	1.011	4.56 (±12.36)	5.50 (±15.16)
Shannon diversity index of habitats (unitless)	1.307	1.81 (±0.59)	1.71 (±0.56)
Dense human settlements (%)	2.961	2.18 (±4.88)	1.91 (±6.15)
Sparse human settlements (%)	2.989	4.59 (±13.47)	4.76 (±14.41)
Distance to human settlements (m)	2.582	166.57 (±130.69)	180.35 (±148.06)
Distance to roads (m)	1.701	221.89 (±257.03)	167.67 (±200.11)

To avoid multi-collinearity amongst predictors, the Variance Inflation Factor (VIF) was calculated. Following Zuur et al. (2010), predictor variables with VIF values > 3 (highly related to other predictors) should be removed.

The ESMs approach is based on the development of all the possible bivariate models (only two predictors at a time out of a larger set of predictors), followed by their combination into an ensemble (Lomba et al. 2010, Breiner et al. 2015). ESMs were developed with seven widely used modelling techniques (see Appendix S.1). Similar to Breiner et al. (2015), ESMs were developed using locations where the species was recorded and 10,000 random points were generated to serve as pseudo-absence data (including sampled sites where the species was not recorded).

20-fold split sampling (90% training data and 10% test data) was used to evaluate the bivariate models and the resulting ESMs. Similarly to Breiner et al. (2015), all bivariate models with Somers’ D (also known as Gini coefficient; i.e. 2 × AUC-1) values equal to or less than 0 (i.e. AUC ≤ 0.5 and thus worse than a random model) were not included in the ESM. The area under the receiver operating characteristic curve (AUC) and the continuous Boyce index were used to evaluate model performance. ESM outputs in an ensemble prediction (EP) weighted by Somers’ D values were averaged and thus the resulting map (ranging between 0 and 100) was used in the further analyses. Finally, considering the cells with values above the cut-off point identified by the EP, areas of potential presence for stag beetles were identified. The R package ‘ECOSPAT’ (Di Cola et al. 2017) in R version 3.2.4 (R Core Team 2013) was used to develop, evaluate and project ESMs.

A resistance map was derived as the inverse (1 – probability of occurrence) of the resulting map from ESMs and thus it was combined with circuit theory to explore landscape connectivity (McRae et al. 2008). Circuit theory was chosen over other commonly used connectivity methods as circuit theory models movement ecology via random walk pathways across all available movement possibilities (McRae et al. 2008). The software Circuitscape v.4.0.5 (McRae et al. 2008) was used, implementing circuit theory principles and assuming that the intensity of current flow between two nodes (locations) is proportional to the number of times an individual goes from one node to another moving through the particular cell under consideration (McRae et al. 2008). Thus, current flow can be used as a measure of the probability of movement between the two nodes to provide continuous estimates of landscape connectivity for the species considered (McRae et al. 2008). To estimate current flow, Circuitscape simultaneously considers all possible pathways connecting nodes (locations) accounting for an underlying resistance map to movement, producing a map where current density varies across pixels (analogous to the probability of use by random walkers; Koen et al. 2014).

Since high current is produced near nodes, using species locations as nodes could lead to a biased estimation of landscape connectivity (Koen et al. 2014) and thus, similarly to Pitman et al. 2016, a landscape-scale permeability map was developed independent from the former source or destination locations. Specifically, 10–300 random nodes were placed at intervals of 10 around the perimeter of the study area and then Circuitscape was run to identify the optimum number of random nodes required to generate an unbiased landscape connectivity map (Koen et al. 2014). After each iteration, the Pearson correlation was calculated to compare the resulting current density to those estimated, placing 1,000 random locations within the study area; a current density map was considered independent of node placement when the curve comparing correlation coefficients to the number of node pairs reached an asymptote (Koen et al. 2014; Pitman et al. 2016). Subsequently, the final current map was rescaled between 0 and 1 (indicating low and high current, respectively).

SDMs and landscape connectivity based on resistance surfaces are not without caveats. On one hand, the predictive accuracy of SDMs decreases considerably when rare species are considered in the models, as few species’ occurrences and many predictor variables lead to model overfitting and thus reduced generalisation and applicability of the models (Lomba et al. 2010, Breiner et al. 2015). On the other hand, empirical knowledge about the factors facilitating (or hindering) individuals from moving across landscapes is often limited (Graf et al. 2007) and thus resistance surfaces are often derived from expert opinions. These can be potentially biased by the experience of the expert(s), scarce knowledge about the study area and limited or lacking information and published literature on the ecology of species (Zeller et al. 2012, Stevenson-Holt et al. 2014).

However, ESMs can overcome model-overfitting for rare species and thus provide more accurate predictions compared to standard SDMs (Breiner et al.2015) and it has been shown that ensemble predictions perform better compared to single modelling techniques (Marmion et al. 2009). Moreover, SDMs are increasingly and successfully used to derive resistance surfaces, avoiding subjective expert opinions (Wang et al. 2008, Milanesi et al. 2016a, b) and, recently, attention has shifted from single ecological corridor identification between species locations to independent node-based models (Koen et al. 2014).

Despite their high predictive accuracy, ESMs have not been often used to model rare species distribution and they have never been used to derive resistance surfaces to model landscape connectivity. Thus, in this study, a novel application of ESMs has been provided and it is inferred that they might be valuable tools for estimating unbiased landscape connectivity. Since the presence of rare species and landscape connectivity are amongst the most frequently cited criteria for site selection by conservationists (Tuomisto 2010), by combining these two important aspects of biodiversity, ecologists are strongly encouraged to implement this framework in order to develop successful management policies and conservation actions.

Suitable habitats for the stag beetle in the study area were mostly located in the largest patch of contiguous forest along rivers and lakes within the two natural parks, the Ticino Valley Regional Parkand Campo dei Fiori Park. Less suitable habitats were identified in the south-eastern part of the study area, which includes an intensive agricultural matrix, the Milano Malpensa airport and several cities surrounding Milan.

The relationship between the probability of occurrence predicted by ESMs and the predictor variables considered in this study is consistent with the ecological requirements of the species (Thomaes et al. 2008, Harvey et al. 2011). Specifically, response curves showed that forests and grasslands were the most suitable habitats for the stag beetles, while shrub-lands, croplands and human settlements (both sparse and dense) were largely avoided. On a landscape level, woodlands were the primary habitat type for this saproxylic species, since they are the main sources of high amounts of coarse woody debris, (Ranius and Kindvall 2006, Davies et al. 2008). In fact, in many European countries, the stag beetle is mainly associated with oak woodlands and, secondarily, with other tree species such as beech, chestnut and willow (Harvey et al. 2011). As expected, oak and mixed forests were found with the highest probability of occurrence followed by chestnut forests. Some invasive woodland, such as black locust and lack cherry, also seemed to be suitable for the species. In particular, the black locust forests ensured the maximum probability of occurrence for the stag beetle if they extended to about 70% of the landscape. However, over this threshold the stag beetle occurrence decreased rapidly.

This result is of interest as no previous information on the use of invasive woodlands by the stag beetles has been available.

The idea of a negative impact of invasive tree species on native species and ecosystems is generally supported (Maerz et al. 2005; Powell et al. 2011). However, most of these studies refer to herbaceous and shrub invasive plants, while the effect of the widespread invasive tree species on animal communities is still hotly debated. When considering saproxylic beetles, a general avoidance for some exotic coniferous tree species was recorded for Pinus brutia in Israel (Buseet al. 2010) and for Pseudotsuga menziesii in Germany (Müller et al. 2015). A recent study, however, carried out in the temperate forests of northern Italy (Della Rocca et al. 2016), demonstrated that the black locust dead wood, as well as that with oak, hosts a wide range of saproxylic beetles and thus, acts as a surrogate for native deadwood for saproxylic beetles. Some studies supported the hypothesis that the use of exotic plants by native insects increases with extension and time since their introduction (Brändleet al. 2008; Branco et al. 2015). Thus, considering the wide distribution of the black locust and its presence over a long period in the north Italian forests, it is likely that the stag beetle could use this exotic species similarly to other native species, such as oaks.

It was found that the black locust forest can contribute to stag beetle occurrence if its coverage does not exceed 70% of the landscape. This means that the remaining 30% should be represented by other tree species or forest types, especially oak trees and broad-leaved forests. It is possible that, despite the extensive spread of invasive forests, the occurrence of the stag beetle is assured by suitable natives tree species that remain in a small proportion scattered in the invaded area. This occurred in northern Spain, where the stag beetle was found in a Eucalyptus plantation due to the presence of very old chestnut trees remaining within those plantations (Marco Mendez personal comm).

These results showed that the stag beetle remains within the proximity of urban settlements and is positively affected by the presence of roads although it does not seem to be as anthropophilic as in other European countries such as Belgium or Great Britain (Thomaes et al. 2008, Harvey et al. 2011). As suggested by Thomaes et al. (2008), the thermophilous character of the species (Napier 2003) could drive it to the warmer microclimate of cities. Thus, since the Italian mean summer temperatures are considerably higher than those of central and northern Europe (5^th assessment of the Intergovernmental Panel for Climate Change, IPCC AR5WG1 2013; http://www.ipcc.ch), the species remains distant from urban centres and only uses open human infrastructure, such as roads, to facilitate flight and to warm themselves before flight (Harvey et al. 2011).

Considering landscape connectivity, the stag beetle was directly related to the high coverage of woodlands. However, oak forests did not seem to be as important for the species during movement and dispersion compared to other forest types. This phenomenon is probably due to the difference that naturally exists between the habitats used for dispersion and those used for reproduction, as has already been reported for other species (Blàzquez-Cabrera et al. 2016). Actually, while the role of oaks as primary resources for the reproduction of the stag beetle is widely recognised (Thomaes et al. 2008), it is less critical for its dispersion.

Finally, human settlements (both sparse and dense) were the main barriers for movement and dispersion of the stag beetle, while open habitat, especially grasslands and roads, represented important components for the species movement.

A successful management policy for the protection of threatened animal species in human-dominated landscapes should take into account strategies for ensuring the persistence of good-quality habitats and landscape connectivity (Rudnick et al. 2012). In fact, for those forest-dependent species which barely survive in fragmented landscapes dominated by human settlements and intensive agriculture, it is fundamental to plan precise and effective management actions for ensuring suitable breeding sites and minimising costs and efforts (Laycock et al. 2009). This is the case for stag beetles living in the highly urbanised Po Plain of northern Italy.

According to a recent report on the distribution and conservation status of species and habitats of Community Interest in Italy (Genovesiet al. 2014), for almost all saproxylic species, including the stag beetle, the main pressure factors are represented by forestry and silvicultural practices, such as removal of dead wood and stumps for economic and civic purposes and felling old-growth trees in both forests and urban parks. Due to the development of ecological corridors for the stag beetles in this study area, sites were identified along the corridors to maintain ancient woods, forest remnants, hedgerows and old deciduous trees, represented by both native and exotic species. In between this network of strict conservation areas, tree retention and islands of leftover parts of trees and snags after harvest (Heikkala et al. 2014; Vandekerkhove et al. 2013) can be implemented. In this way, larval habitat, represented by deadwood on the ground or rotting stumps, is provided and the survival of stag beetle populations is assured (Van der Sluijs et al. 2004).

These findings also showed that black locust deadwood can be considered in forest management operations aimed at restoring habitats for the reproduction of the stag beetles and which can serve as a temporary food source for the larvae. However, in light of a conservation strategy for this saproxylic beetle, it is essential to ensure the presence of other forest types, especially mixed broadleaved forests and to use black locust only in combination with oak deadwood.