Conservation biogeography of high‐altitude longhorn beetles under climate change

High‐altitude insects are expected to be strongly affected by climate change because of their limited range. Phytophagous species will be subject to further threats because of their dependence on host plants. We investigated the impact of climate change on the distribution of Italian high‐altitude longhorn beetles (Cerambycidae) using a maximum entropy approach based on bioclimatic variables. We used 510 presence records for 15 species distributed throughout the Italian Alps and Apennines. Then, we combined climate‐based predictions with vegetation data to predict the future changes in the extent of suitable areas. All species but two will move uphill to track suitable climates and will face a range contraction (with an average loss of 44%) under both climatic change scenarios considered. Suitable vegetation covers, on average, only 56% of the estimated current species ranges, which means that the future distribution will be even more limited. Given the importance of Italian mountains as hubs of diversity in the Mediterranean hotspot, these results are particularly alarming. Conservation actions that can mitigate the effects of climate change on high‐altitude cerambycids should be focused on contrasting habitat loss and degradation through land preservation and the adoption of appropriate forest management practices.


INTRODUCTION
Human-induced climate change is one of the most dramatic threats to biodiversity (Brodie et al., 2012;Cowie, 2012;Hodkinson et al., 2011;Wagner et al., 2021), possibly surpassing the role of habitat destruction as first cause of extinctions (Bellard et al., 2012).
Insects are predicted to be significantly affected by climate change, with a 1%-2% decline in abundance per year because of various causes, including climate change as one of the principal stressors (Wagner et al., 2021). In principle, insects might respond to climate change through adaptation. However, their capacity to adapt to climate change will depend on the speed of climate change itself, species' life-history characteristics, genetic architecture of key traits and the speed at which a species can change these key traits in response to climate change (Kellermann & van Heerwaarden, 2019), which makes unlikely that most species will adapt to the change in progress. This will harm especially the species living near their physiological limits, like high-altitude insects, since this circumstance reduces plastic responses (Dahlhoff et al., 2019;Hoffmann & Sgro, 2011;Yadav et al., 2021).
Species might track climate change and avoid extinction by changing their distribution. In fact, the fossil record of Quaternary insects shows no high extinction rates in response to climatic oscillations probably because of their ability to track the geographic shifts of tolerable climates (Coope, 1994). Yet, insect species might be unable to change their ranges as fast as needed to cope with the velocity of the current climate change and, at high altitudes, their ability to track a suitable climate might be further limited by geographical constraints.
Altitudinal range shifts are a common response of montane species to chase climate change (Crozier & Dwyer, 2006;Scalercio et al., 2006;Wilson et al., 2007), but high-altitude species living on mountain peaks cannot further shift their range upward (Dirnböck et al., 2011;Shah et al., 2020), which makes them extremely vulnerable to climate change (Halsch et al., 2021;Haslett, 1997). Even when this would be possible, uphill retreats of species with past larger and more continuous distributions at lower elevations will lead to smaller and fragmented distributions (Hodkinson, 2005). Adaptation and survival of such populations isolated on mountain tops will be strongly reduced by genetic drift caused by genetic isolation and founder-effect or bottleneck episodes during range shifts (Arenas et al., 2011). The impact of such changes on insect biodiversity is expected to be particularly momentous, since mountains, in virtue of their isolation, history and environmental heterogeneity, are home to many endemic species (Fattorini et al., 2019).
The Mediterranean basin is recognised as one of the 36 global hotspots of biodiversity (Myers et al., 2000;Fattorini, 2021), and, at the same time, one of the most significantly altered (Pascual et al., 2011). The Italian peninsula, placed in the centre of the Mediterranean hotspot, has one of the richest faunas in Europe, with endemic species mainly concentrated in the two main mountain systems, the Alps and the Apennines (Dapporto, 2010;Menchetti et al., 2021;Ruffo & Stoch, 2005;Urbani et al., 2017).
Here, we investigate the impact of climate change on Alpine and Apennine high-altitude cerambycids (Coleoptera Cerambycidae).
Cerambycids are one of the largest groups of beetles, including 36,000 known species worldwide (Wang, 2017). The Italian fauna includes some 300 species (Sama & Rapuzzi, 2011), which account for about one-third of the whole European fauna (Sama & Löbl, 2010).
Cerambycids, commonly known as longhorn beetles, are one of the most popular groups of insects among both amateur and professional entomologists (Vitali & Schmitt, 2017), so that their distribution in Italy, both past and present, is particularly well known (Sama, 1988(Sama, , 2005. Since cerambycids are also frequently spotted by nature photographers and amateur naturalists (Malmusi et al., 2017;Sama & Rapuzzi, 2011), the recent diffusion of online forums or citizenscience platforms has further refined our knowledge of their distribution in Italy. This large amount of detailed distributional data makes the Italian cerambycids ideal organisms to predict the impact of future climate change on species distribution through habitat suitability models (HSMs) (Elith & Leathwick, 2009;Guisan et al., 2017;Zuur et al., 2009).
For phytophagous insects, climate-based HSM can identify places climatically suitable for a given species, but real occurrence is limited by the presence of the host plant(s). At the same time, the presence of the host plant does not imply the presence of the insect species, because insects and plants may have non-overlapping ranges due to ecological or historical reasons. Therefore, both climatic suitability and host plant presence should be considered to obtain more reliable models of current and future suitable areas.
Since all cerambycids are phytophagous, with larvae developing in dead or living plants, with different degrees of trophic specialisation (Sama, 2002;Şvacha & Danilevsky, 1987), their distribution is also conditioned by that of the host plants, and insect-plant relationships can be used to refine climatic HSM predictions.
In this study, we use HSMs techniques to assess the future conservation status for Italian high-altitude cerambycids under climate change by predicting the extent of future suitable areas. In particular, we test whether cerambycids will show one or more of the following responses: (1) altitudinal range shifts, (2) local extinction, and (3) reduction of the suitable areas.

Data collection and study area
Based on the elevational ranges reported for the Italian cerambycids (Sama, 2002(Sama, , 2005, we selected 25 cerambycid species distributed mainly above 1000 m a.s.l. We chose this elevational threshold as roughly corresponding to the lower limit of the alpine belt environment in the study area (see Fattorini, 2013Fattorini, , 2014Fattorini et al., 2020;Körner & Ohsawa, 2005;Marta et al., 2013). The study area covered the Italian peninsula and Sicily island (Biondi et al., 2013). Sardinia and Corsica were not considered because none of the selected species is present in these islands.
Occurrence data were obtained from literature, collection specimens, entomological forums and citizen-science platforms. Overall, 282 point records were obtained from the literature search. The main source was Sama (2005), from which 198 records were obtained, supplemented by further 84 records from Dioli et al. (1995), Grottolo and Pedersoli (2015), Hellrigl (2010) and Malmusi et al. (2017). Additional unpublished records were obtained from Pierpaolo Rapuzzi's personal database (42 records) and other private collections (94 records, see Acknowledgements). We also scrutinised data present in the two largest Italian entomological forums: Forum Entomologi Italiani (www. entomologiitaliani.net: last accessed 28th April 2019: 28 point records) and Forum Natura Mediterraneo (www.naturamediterraneo.com/ forum: last accessed 28th April 2019, 7 records). Finally, we retrieved data from Global Biodiversity Information Facility (GBIF, 2019a(GBIF, , 2019b(GBIF, , 2019c(GBIF, , 2019dwww.gbif.org;last accessed 19th June 2019;73 records) and in the citizen-science platform iNaturalist (www.inaturalist.org; last accessed 12th May 2019; the data from iNaturalist were all already present in GBIF database, but each observation has been validated on the website of iNaturalist). All data were critically filtered: ambiguous or incomplete data (e.g. uncertain identification, lack of locality names, etc.) were excluded. We considered only data after the year 1960, for consistence with WorldClim data, which refer to the period 1960-1990 (www.worldclim.org, v1.4).
After data collection, checking and filtering, overall, 526 georeferenced point records belonging to 16 species were used for HSM analyses ( Figure S1 and Data S1), which in turn provided reliable results for the following 15 species (510 records (Carpaneto et al., 2015) as 'near threatened' in all cases but for A. pratensis and O. cursor, which are evaluated as 'least concern'. When not directly provided by collector(s), geographical coordinates were obtained from Google Earth Pro 7.3.2.5776 (64-bit) as those of the place names reported on the specimen labels or in the original publications. Coordinates were recorded in WGS84 format.
For each species, a distribution map was built using QGIS 2.18.4 (QGIS Development Team, 2019).

Model building
Current and future species suitable areas were modelled using MaxEnt software (ver. 3.4.1), which is based on a maximum entropy approach coupled with machine learning techniques (Phillips, 2017;Phillips et al., 2006). When properly calibrated, MaxEnt outperforms other presence-only HSM algorithms in predictive accuracy (Merow et al., 2013), especially with small, sparse and irregularly sampled data (Urbani et al., 2015). This is particularly important for studies dealing with rare and endangered species, which are supposed to be present in few localities. The key feature of the programme is to process presence-only data, which avoids absence data that are usually difficult to obtain and can produce misleading results, since they can indicate either a true absence or a low detectability (Elith & Leathwick, 2009;Iannella et al., 2019).
To calibrate the model for current climatic conditions, we obtained an initial set of 19 climatic variables from the WorldClim -Global Climate database (Hijmans et al., 2005; http://www.worldclim. org, ver. 1.4) at 30 arc-sec resolution ($1 Â 1 km UTM cell, temperature values are expressed in C Â 10 and the precipitation in mm): T A B L E 1 Values of area under the curve (AUC) and percentages of contribution of the three most important variables (VAR 1, VAR 2 and VAR 3) resulting from MaxEnt models fitted for high-altitude cerambycids in Italy • BIO1 = annual mean temperature To assess models' discrimination performance, the receiver operating characteristic (ROC) values were considered by means of the corresponding area under the curve (AUC, calculated as sensitivity vs. [1 À specificity]). As an additional indication of model fitting, we also considered the omission rate curve.
MaxEnt parameter settings used in our analyses were as follows: convergence threshold = 0.00001, replicates = 5, replicate run type = cross-validate, regularisation multiplier = 1, maximum number of iterations = 500; other parameters were retained with their default values (see Merow et al., 2013;Morales et al., 2017). We chose the jackknife test option to see the contribution of each variable; the background was created using 10,000 random points. A tenth-percentile training threshold (Freeman & Moisen, 2008) was used to binarize predictions (Elith et al., 2011;Lahoz-Monfort et al., 2014). This threshold omits all areas with a habitat suitability lower than the suitability values for the lowest 10% of records, thus assuming that the 10% of records in the least suitable habitat are not occurring in regions that are representative of the species overall habitat. This threshold is more cautious than a threshold based on the lowest predicted suitability, which may be strongly affected by outliers represented by inconsistences in georeference or identification of species (e.g. Escalante et al., 2013). scenario is the less optimistic one (Sanford et al., 2014). Among the models available in Worldclim, we chose the CNRM-CM5 climate model (Voldoire et al., 2013), since it is particularly appropriate for the study area (Urbani et al., 2017).

Range changes
The binary maps generated by processing MaxEnt outputs were han-  (1), loosing suitability (2), remaining suitable (3), gaining suitability (4), occupied with future climate (5), occupied with current climate and with suitable land cover (6), occupied with future climate and suitable land cover (Corine 1) (7), occupied within suitable land cover (Corine 2) (8). Central insets with the arrow indicate the vertical shifts (in metres) of species rangess a suitability status change index (SSCI) was calculated for each cell as the intersection between present and future distributions (Ceccarelli & Rabinovich, 2015). SSCI index can assume the following values: À1 (loss of suitability), 1 (stable suitability) and 2 (gain of suitability).

Relationships between species distributions and habitats
Since vegetation plays a critical role in determining cerambycid spe-  (3), gaining suitability (4), occupied with future climate (5), occupied with current climate and with suitable land cover (6), occupied with future climate and suitable land cover (Corine 1) (7), occupied within suitable land cover (Corine 2) (8). Central insets with the arrow indicate the vertical shifts (in meters) of species ranges IUCN, 2016). To maintain a conservative approach, these analyses were applied only to the 4.5 scenario.

Elevation shifts in species ranges
Since uphill shift in elevational range is one of the most commonly observed patterns of species response to climate change (Shah et al., 2020;Wilson et al., 2007), maps of current and future distributions were interpolated with shapefiles of elevation, and differences between current and future average, minimum and maximum elevations were tested using paired Student's t-tests (one-tailed, with the null hypothesis being no shift, and the alternative hypothesis being an uphill shift).

Influence of climatic variables and model evaluation
The models showed good discrimination performance for 15 out of 16 species, with AUC > 0.750 (mean AUC AE SE = 0.853 AE 0.016); the only species with a poor fit, Callidium coriaceum, was excluded from further analyses ( Table 1). The most influential variables were BIO19 (precipitation of the coldest quarter, for 10 species), followed by BIO15 (precipitation seasonality, 8 species) and BIO1 (annual mean temperature), BIO9 (mean temperature of driest quarter), and BIO18 (precipitation of warmest quarter) (5 species) (Table 1, Figure S2). All species whose distribution is restricted to the Alps were influenced by precipitation, namely BIO15, with exception of Brachyta interrogationis. The Apennine species were mainly influenced by temperature, in particular BIO2 (mean diurnal range), BIO8 (mean temperature of wettest quarter), BIO7 (temperature annual range), and by winter precipitation (BIO19). The species occurring in both the Alps and the Apennines were influenced by both temperature and precipitation.
High probability values were predicted for BIO8 in the interval À5 C to 5 C, steeply declining to zero at 10 C. All the 10 species influenced by BIO19 (Table 1)   showed a range decline greater than 50% and four a decline greater than 30%, with an average loss of 44% and an average gain of 6% (Table 2). Under the 8.5 scenario, declines were larger; 10 species showed a decline greater than 50% and two species a decline greater than 30%, with an average loss of 62% and, again, an average gain of 6% (Table 2).
With the 4.5 scenario, patterns of distribution changes indicated that most Alpine species will be affected by a range contraction especially along the peripheral mountains (e.g. the Prealps) and in mountains near the sea, whereas loss in the central Alps will be less marked (Figures 1 and 2). For these species, overall loss of suitable area ranged about 30%-73% (see insets in Figures 1 and 2). Judolia sexmculata (Figure 1d) (Figures 3 and 4a), range contractions will particularly affect the Northern Apennine sector, with overall loss ranging between 13% and 52% (see insets in Figures 3 and 4). Three species living both in the Alps and in the Apennines would experience gains of suitable areas in the Alps but virtually no gain in the Apennines: Evodinus clathratus (Figure 3a, with an overall gain of 25%), Oxymirus cursor (Figure 3b, with an overall gain of 4%) and Tetropium gabrieli ( Figure 4a, with an overall gain of 6%). For these species, suitability loss was 33%, 52%, and 40%, respectively. The overlap between climate-based models and suitable vegetation showed a lack of suitable habitats in a large fraction of the species' climatically suitable area (see Table 4 and values of categories 6-8 in insets of Figures 1-4). On average, only 56% of the area T A B L E 4 Suitable areas covered by forests and protected areas for high-altitude cerambycids in Italy and areas covered by protected areas climatically suitable is occupied by suitable vegetation. Even assuming that the fraction of suitable land cover will remain stable, a future negative trend is observed in the extent of predicted distributions for all species (Table 4). On average, only 39% of the future climatically suitable areas will be included in current protected land (Table 4).

DISCUSSION
Several studies using HSMs have shown the importance of climatic factors in constraining species ranges of cerambycids (e.g. Aguilar et al., 2016;Bosso et al., 2018;Kadej et al., 2017;Lachat et al., 2013;Peterson & Scachetti-Pereira, 2004;Rukavina et al., 2018;Silva et al., 2016). We found that the distribution of high-altitude cerambycids on Italian mountains is strongly constrained by climatic variables. In particular, Alpine species are mainly limited by very high values of winter precipitation, which corresponds to a high abundance and long persistence of snow cover (Harris et al., 2019). The only exception is Brachyta interrogationis, that is the only non-forest species, living in open grasslands. A possible explanation is that the influence of snowfall and precipitation is less important for herbaceous vegetation and for the larva of this species, which feeds in rhizomes of Geranium sylvaticum (Sama, 2002). In contrast, Apennine species, which live in a drier and warmer climate, are mainly limited by the temperature in the wettest quarter, a parameter that also influences the presence of snowfall in Mediterranean regions, where winter is the wettest season (Mooney et al., 2001). Therefore, future climate changes affecting snow cover persistence will have important consequences on their distribution.
In response to increasing temperature, species may track more favourable climates by moving latitudinally (polewards) and upwards (towards higher elevations). As the vertical temperature gradient is more rapid than the latitudinal one, species can compensate for warming more easily moving upwards than northwards, and uphill movements in response to climate change are now well documented for a variety of taxa (e.g., Cerrato et al., 2019;Chen et al., 2011;Doak & Morris, 2010;Harsch et al., 2009;Merckx et al., 2013;Moritz et al., 2008). However, elevational shifts cannot be used by all species to track climate change; for example, species that are already on mountain tops have no place to go. Moreover, because of the roughly conical shape of mountains and the isolation of mountain tops, upslope range shifts, even when possible, will transform previously large and continuous ranges into smaller and fragmented distributions (Fattorini et al., 2020).
In the case of strong selection, sufficient genetic variation and heritability of the relevant traits, some species might undergo rapid evolutionary changes generating evolutionary rescue (Catullo et al., 2019). However, even in the hypothesis of high heritability for the traits associated with climate adaptation, given the long generation time (2-3 years for most cerambycid species, Şvacha & Danilevsky, 1987Şvacha & Danilevsky, , 1988 and the scarce connectivity that many species show in the model, such a rescue seems unlikely. In addition, a recent meta-analysis (Diamond, 2017) of experimental values of heritability of heat-tolerance in insects (which is probably a relevant character for climate change adaptation) reported modest values of heritability, with higher values in species living in areas with broader temperature variations.
In the Alps and Apennines, many high-altitude species occur in mountain tops and therefore cannot shift upwards; thus, the risk of local extinctions is high. This expectation is largely confirmed by our findings. All analysed species are predicted to shift towards higher altitudes and hence to reduce their ranges. Although our study was restricted to a single group of beetles, we expect that other groups of high-altitude insects would show the same patterns. As Alpine and Apennine high-altitude faunas are rich of endemics Urbani et al., 2017), the impact of climate change on the Mediterranean biodiversity hotspot appears particularly alarming.
Climatic suitability is a necessary, but not sufficient condition for the presence of phytophagous beetles, because of their dependence on the co-occurrence of their host plant(s). Some insights from leaf beetles (Chrysomelidae) suggest that host plants are a potentially important predictor of species ranges, although less than climate (Cerasoli et al., 2019) and these two factors together (climate suitability and habitat availability) proved to be the most effective to predict the northwards shifts in response to climate change (Platts et al., 2019). Thus, whereas the predicted loss of suitable area is realistic and probably underestimated, the gain is more uncertain, because in the gained areas the host plants might not be present.
To take into account the dependence of cerambycid species on their host plants, we used the extent of suitable vegetation based on CLC maps to refine our climate-based predictions. We found that roughly half (56%) of the climatically suitable area also presents a vegetation which likely includes the host plants. As there is no available projection of the future distribution of CLC, to calculate species future areas of presence we optimistically assumed that this percentage will remain stable. Even if the most common response of vegetation to climate change is shifting upwards to find cooler temperatures, there is no guarantee that this will happen. Firstly, several studies did not find a significant upward shift of forest vegetation over the last decades (Scherrer et al., 2020). Secondly, even when an upslope movement can occur, this might be more limited by tree demography (as host tree may have a long-life cycle that does not allow a rapid colonisation), and competition with other species (more plastic and with shorter life cycles), than by dispersal possibilities (Scherrer et al., 2020). Currently, observed shift rates in most plant species seem to be insufficient to keep up with climate change (Chen et al., 2011;Corlett & Westcott, 2013;Loarie et al., 2009). Actually, on the Alps, there is evidence of ongoing upward shifts in plant elevational ranges, which is interpreted as a result of both climate change and decreasing grazing pressure (Frei et al., 2010;Leonelli et al., 2011;Vitasse et al., 2021;Wieser et al., 2019). However, woody plants show, on average, a positive shift in optimum elevation of about 33 m per decade, which appears too slow to track isotherm shifts induced by climate warming (Vitasse et al., 2021).
Also, there is evidence that some plants could shift downwards instead of upwards in more humid habitats to track water availability which is a limiting factor Lenoir et al., 2010).
This can also explain our findings that for some cerambycids the upslope movement will be less pronounced with the 8.5 than with the 4.5 scenario, probably because of nonlinear relationships with climate or because of an important role of other factors, such as water availability.
Although there is indication that some terrestrial insects might have an upward shift of their leading edge within the range of the pace of climate warming, or higher (Vitasse et al., 2021), responses are species specific and are not necessarily valid for the cerambycids as a whole. Even assuming that both plants and beetles will move upwards with a sufficient velocity, this does not mean that beetles will be really able to use the newly formed forests. The development of appropriate new habitats uphill needs long time periods, especially for species adapted to mature forests with abundance of dead wood, such as xylophagous cerambycids. In fact, the simple presence of the host plants is not a sufficient condition for most of the investigated species, which require old forests with abundance of dead wood (Sama, 2002;Şvacha & Danilevsky, 1987Şvacha & Danilevsky, , 1988, a resource that is clearly lacking in a recently grown forest. This highlights the importance of having pristine or properly managed forests for the conservation of cerambycids. Small gains of climatic suitability can be also virtually null, if the climatically gained areas do not have the host plants or if the habitats are not properly managed. This becomes even more evident considering that, on average, only 39% of the predicted suitable species ranges will be included in protected areas, which could guarantee the protection of habitats and a careful management of forest resources. It is important to stress that, in addition to climate change, many other forms of anthropogenic pressures affect mountain areas (Fattorini et al., 2020), with possible detrimental effects on cerambycids and their host plants. For example, like other xylophagous insects, cerambycids are threatened by deforestation, invasion of exotic plants, wildfire, inappropriate forest management, pollution, and so on (Cálix et al., 2018). Thus, in addition to range contraction and population fragmentation due to climate change, high-altitude cerambycids might suffer from various forms of habitat loss and degradation. For these reasons, it is important a more widespread adoption of appropriate forest management practices. For example, traditional silviculture in Italy considers negatively the presence of dead trunks and wood (Carpaneto et al., 2015). As a result, in many managed forests, the undergrowth is systematically 'cleaned' from dead wood, even in protected areas, which has serious negative effects on cerambycids and other xylophagous beetles.
This study was restricted to the Italian territory, although many species have much wider distributions. This means that they might be more variable in habitat requirements than estimated from our data.
However, species with disjunct or fragmented distributions tend to develop local adaptations when the dispersal capabilities of the species do not overcome the distance (Savolainen et al., 2013;Storz, 2005). Therefore, including data from the entire species ranges might result in the opposite problem of overestimating niche sizes, and hence species plasticity of local populations. In fact, data on butterflies indicate that several Alpine and Apennine populations represent distinct lineages and that the Alps and the Apennines are to be considered different functional refugia during climatic cycles . On this basis, it might be suggested that not differentiating Alpine and Apennine populations might lead to an overestimation of true species plasticity. We have decided to collectively use all Italian data as a compromise between the need of avoiding the risk of mixing populations with potentially too different adaptations (entire ranges) and that of excessively underestimating niche size (separate analyses for the Alpine and the Apennine populations). Finally, we have overlooked the possibility that species might change habitat requirements under global change. For example, species might be less demanding in substrate conditions (such as the diameter of dead wood, its state of degradation, or the range of suitable tree species). However, these adaptations would require important changes in morphological, physiological and developmental characteristics that are unlikely to evolve in short time.

CONCLUSIONS
Our results indicate that, even under an optimistic scenario, climate change will produce strong contractions in high-altitude cerambycids in Italian mountains. In response to increasing temperature, mountain species may track more favourable climates by moving upwards.
However, while tracking favourable climatic conditions by uphill movements, mountain cerambycids will experience a substantial reduction of their climatically suitable areas. Moreover, since not all the climatically suitable area of a given species is occupied by its host plants, the true suitable area will be even more reduced. Given the importance of Italian mountains as hubs of diversity in the Mediterranean hotspot, these results are particularly alarming for the conservation of montane biotas. Conservation actions able to mitigate the effects of climate change on high-altitude cerambycids should be focused on contrasting habitat loss and degradation through land preservation and the adoption of appropriate forest management practices.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website.
Data S1 Distributional records used in the study.