Conservation Biogeography of Tenebrionid Beetles: Insights from Italian Reserves

: The species-area relationship (SAR), the latitudinal gradient, the peninsula e ﬀ ect, and the elevational gradient are widespread biogeographical patterns. Using data from Italian reserves, these patterns were tested for tenebrionids and used as a framework to calculate expected extinction rates following area loss. Area was an important determinant of overall tenebrionid species richness, but not for xylophilous and endemic species. Thus, focusing on reserve areas is not the best approach for conserving insects with specialised ecology and restricted distribution. In general, species richness declined northwards, which contrasts with the peninsula e ﬀ ect, but conforms to the European latitudinal pattern observed in most taxa because of current and past biogeographical factors. Minimum elevation had an overall negative inﬂuence, as most tenebrionids are thermophilic. However, xylophilous tenebrionids, which are mainly associated with mesophilic forests, did not decline northwards, and were positively inﬂuenced by higher elevational ranges that allow more forms of vegetation. SAR-based extinction rates reﬂect species dispersal capabilities, being highest for geophilous species (which are mainly ﬂightless), and lower for the xylophilous species. Extinction rates based on multiple models indicate that the use of area alone may overestimate extinction rates, when other factors exert an important role in determining species richness.


Introduction
Italy is located in the centre of the Mediterranean basin, one of the world's biodiversity hotspots [1][2][3]. Italy shows an extraordinary diversity of species for the most disparate plant and animal groups [4,5], as a consequence of both current factors, such as the variety of climatic conditions and landscapes [6], and the past history of its territory, in particular its role as a crucial refuge and evolutionary centre during Pleistocene glaciations [7][8][9][10].
Italy hosts some 900 protected areas [11], and another 400 areas that benefit from some form of protection [6], accounting for a total surface area of about 3.5 million hectares (more than 11% of Italian land). For most areas, however, faunal knowledge is still very poor, especially for invertebrates, which, of course, represents a serious limit to effective conservation actions.
Among the few invertebrates for which there is information about the number of species present in some reserves, beetles play a prominent role, with relatively abundant data on their distribution, at least for the best investigated families [4].
Tenebrionid beetles (Coleoptera Tenebrionidae) are a highly diversified group of mainly saprophagous beetles, although there are many species that feed on fungi and lichens and others that are predators or semi-predators [12]. The Italian fauna comprises about 400 taxa, i.e., species or subspecies, as the current taxonomic dividing line between species and subspecies is arguably arbitrary in most cases [13]. high altitudes; as many tenebrionids are thermophilic species [12,35], a negative relationship can be expected. This negative effect should be particularly apparent if maximum and/or minimum elevation are used as predictors of species richness. It is known that tenebrionid richness tends to decline with elevation, as the cold climatic conditions of high altitudes are not favourable for them [31,33]. Thus, we expect that, if most of the species are Mediterranean elements mainly associated with lowland areas, species richness should decrease rapidly with increasing minimum elevation. If most of the fauna is composed of more tolerant species vertically distributed from lowland to middle elevations, only very high altitudes would impact profoundly on species richness, and hence a negative correlation with maximum elevation should be detected. Finally, the average elevation of each study area can be a rough indication of its overall orographic physiognomy, since a high average elevation, especially if associated with a small elevational range, means that most of the area is at high elevation, and species richness is expected to be negatively affected by this measure.
To explore how diversity patterns can be influenced by species' ecology, I conducted analyses not only for the family as a whole, but also for two main ecological groups separately: (1) species that occur in the soil (geophilous species), and (2) those associated with wood (xylophilous species). As endemics are considered of higher conservation priority [36], based on their geographical distribution, I also divided the species into widespread (occurring also out of Italy) and endemic (occurring only in Italy).

Materials and Methods
Tenebrionid species richness was determined by reviewing all available published papers dealing with the Italian tenebrionids. These references were checked to extract occurrence records from Italian preserved areas. This led to the selection of 18 reserves or reserve assemblages, since contiguous reserves were considered to a be a single reserve [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54] (Figure 1, Table 1). All selected reserves were the object of intense sampling, in most cases through several years or even decades of insect collection, by expert taxonomists using a variety of sampling techniques. In most cases, species lists were explicitly indicated as virtually complete by the author(s) of the study. All occurrence data present in the sources were checked in light of present knowledge, and the nomenclature was standardised to a current taxonomy. I considered both species and subspecies, but the word "species" will be used for simplicity. Table 1. Number of species for various beetle groups in Italian reserves. Location, area, and elevation of reserves are also given. Res: Reserve number (reserves are numbered as in Figure 1), Lat: Latitude (decimal degrees), Long: Longitude (decimal degrees), Area: reserve size (km 2 ), Min: Minimum elevation (m), Max: Maximum elevation (m), Range: Elevational range (m), Mean: Mean elevation (m). Tenebrionid groups: S: Total species richness, G = Geophilous species, X = Xylophilous species, W = widespread species, E = Endemic species. I omitted from the analysis the genus Lagria and the subfamily Alleculinae (formerly considered as a separate family), because of their peculiar ecology (in contrast to other tenebrionids, they are flower-visiting insects) and paucity of distributional data [13,33]. Each species was classified according to the main lifestyle of the adults either as geophilous (i.e., adults that occur in the soil) or xylophilous (e.g., adults that occur in rotten wood or on living trees) using information reported in Aliquò et al. [55] and personal observations. Endemic status was assessed using data reported in Aliquò et al. [55] and Löbel and Smetana [56]. I considered as endemic the species with distribution restricted to Italy and possibly with extremely small extensions to adjacent areas (e.g., small, non-Italian islands close to Italy). The complete dataset assembled for the present study is reported in Supplementary Materials Table S1.

Res
When the study areas perfectly matched the borders of the reserves, I used the reserve size as the size of the study area. In some circumstances, study areas did not correspond exactly to the borders of a reserve, as they were comprising only a part of a reserve (reserve 16), including adjacent contiguous areas (reserves 6 and 18) or arising from the aggregation of continuous reserves (reserves 9 and 10). In such cases, I used as size of the study area the surface explicitly reported in the references used to extract the faunal data, or obtainable from the maps provided in these references.
When authors of faunal studies reported the elevation of their sampling sites, I used these values to calculate minimum elevation, maximum elevation, elevational range (maximum minus minimum), and mean elevation for the respective study areas (reserve 9 [57], reserve 13 [58], reserve 16 [59], reserve 17 [60]). For reserve 4, I referred to the minimum, maximum, and mean elevation reported by Hardersen [51]. For reserve 18, I calculated the average elevation using minimum and maximum values as given by the authors of the faunal study [37], integrated with topographic maps. In other cases, I calculated these values from topographical maps.
Selected reserves were distributed from the Alpine Region to Sicily, therefore being distributed along virtually all the Italian latitudinal gradient. Reserve size spanned three orders of magnitude, from 1 to 1700 km 2 , which ensures a sufficiently large variation for the power function model (if the threshold in area size at which the slope of the SAR begins to decline is not met by the largest of the sampling areas, the relationship might be best modelled by a linear relationship instead of a power function) [61]. Reserves also spanned three orders of magnitude in elevational range, as well as in minimum, maximum, and mean elevations, including both strictly coastal and mountainous reserves ( Figure 1, Table 1).
The species-area relationship is typically best modelled by the Arrhenius power function: where S is the number of species, A is the area size, and c and z are fitted parameters [62,63]. An alternative model is the Gleason semilogarithmic function: The power function is frequently applied in its linearized (log-log) form [62,63]: In this form, c is the expected number of species per area unit, and z is the slope of the function [64]. To investigate the SAR, I tested both the Arrhenius and the Gleason models. The Arrhenius model was fitted using its linearized form with values of species richness and area log 10 -transformed. In all cases, the Arrhenius function showed values of the Akaike Information Criterion (AICc) largely inferior to those of the Gleason function (with ∆AICc as follows: 144.1 for all species, 117.1 for geophilous species, 113.3 for xylophilous species, 140.6 for widespread species, and 80.4 for endemic species). Therefore, only results from the Arrhenius power function model were considered. Values of c and z of the Arrhenius function for the various groups were compared using analysis of covariance (ANCOVA) [64]. Because of the presence of zero values, number of geophilous species, number of xylophilous species, and number of endemics were log 10 (x+1)-transformed.
To evaluate the relative importance of the other environmental variables in determining species richness, I used a multi-model selection procedure based on the corrected Akaike Information Criterion (AICc) [65]. In multi-model selection, all variables were log 10 -transformed before analysis (minimum elevation was log 10 (x+1)-transformed because of the presence of areas with a minimum elevation of 0 m), and then tested individually and in all their possible combinations. Use of logarithms assumes that elevation effects combine with area in a multiplicative way [66]. Models were ordered by decreasing AICc values, and the model with the lowest AICc was selected as the best model, but alternative models with ∆AICc values ≤ 2 were also considered [65].
According to the SAR, if the original area A 0 is reduced to A 1 , the original number of species S 0 is expected to decline to S 1 [22,67,68], following the equation: This approach allows the calculation of expected local extinction rates in isolated blocks of fragmented habitats (such as protected areas), following area loss [67,68]. I used this approach to calculate the number of species lost following area reduction from the A 0 down to zero.
Although the standard way to calculate extinction rates is the application of Equation (4) with the z-value obtained from the power function of the SAR (Equations (1) and (3)), variables other than area may exert a relevant role in determining species richness. Thus, I also calculated extinction rates using, as z-values in Equation (4), the coefficients for area of the best-fit multiple models, which take into account the influence of elevation and latitude. This approach allows the calculation of expected species loss by area reduction (holding stable the influence of all other variables) with an exponent z that takes into account the influence of other variables on species richness. This assumes that the suitable surface of a certain study area can be reduced by habitat fragmentation, loss, or alteration, but its geographical position and orography remain unchanged.
Analyses were performed in R version 3.5 [69]. The library MuMIn [70] was used for the multi-model selection. Collinearity was assessed with the Variance Inflation Factor (VIF) using the package mctest [71]. VIF values > 5 were considered as indicative of collinearity.

Results
SARs for the whole dataset and the various groups (Figures 2 and 3, Table 2) showed similar z-values (ANCOVA: f = 0.348, p = 0.845). c-values decreased in the order: total species (9.5 species per area unit) > widespread species (8.0 species per area unit) > xylophilous species (5.3 species per area unit) > geophilous species (4.1 species per area unit) > endemic species (2.3 species per area unit). On average, the average number of xylophilous species per reserve (mean = 14.167) was similar to that of the geophilous ones (mean = 12.778) (Student's t-test for paired samples: t = 0.721, p = 0.481). By contrast, the number of endemics (mean = 4.222) was always lower than that of the widespread species (mean = 22.722) (t = −5.647, p << 0.001).
Pairwise correlations (Pearson coefficient) indicated a significant correlation of species richness with area in all cases, while latitude was always correlated negatively with richness (albeit not significantly for the xylophilous species) (Figure 4). Elevational range was always correlated positively with richness (albeit not significantly for the endemics).
For the whole dataset, the best-fit models (Table 3) included: (1) area, latitude, and minimum elevation, and (2) latitude, minimum elevation, and elevational range. Area exerted a positive, albeit weak, influence. Latitude and minimum elevation had negative influences, whereas range had a positive influence. This approach allows the calculation of expected local extinction rates in isolated blocks of fragmented habitats (such as protected areas), following area loss [67,68]. I used this approach to calculate the number of species lost following area reduction from the A0 down to zero.
Although the standard way to calculate extinction rates is the application of Equation (4) with the z-value obtained from the power function of the SAR (Equations (1) and (3)), variables other than area may exert a relevant role in determining species richness. Thus, I also calculated extinction rates using, as z-values in Equation (4), the coefficients for area of the best-fit multiple models, which take into account the influence of elevation and latitude. This approach allows the calculation of expected species loss by area reduction (holding stable the influence of all other variables) with an exponent z that takes into account the influence of other variables on species richness. This assumes that the suitable surface of a certain study area can be reduced by habitat fragmentation, loss, or alteration, but its geographical position and orography remain unchanged.
Analyses were performed in R version 3.5 [69]. The library MuMIn [70] was used for the multi-model selection. Collinearity was assessed with the Variance Inflation Factor (VIF) using the package mctest [71]. VIF values > 5 were considered as indicative of collinearity.

Results
SARs for the whole dataset and the various groups (Figures 2 and 3, Table 2 Table 2. Parameters of species-area relationships (SARs) for tenebrionids in Italian reserves. SARs were modelled using the linearized version of the power function: log(S) = log(c) + z log (A), where S is the species number, A is area, and c and z are fitting parameters. SE: standard error, t = Student's t, p = probability, R 2 = goodness of fit. Pairwise correlations (Pearson coefficient) indicated a significant correlation of species richness with area in all cases, while latitude was always correlated negatively with richness (albeit not significantly for the xylophilous species) (Figure 4). Elevational range was always correlated positively with richness (albeit not significantly for the endemics).  Table 2. Parameters of species-area relationships (SARs) for tenebrionids in Italian reserves. SARs were modelled using the linearized version of the power function: log(S) = log(c) + z log (A), where S is the species number, A is area, and c and z are fitting parameters. SE: standard error, t = Student's t, p = probability, R 2 = goodness of fit.  For the whole dataset, the best-fit models (Table 3) included: (1) area, latitude, and minimum elevation, and (2) latitude, minimum elevation, and elevational range. Area exerted a positive, albeit weak, influence. Latitude and minimum elevation had negative influences, whereas range had a positive influence.
For the geophilous species (Table 3), the best-fit models included: (1) latitude, maximum elevation, and mean elevation, (2) area, latitude, and minimum elevation, and (3) area, latitude, maximum elevation, and mean elevation. Area exerted a positive influence, and latitude and minimum elevation negative effects, as in the best-fit models with all species. In the geophilic species, however, there was a positive effect of maximum elevation and a negative effect of mean elevation. For the xylophilous species (Table 3), the best-fit model included only elevational range, with a positive influence. VIF values indicated collinearity for maximum elevation and mean elevation in the first and third models. These were the only two models affected by collinearity.
For the widespread species (Table 3), the best-fit model included area, latitude, and minimum elevation. Area exerted a positive, albeit weak, influence, whereas both latitude and minimum elevation had negative effects. Finally, for the endemic species (Table 3), the best-fit models included: (1) latitude (negatively), minimum elevation (negatively), and mean elevation (positively), and (2) latitude (negatively), minimum elevation (negatively), and maximum elevation (positively). Table 3. Best-fit models for the influence of area, latitude, and elevation on tenebrionid richness in Italian reserves. DF: degrees of freedom, R 2 adj: adjusted R 2 , AICc: corrected Akaike Information Criterion. Errors refer to standard errors. * = p < 0.05, ** = p < 0.01, *** p < 0.001.   Table 3. Best-fit models for the influence of area, latitude, and elevation on tenebrionid richness in Italian reserves. DF: degrees of freedom, R 2 adj : adjusted R 2 , AICc: corrected Akaike Information Criterion. Errors refer to standard errors. * = p < 0.05, ** = p < 0.01, *** p < 0.001.  For the geophilous species (Table 3), the best-fit models included: (1) latitude, maximum elevation, and mean elevation, (2) area, latitude, and minimum elevation, and (3) area, latitude, maximum elevation, and mean elevation. Area exerted a positive influence, and latitude and minimum elevation negative effects, as in the best-fit models with all species. In the geophilic species, however, there was a positive effect of maximum elevation and a negative effect of mean elevation. For the xylophilous species (Table 3), the best-fit model included only elevational range, with a positive influence. VIF values indicated collinearity for maximum elevation and mean elevation in the first and third models. These were the only two models affected by collinearity.

Intercept
For the widespread species (Table 3), the best-fit model included area, latitude, and minimum elevation. Area exerted a positive, albeit weak, influence, whereas both latitude and minimum elevation had negative effects. Finally, for the endemic species (Table 3), the best-fit models included: (1) latitude (negatively), minimum elevation (negatively), and mean elevation (positively), and (2) latitude (negatively), minimum elevation (negatively), and maximum elevation (positively).
Use of the SAR to predict expected species loss following area reduction of 25% indicates an overall species loss (local extinction) of about 6% ( Figure 5, solid line). The most sensitive group is the geophilous species, with an expected extinction of more than 7% of species with a loss of 25% of area (Figure 6a, solid line). The xylophilous component would lose about 5% of species (Figure 6b, solid line), and the widespread one about 5% (Figure 6c, solid line), with 25% of area lost. Use of the SAR to predict expected species loss following area reduction of 25% indicates an overall species loss (local extinction) of about 6% ( Figure 5, solid line). The most sensitive group is the geophilous species, with an expected extinction of more than 7% of species with a loss of 25% of area (Figure 6a, solid line). The xylophilous component would lose about 5% of species (Figure 6b, solid line), and the widespread one about 5% (Figure 6c, solid line), with 25% of area lost.  The endemic species would lose about 4% of species (Figure 6d, solid line), but it is important to note that the SAR was poorly fitted in this case. For an area loss of 50%, the expected extinction rates would be: about 11% for the total tenebrionid diversity, 17% for the geophilous component, 11% for the xylophilous component, 14% for widespread species, and 10% for endemics. The coefficients for area obtained from multiple models were slightly lower than those from the power function, thus reducing the expected extinction rates (Figures 5 and 6, dotted lines).  The endemic species would lose about 4% of species (Figure 6d, solid line), but it is important to note that the SAR was poorly fitted in this case. For an area loss of 50%, the expected extinction rates would be: about 11% for the total tenebrionid diversity, 17% for the geophilous component, 11% for the xylophilous component, 14% for widespread species, and 10% for endemics. The coefficients for area obtained from multiple models were slightly lower than those from the power function, thus reducing the expected extinction rates (Figures 5 and 6, dotted lines).

Discussion
The z-values obtained for the tenebrionids of Italian reserves were between z = 0.15 and z = 0.26, thus being consistent with the range found for isolates (0.2-0.4), and very close to those found for habitat islands (z = 0.22, on average) [62,63]. The highest c-value was found for the widespread species, and the lowest for the endemic species; namely, the number of species per area unit among the widespread species is about four times that of the endemic species. This is a reflection of the fact that widespread species are always more numerous than endemic species. The number of xylophilous species per area unit was slightly superior to that of geophilous ones. This can be explained by the lower dispersal ability of geophilous species (which belong mostly to the subfamily

Discussion
The z-values obtained for the tenebrionids of Italian reserves were between z = 0.15 and z = 0.26, thus being consistent with the range found for isolates (0.2-0.4), and very close to those found for habitat islands (z = 0.22, on average) [62,63]. The highest c-value was found for the widespread species, and the lowest for the endemic species; namely, the number of species per area unit among the widespread species is about four times that of the endemic species. This is a reflection of the fact that widespread species are always more numerous than endemic species. The number of xylophilous species per area unit was slightly superior to that of geophilous ones. This can be explained by the lower dispersal ability of geophilous species (which belong mostly to the subfamily Pimeliinae and are typically flightless) compared with xylophilous species (which include many flying species) [13,72].
Contrary to expectation, area was not always a strong determinant of tenebrionid diversity in Italian reserves. It was included in models dealing with the total species, the geophilous species, and the widespread species, but not in those for the xylophilous and the endemic species. These findings suggest that area may be a relatively good predictor for generalist and widespread species, whereas more specialised and range-restricted species are more influenced by other factors. In terms of conservation planning, these findings suggest that large reserves tend to have more species but might not automatically preserve a good number of the most specialised and range-restricted, which contrasts with the idea that few large reserves would be better than several small [16].
Latitude had a negative effect in all cases except for the xylophilous species, where it did not enter the best-fit model. The impact of latitude is particularly evident in reserve 1, which shows a very low richness and is much below the SAR regression line. This is the northernmost area and was used in the past as a willow plantation. While it is important as a wet area, its vegetation has been profoundly altered by man, which may have contributed to decrease its diversity. The negative influence of latitude is in accordance with the overall latitudinal gradient documented for most taxa in the Northern hemisphere [24][25][26][27][28], including tenebrionids [8][9][10], and with the latitudinal pattern already documented for the Italian tenebrionids [73]. This latitudinal gradient is opposed to what can be expected according to the peninsula effect evoked to explain a southward decrease in species richness along the Italian peninsula observed in other taxa [29]. The inverse relationship between latitude and richness can be explained by two non-mutually exclusive explanations: (1) most tenebrionids are thermophilic insects, and hence they are filtered northwards by decreasing temperatures, and (2) the southern parts of the Italian peninsula and Sicily acted as refugial and possibly speciation centres during Pleistocene glacials, where most of the fauna survived [73]. Interestingly, latitude was not an important predictor for xylophilous species. These species include many tenebrionids associated with mesophilic forests and are therefore distributed through the Italian peninsula, occurring in the southern parts, especially on mountain areas, where temperatures are similar to those found at lower elevations in northern areas.
Minimum elevation had a negative influence in all cases except for the xylophilous species. This is consistent with the thermophilic preferences shown by most tenebrionids, and their decline in species richness with elevation [31,33]. Mean elevation also had a negative influence on geophilous species. High-altitude areas are typically considered important centres of speciation [32]. In the case of Italian tenebrionids, the negative correlation between minimum elevation and number of endemics suggests that, for these mostly thermophilic beetles, the filtering effect exerted by lower temperatures overwhelmed the role of isolation as a factor promoting speciation in mountain areas [33], although there are some cases of endemic taxa restricted to high altitudes on the Apennines [73]. The presence of high-altitude endemics explains the positive influence of maximum and mean elevation on endemics, and the influence of maximum elevation on geophilous species, because the few high-altitude endemics are also geophilous species. Elevational range positively influenced total species richness and xylophilous species. A greater elevational range means a greater variety of habitats, and this can be particularly important for xylophilous species, because a higher elevational range allows the presence of more forms of vegetation (vegetational belts) [74].
Expected extinction rates as predicted by the SAR were highest for the geophilous species, probably due to their low dispersal capabilities. The comparatively low extinction rates predicted for the xylophilous species can be explained by their higher dispersal power, which reduces habitat confinement and hence makes them less sensitive to area loss, despite their dependence on the presence of forest fragments. The very low extinction rate predicted for endemics may seem like good news, but this result must be regarded with great caution, because the SAR for these species was poorly fitted. When other variables besides area were considered in multiple models, extinction rates following area reduction decreased and became more similar among groups. This suggests that considering only area may inflate expected extinction rates; when other factors are introduced, the role of area in explaining species richness may diminish, and hence the expected extinction rates following area loss. This should be considered when using the SAR to predict extinction rates to avoid overestimates. The strong negative effect of latitude highlights a prominent role of the southern areas for tenebrionid conservation in Italy. Thus, area loss is expected to act more dramatically in southern regions.

Conclusions
Area was an important determinant of overall tenebrionid species richness in Italian reserves, as well as for geophilous and widespread species, but not for xylophilous and endemic species. This indicates that focusing on reserve size in conservation planning may be important for generalist and widespread species but is not the best approach to protect insects with more specialised ecology and restricted distribution. The inverse relationship between species richness and latitude found in this study conforms to the well-known northward decrease in species richness recorded in Europe for a variety of taxa and can be explained by both the thermophilic preferences shown by most tenebrionids and the role of southern Italy as a refugial centre during Pleistocene glaciations. This pattern, however, was not found in xylophilous tenebrionids, which include many species associated with mesophilic forests. This is confirmed by the fact that minimum elevation had a negative influence on tenebrionid richness, except for xylophilous species, which are mainly associated with forests that are mainly distributed at middle elevations. Moreover, xylophilous species were positively influenced by elevational range, since a higher elevational range allows the presence of more forms of vegetation. Thus, the response of tenebrionids to area, latitude, and elevation did not follow the expected biogeographical patterns in all cases but depended on species' ecology. SAR-based extinction rates reflect species dispersal capabilities, being highest for the geophilous species (which area mainly flightless) and lower for the xylophilous species. These results highlight the importance of reserves in Southern Italy for the conservation of Mediterranean and South European insects, especially those with limited dispersal capabilities and with thermophilic preferences. In this regard, conservation of open habitats at low and middle elevations can be particularly important, not only for tenebrionids, if the results of the present research are supported by studies involving other animal groups. Forest environments are important for conservation of xylophilous insects through the Italian peninsula, even at the highest latitudes, and especially on mountains. Extinction rates based on multiple models that incorporated the effect of other variables besides area were lower than those based on area alone. This suggests that the use of area alone may overestimate extinction rates, when other factors (such as latitude and elevation) exert an important role in determining species richness. Further research involving other taxa in other regions would shed light on the generality of this pattern.
Funding: This research received no external funding.