A higher taxonomic richness does not ensure the functional resilience of saproxylic beetle communities in evergreen Quercus forests

1. Despite the fact that insects are suffering a global decline, little is known about the extent to which species loss affects functional diversity. Thus, to understand the relationship between taxonomic and functional diversity metrics, we focused on saproxylic beetles, which perform vital functions in forest ecosystems.


Introduction
Insects are facing a significant decline at a planetary level (Homburg et al., 2019;Sánchez-Bayo & Wyckhuys, 2019). The ecological consequences of this drop, however, are still different facets of insect diversity, such as taxonomic and functional dimensions, are related (Santini et al., 2017).
Many studies about functional diversity in plant communities have linked specific functional traits to ecosystem functions and services, thus allowing improvement management and conservation practices (see references in Drenovsky et al., 2012;Lozanovska et al., 2018;van der Plas, 2019). Conversely, the current state of knowledge of the functional diversity of insects seems to be still insufficient to fulfil such purposes. This has resulted in a poor understanding of the relationship between taxonomic diversity and functional diversity, which may be key to defining the consequences of species loss for ecosystems, an issue that is yet to be resolved (Cadotte et al., 2011). Numerous studies on insects have found that perturbations or habitat fragmentation can alter both ecosystem functioning and taxonomic diversity equally, because coexisting organisms have different functional characteristics (Audino et al., 2014;Luiza-Andrade et al., 2017;González et al., 2018;Martello et al., 2018). However, other studies have found different trends in taxonomic diversity and functional diversity, linking this decoupling to the functional identity of the species (Correa et al., 2019;Santoandré et al., 2019). That is, the stability and resilience of ecosystem services in the face of loss or changes of species may be more influenced by the species' function than their numbers (Tilman et al., 1997;Oliver et al., 2015). Moreover, species loss in natural communities is not random; it can depend on a range of factors such as population size or sensitivity to environmental stress Engel et al., 2017). All this makes it necessary first to analyse the relationship between different biodiversity metrics, such as taxonomic and functional measurements, so as to help calibrate the extent to which a large number of species may ensure ecosystem functions or a small number may jeopardise them. And second, we need to understand how species loss (both random or directed) can affect functional metrics, in order to determine how ecosystem functioning may react to future perturbations or changes (Engel et al., 2017;Leclerc et al., 2020).
One target group that allows to analyse the relationship between functional and taxonomic diversity metrics is saproxylic beetles, which depend on dead wood during some part of their life cycle (Speight, 1989;Alexander, 2008). Not only do they make up about one-third of the world's insect diversity (Ulyshen & Šobotník, 2018), but they also are a highly diverse, complex, and functionally important group for forest ecosystems. In this way, saproxylic beetles contribute to several key ecosystem functions in forests, as they play a significant role in the maintenance of nutrient cycling and trophic complexity. In addition, saproxylic predators naturally control forest pests and some adults are pollinators (Speight, 1989;Micó et al., 2011;Sánchez-Galván et al., 2014). Due to the significance of saproxylic insects, and particularly of their ecosystem services, many attempts have been made to understand the environmental factors that affect their biodiversity. Most studies have focused on taxonomic diversity (e.g. Müller & Bütler, 2010;Bouget et al., 2012;Quinto et al., 2014;García-López et al., 2016), whereas little is known about the relationship between the taxonomic and functional diversity of saproxylic beetles or how species loss can affect functional metrics.
Within this framework, we used the data of saproxylic beetle communities of three protected areas in the Iberian Peninsula dominated by mature evergreen Quercus species. The objective was two-fold: on the one hand, to assess the relationship between the taxonomic and functional diversity patterns, and on the other, to test the sensitivity of saproxylic beetle functional diversity to possible species loss (Petchey & Gaston, 2002;Cadotte et al., 2011;Parisi et al., 2018).
Taxonomic diversity may not necessarily be a reliable surrogate for the functional diversity of saproxylic groups, so both approaches should be used as complementary tools in conservation studies. In fact, since saproxylic beetles are strongly affected by environmental filters relating to habitat and microhabitat (e.g. Okland, 1996;Gibb et al., 2006;Micó et al., 2015;García-López et al., 2016;Kissick et al., 2018;Parmain & Bouget, 2018), we believe that the various taxonomic and functional diversity metrics would not necessarily change in the same direction. For example, communities in a preserved area, with a small number of species, may be able to maintain their functional diversity, while at the same time, other metrics, such as functional redundancy or rarity, may reflect the vulnerability of these communities to future changes (Scheffer et al., 2015;Kissick et al., 2018;Martello et al., 2018;Thorn et al., 2018). Moreover, we expect that the loss of species (whether random or directed) would affect the functional richness and resilience of communities in different ways, since species can have different functional roles in communities (Mason et al., 2005;Carmona et al., 2017;Violle et al., 2017). In addition, we anticipate that the analysis of trait patterns would help to better understand the diversity of each assemblage and to identify potential vulnerabilities regarding preservation functions (Watts & Mason, 2015;Micó et al., 2020).
Achieving our objectives would allow us not only to highlight the importance of using both taxonomic diversity and functional diversity when proposing appropriate conservation measures, but also to detect the vulnerability of functionally relevant communities, such as saproxylic beetles, even in protected areas, despite the number of species.

Study area
Evergreen Quercus forests were surveyed in three protected areas in the Spanish Mediterranean region separated by a distance of between 200 and 480 km ( Fig. 1): Cabañeros National Park (hereinafter, Cabañeros), Sierra Espadán Natural Park (hereinafter, Sierra Espadán), and Font Roja Natural Park (hereinafter, Font Roja). The dominant trees in all three parks are evergreen mature Quercus species; a description of the characteristics of each park is given in Table 1.

Sampling strategy
Fifteen sites were selected within the three protected areas, which best represented the vegetation's structural variability within each area (see Table S1 for details). In Sierra Espadán, three sites (separated by a maximum distance of 6 km) in Quercus suber forests were selected: two of the sites were located in a valley and had a highly complex undergrowth; and in the third, the undergrowth was sparser, and the cork oaks undergo occasional management.
In Font Roja, the three selected sites in Quercus rotundifolia forests were located on an altitudinal gradient (from 1087 to 1275 m) and separated by a maximum of 1 km. They presented the same forest composition with homogeneous tree size, shallow, and sun exposure, and little undergrowth.
Cabañeros presented the highest woodland heterogeneity because of variations in tree species composition and had both open and closed landscape structures. So, we selected nine sites, separated by a maximum distance of 24 km: three in a Q. suber forest with the trees of similar size, undergrowth, and deadwood in soil; and six in a Q. rotundifolia forest, of which four were in 'la Raña', typified by savanna-type open vegetation with sparse undergrowth, and two were on a private land with denser vegetation.
Since the samples were collected in protected areas, three flight interception traps were placed on each site, which ensures the minimum number of replications. Thus, nine traps were installed in both Sierra Espadán and Font Roja, and 27 traps were set up in Cabañeros. All the traps were kept in the field for one year, from April 2015 to June 2016. The collector bottle in each trap was replaced monthly, using propylene glycol as a preservative. The beetle specimens were stored in alcohol until their species were identified, with the collaboration of several taxonomists (see Acknowledgements). The nomenclature was based on the Fauna Europaea (http://www.faunaeur.org/), the study by Bouchard et al. (2011), and the Catalogue of Palaearctic Coleoptera (Löbl & Smetana, 2004, 2006, 2008, 2010. The specimens were added to the Entomological Collection of Alicante University (CEUA) and deposited at the CIBIO Research Institute (University of Alicante, Spain).

Functional trait selection
Four different types of functional traits (morphological, trophic, phenological, and a surrogate for physiological) were chosen to characterise the assemblages (Table 2).
Of the five morphological traits, three were directly measured for each specimen (body length, antenna length, and eye surface), whereas two (robustness and elytra ratio) were calculated based on a combination of other morphological measures (Table 2, Fig. 2). The morphological measurements were taken from up to 11 individuals per species, when possible, to capture within-species variations (Fountain-Jones et al., 2015). A stereomicroscope (Leica M205C) and leica Application Suite software version 4.8 were used for the measurements. After calculating each trait's mean per species, it was standardised using body size, following the study by Barton et al. (2011). The single phenological trait (the beetles' active period or active months) was calculated as the number of months that the adult beetles of each species were captured; it was thus used as a surrogate for the number of months that the beetles were active (Moretti et al., 2017). The trophic trait was the abundance of each species in the feeding guild (Watts & Mason, 2015;Moretti et al., 2017; Table 2).
Finally, regarding the surrogates for the physiological traits, we measured the parks' air temperature (Table 2). We installed one data logger (HOBO U23 Pro v2 Temperature/RH Data Loggers with U23-001 sensor) in each site in Font Roja and Sierra Espadán. In Cabañeros, the temperature data were obtained from the 'Alcornoquera' meteorological station located inside the park, which belongs to the Global change Monitoring Network in National Parks of Spain's Ministry for Agriculture, Food and Environment (MAPAMA, 2017).
In order to simplify our analytical model, we performed three successive principal component analyses (PCAs), removing the previously discarded traits at each new analysis. The first was performed using only the residuals of the morphological traits. Its first two axes explained nearly 99% of the variation, body length having the main coefficient, whereas the second axis captured the rest of the morphological traits (Table S2). For the second PCA, to which we added the surrogates for the physiological traits, the first three axes explained nearly 90% of the variance and, therefore, we decided to keep TRangeAct and TmeanAct (see Table 2 for the definitions of the abbreviations) as independent physiological traits (Table S2). In the last PCA, which included the added phenological trait, the first three axes explained nearly 96% of the variance. Since TRangeAct and active months were closely correlated, we decided to remove TRangeAct (Table S2). The final traits that we used for the rest of the analysis were body length, robustness, elytra ratio, both sensory traits (body length-corrected eye area and antenna length as residuals), TmeanAct, active months, and all the trophic guilds (Table 2).

Data analysis
To verify the inventory completeness of our sampling, both for each park and for all the parks together, we calculated the sample coverage estimator, which measures the proportion of the total number of individuals in a community that belongs to the species represented in the sample (Chao & Jost, 2012). We used the 'SPADE' package in r version 3.4.3 (Chao & Shen, 2010;Chao et al., 2016;r Core Team, 2017).
To analyse the differences, between the parks, in the patterns of diversity metrics for taxonomic diversity and functional diversity, we calculated two and six metrics for each park, respectively. Due to the data's non-normality, we looked for  (Forsythe, 1987;Ribera & Nilsson, 1995;Cunningham & Murray, 2007;Barton et al., 2011;Fountain-Jones et al., 2017;Gillespie et al., 2017) Antenna length* Maximum length from the base of the antenna to its apex (residual when body length-corrected) Both sensory traits can vary between trophic or taxonomic groups (Fountain-Jones et al., 2015) and they may indicate • Microhabitat use and structure • Lifestyle of species • Ability to locate new microhabitats or prey (Ribera et al., 1999b;Woodcock et al., 2010;Talarico et al., 2011) Eye surface (mm 2 )* Area contained within the outline of the eye perimeter (residual when body length-corrected) Robustness* Combination of the body length-corrected maximum dorsal width of the head and the body length-corrected maximum width, length, and lateral depth of the pronotum (Barton et al., 2011;Micó et al., 2020;Pérez-Sánchez et al., 2020) It can be used as • An indicator of the microhabitat • As a protective trait (Barton et al., 2011;Fountain-Jones et al., 2017) Elytra ratio* Body length-corrected elytra ratio: ratio of maximum lateral length of the elytra to the measure of maximum lateral length from meso-metathorax to the apex of the abdomen It measures the extent to which the elytra cover the abdomen (lower ratios mean exposed abdomen) It may provide information about • Flight capacity of the beetles (in the Carabidae short elytron length is linked to better dispersal efficiency) • Microhabitat use (longer elytra may provide protection for the wings in harsh environments; Forsythe, 1987;Ribera et al., 1999a;Barton et al., 2011;Johansson et al., 2012;Fountain-Jones et al., 2015) Trophic Difference between the environmental maximum and minimum temperatures of the months when the adult beetle was active ( ∘ C) TmeanAct* Environmental mean temperature of the months when the adult beetle was active ( ∘ C) WSD Weighted mean of the environmental temperature of the months with maximum abundance of adult beetles * Final traits used for statistical analysis.
significant differences using the nonparametric Kruskal-Wallis test (Ostertagová et al., 2014;De La Riva et al., 2017). For taxonomic diversity, we calculated species richness (S) and the number of effective species using Hill numbers, or true diversity of order 1 ( 1 D), since it weights each species according to its frequency in the sample (Jost, 2006;Scheiner, 2012). True diversity was calculated as the exponent of the Shannon index using the 'SPADE' package in r (Chao & Shen, 2010).
For functional diversity, we calculated functional richness, functional richness after removing the effect of the number of species, functional evenness, functional divergence, functional redundancy (after removing the effect of the number of species), and functional rarity. Functional richness is the amount of trait space occupied by the species present within a community (Mason et al., 2005;Laliberté & Legendre, 2010). Functional evenness is the regularity of the distribution of species abundances in the occupied functional trait space Maximum dorsal width of the head and the body and maximum width, length, and lateral depth of the pronotum were used to create Robustness, whereas Elytra length and meso-metathorax and abdomen length were used to create the Elytra ratio. Information and definition of traits in Table 2. Illustration by Sara Pérez-Sánchez (Mason et al., 2005;Laliberté & Legendre, 2010). Functional divergence measures the degree to which species abundances in a community are distributed towards the edges of the occupied trait space. Functional redundancy refers to the diversity of species performing similar functions in the ecosystem and represents community resilience (Mason et al., 2005;Villéger et al., 2008;Laliberté & Legendre, 2010;Pillar et al., 2013;Carmona et al., 2016;McPherson et al., 2018). Finally, the functional rarity of species in each park was calculated as the combination of functional distinctiveness (whether a species is more or less functionally close to the rest of the community) and taxon scarcity (a local-scale feature of a species with low relative abundance in the community; Grenié et al., 2017;Violle et al., 2017). We used 'FD' (Laliberté et al., 2015), 'SYNCSA' (Debastiani & Pillar, 2012), and 'funrar' packages (Grenié et al., 2017) in r (R Core Team, 2017).
Since both functional richness and redundancy are strongly dependent on the taxonomic richness, we divided the functional redundancy by its upper bound (species richness -1) (Carmona et al., 2017), while we applied a null model for functional richness (Mason et al., 2013). We generated 999 randomised abundance matrices, through the 'trialswap' randomisation algorithm (Miklós & Podani, 2004). For each park, we calculated the standardised effect size (SES) as the difference between the observed functional richness and the mean of the 999 matrix simulations divided by the standard deviation of the simulated values (de Bello, 2012); the P-value was calculated as the proportion of expected functional richness values that were less or more extreme (i.e. ≥ or ≤) than the observed value (Götzenberger et al., 2016). For the null model, we used the 'randomizeMatrix' function of the 'picante' package in r (Kembel et al., 2010;r Core Team, 2017).
In order to determine the strength and direction of the relationships between the two taxonomic metrics and the five functional metrics (without considering the functional richness dependent on taxonomic richness), we performed Spearman rank correlation tests (Vieira et al., 2013;Morelli et al., 2018). They were conducted for each protected area separately and also for the three protected areas together (at a regional level).
In addition, to uncover how the loss of species (reduction in species richness) may influence the functional richness and redundancy of the communities for each protected area, we applied three species extinction scenarios, the first with a random loss of species. However, since extinctions in nature are not random, we also modelled a directed loss of species in the other two species extinction scenarios. To do this, we focused on the fact that rare species (understood as those represented only by few individuals or restricted to particular habitats -functionally rare species) are considered to be more vulnerable to loss and affected by changes or perturbations Violle et al., 2017). Therefore, for one of the directed scenarios, we considered the functional distinctiveness or rarity of each species (the uncommonness of a species' traits compared to other species' traits in an assemblage; Grenié et al., 2017;Violle et al., 2017), while for the second scenario, we considered the abundance of each species (by removing the less abundant species). Species extinction was applied to all scenarios by removing intervals of 5% (from 5 to 70% species loss) from each protected area's total richness. In addition, to detect at which point significant differences of functional richness and redundancy appeared between the loss intervals, we applied the Kruskal-Wallis test between the species loss intervals.
Finally, to analyse trait patterns between parks, the selected traits were first weighted according to the number of individuals per trap (community weighted mean; Watts & Mason, 2015). Due to the data's non-normality, we then looked for significant differences between parks using the nonparametric Kruskal-Wallis test (Ostertagová et al., 2014;De La Riva et al., 2017).

Results
The sample coverage estimator for the group of parks showed that our inventory completeness was 82.2% accurate. It was also high for each park, with 89.4% for Cabañeros, 82.2% for Font Roja, and 91.3% for Sierra Espadán. The total number of specimens collected during the sampling year in the three parks was 14 127. It spanned 288 species and 44 families (Table S3). Sierra Espadán presented the highest number of species (189 species), followed by Cabañeros (178 species) and Font Roja (94 species), with the lowest species richness. Of the 288 identified species, the 242 species (13 616 individuals, belonging to 43 families) for which all traits were measured were included in the analysis (Table S4).
Sierra Espadán presented a greater taxonomic richness than Font Roja and Cabañeros (Kruskal-Wallis: X 2 2 = 21.20; P < 0.001; Fig. 3). However, when using the true diversity of order 1 (X 2 2 = 5.36; P = 0.07), Sierra Espadán had the lowest values and Cabañeros the highest. Sierra Espadán was the only protected area to present a different functional richness (X 2 2 = 9.57; P = 0.008), with higher values than the other areas. However, when the effect of the number of species was removed from functional richness (using 'trialswap' randomisation), significant differences were found between Font Roja and Cabañeros (X 2 2 = 12.48; P = 0.002), Font Roja presenting the highest values and Cabañeros the lowest (Fig. 3). Moreover, the results from null models showed that the observed functional richness of Sierra Espadán was significantly lower than the expected functional richness (SES = −0.21, P = 0.05). Finally, although functional evenness (X 2 2 = 0.89; P = 0.64) did not show variations between parks, the rest of the metrics showed significant differences between Sierra Espadán and the other parks (Fig. 3). Sierra Espadán presented the highest functional divergence (X 2 2 = 10.22; P = 0.006) and functional rarity (X 2 2 = 15.54; P < 0.001), whereas an opposite pattern was found for functional redundancy (X 2 2 = 21.49; P < 0.001), Sierra Espadán having the lowest (Fig. 3).
Results of the relationship between the taxonomic and the functional metrics showed that at a regional level (considering all protected areas together), both species richness and true diversity were positively related to the functional richness and the SES.Functional Richness (after 'trialswap' randomisations), but negatively related to the functional redundancy (Table 3, Fig. 4). In addition, true diversity was also negatively related to functional rarity and functional divergence (Table 3, Fig. 4). Conversely, when analysing the relationship between the taxonomic and functional metrics for each protected park, the results changed, and only two relationships common to the three areas were significant: the negative association between the taxonomic richness and functional redundancy, and the negative association between true diversity and functional rarity (Table 3, Figure S1).
The three species loss scenarios at a regional level (across all three parks) showed the functional redundancy pattern, which suffered a strong decline in the first percentage of species loss (5%) (Fig. 5b). On the other hand, functional richness (using the 'trialswap' randomisation -SES.Functional Richness) hardly changed in the face of species loss, showing only a very slight increase in the scenarios with a species loss of over 50% (Fig. 5a).
In each park, functional redundancy showed the same pattern as at a general level, suffering a strong decline in the first percentage of species loss (5%) and then slightly increasing ( Figures S2-S4). In addition, SES.Functional Richness also showed the same patterns at a general level (without any abrupt change or with only a very slight increase when species loss was over 50%), except in the case of Sierra Espadán, for which SES.Functional Richness significantly increased when 5% of the species were lost (Figures S2-S4).  Table 3. Spearman rank correlation test results between (a) the taxonomic richness and the functional diversity metrics, and (b) true diversity and the functional diversity metrics, for each park separately and for the three parks together (regional).

Trait patterns
Most traits differed between the protected areas (Fig. 6). Regarding the elytra ratio (Kruskal-Wallis: X 2 2 = 10.10; P = 0.006), robustness (X 2 2 = 21.65; P < 0.001), and mean temperature (X 2 2 = 21.59; P < 0.001), Sierra Espadán and Font Roja had the highest and the lowest values, respectively, whereas Cabañeros presented intermediate values. This pattern was reversed for the residual antenna length (X 2 2 = 22.55; P < 0.001), Sierra Espadán presenting the lowest values and Font Roja the highest. In contrast, only Sierra Espadán, with the lowest values for body length, showed a significant difference regarding this trait (X 2 2 = 10.39; P = 0.006). Sierra Espadán had the highest values for residual eye surface (X 2 2 = 13.26; P = 0.001) and active months (X 2 2 = 6.31; P = 0.04), although these values did not significantly differ from that of Cabañeros. In the case of predators (X 2 2 = 25.89; P < 0.001) and the xylomycetophagous trophic guild (X 2 2 = 31.24; P < 0.001), the three parks presented variations: Sierra Espadán had the lowest abundance of predators but the highest of xylomycetophagous, and Font Roja had the highest abundance of predators but the lowest of xylomycetophagous. Sierra Espadán presented the lowest abundance of saprophagous beetles (X 2 2 = 16.56; P < 0.001), a significant difference with the other areas, whereas Font Roja had the highest abundance of beetles in the xylophagous trophic guild (X 2 2 = 13.67; P = 0.001). Finally, Cabañeros had the highest abundance of beetles in the saproxylophagous trophic guild and Sierra Espadán the lowest, this difference being significant (X 2 2 = 11.88; P = 0.002; Fig. 6).

Discussion
It is worth highlighting that although higher taxonomic diversity values generally support higher functional richness values, they can come with low functional redundancy values. In fact, our results show that the predominance of functionally rare or redundant species in a community can play a major role if a loss  Table 2 of species occurs, having an impact on the communities' functional vulnerability to perturbations or changes. In addition, we found differences between taxonomic and functional diversity metrics and trait patterns in each protected area, suggesting that various habitat and microhabitat factors strongly influence both facets of biodiversity. The comparison of the taxonomic and functional diversity patterns of each protected area showed that communities with small numbers of species (i.e. Font Roja) can present high functional diversity, whereas communities with large numbers of species (i.e. Espadán) may be less resilient because functionally rare species predominate. However, when analysing the relationships between taxonomic and functional metrics for each park and at a regional level (considering all the protected areas together), some common relationships were found. An increase in taxonomic diversity was related to an increase in functional richness (after randomisations), but to a decrease in both the functional redundancy and functional rarity. This indicates that although communities with higher numbers of species are using a greater portion of the niche space (Mason et al., 2005;Morelli et al., 2018;Pellissier et al., 2018), they may also be more highly specialised in the use of resources (Bihn et al., 2010;Farias & Jaksic, 2011;Morelli et al., 2018), for what the disappearance of species may have strong effects on the ecosystem functioning, making communities more vulnerable to changes or perturbations (Pillar et al., 2013;Ricotta et al., 2016;Violle et al., 2017).
Various species loss scenarios showed that although functional richness seems not to be affected by a species reduction at a regional level, the functional redundancy is strongly affected even when only 5% of species are lost, regardless of whether the species loss is random or directed (Fig. 5). This may be because functional richness is associated with a high number of redundant species, ensuring communities against species loss by exercising buffer functions (Winfree et al., 2015). However, this may also imply that functional redundancy will be strongly affected by species loss, causing a considerable decrease in the communities' functional resilience Carmona et al., 2016;Engel et al., 2017;Stotz et al., 2021). Conversely, when a community is dominated by functionally rare species with low functional redundancy, as in the case of Sierra Espadán, despite the fact that the communities may be generally more vulnerable, those rare species might be sustaining the functional diversity (Figures S2-S4; Mouillot et al., 2013;Kissick et al., 2018). Thus, when species are lost, a compensation effect may occur, increasing the community's functional richness and explaining a lower degree of the functional redundancy change (Larsen et al., 2005;Mouillot et al., 2013;Winfree et al., 2015).
According to all those results, a community with high values of taxonomic and functional richness together with high values of functional redundancy could be the best in maintaining ecosystem functions and buffering perturbations. However, saproxylic communities are very complex and are composed of numerous functional groups, which not only allow several key roles in ecosystems, but also respond differently to environmental changes (Johansson et al., 2007;Wetherbee et al., 2020). This makes it difficult to predict how communities are going to respond (Kozák et al., 2021). In this sense, our approach has allowed to detect saproxylic communities more fragile to changes or perturbations, as in the case of Sierra Espadán with lower values of functional redundancy and higher values of functional rarity. However, and as we have seen, the predominance of functional rare species may also mean that the functional diversity of that community could be better maintained in the face of a loss of species, avoiding a sharp drop in its resilience.

Functional trait patterns
Generally, saproxylic assemblages of each protected area showed different patterns of functional traits, some of which could be used to better understand what is underlying the taxonomic and functional metrics. In this sense, some functional traits can provide information about the main characteristics of the development of microhabitat species. For example, the smaller eyes and longer antennae, combined with less robust bodies found in Font Roja, may indicate the possible dominance of beetles exploiting confined microhabitats such as those inside wood or under bark ; larger eyes, shorter antennae, and robust bodies, as found in Sierra Espadán (Fig. 6), may indicate that the main microhabitats driving saproxylic diversity in this park are more open, e.g. tree hollows (Ribera et al., 1999b;Talarico et al., 2011;Micó et al., 2020). In addition, some other trait patterns could indicate possible perturbations within an area or a high sensitivity to change, as in the case of Font Roja and Sierra Espadán, respectively. Thus, communities with low dispersal ability (i.e. small body size and high relative length of elytra) may be more sensitive to changes or perturbations (Ribera et al., 1999b;Barton et al., 2011;Talarico et al., 2011;Johansson et al., 2012). However, the role of these traits in saproxylic beetle communities would need to be further explored, in order to properly determine their ecosystem services and to test to what extent they can help to predict vulnerabilities.
To conclude, both approaches (taxonomic and functional) should be used as complementary tools when making management and conservation decisions, since a higher number of species may not ensure the functional resilience of communities.
Stuben P., Verdugo A., Vienna P.P., Viñolas A., and Zapata J.L. In addition, we would like to thank the translation service of 'Centro Superior de Idiomas' of the University of Alicante for their help with the revision and correction of the language. Financial support was provided by the research Project 'Ministerio de Ciencia e Innovación y fondos EUFEDER' (CGL201231669), 'Ministerio de economía, industria y competitividad' (CGL2016-78181-R), and 'Generalitat Valenciana' (PROMETEO/2013/034 andAICO 2020/192). This research is part of Diana Pérez-Sánchez PhD studies, funded by 'Ministerio de Educación, Cultura y Deporte', through the fellowship FPU14/03721. The authors declare no conflicts of interest. [Correction added on 30-July, after first online publication: Author Contribution section has been updated]

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Information about the sites sampling in each protected area, including the name of the site, the tree species sampled in each site, the number of flight interception in each site, and the coordinates of each site. Table S2. Summary of principal component analysis performed on a) the residuals of the morphological traits; b) the morphological and surrogate of physiological traits; and c) the morphological and surrogate of physiological and the phenological traits. High coefficients are shown in bold. Table S3. List of identified species with the families and subfamily they belong, abundances (number of individuals), and trophic group of each one. Table S4. List of species used for the analyses with the families and subfamily they belong and number of measured individuals. Figure S1. Correlation graphs for each park (blue squares represent Sierra Espadán Natural Park; red triangles represent Font Roja Natural Park, and green circles represent Cabañeros National Park) between a. taxonomic richness and the rest of functional diversity metrics; and b. true diversity and the rest of functional diversity metrics. The results of Spearman rank correlation tests are shown in Table 3, only for significantly correlations (P < 0.05) linear trends are drawn. Figure S2. Boxplots representation of SES.Functional Richness (functional richness after removing the effect of the number of species), and functional redundancy as the loss of random species increases by 5 percent (from 0 to 70 % species loss) for each protected area (Cabañeros National Park, Font Roja natural park and Sierra Espadán natural park). Different letters indicate significant differences according to the Kruskal-Wallis test between intervals (P< 0.05). The central rectangle indicates the second and third quartiles, and the black line the median value. Red points indicate mean values and "whiskers" lower and upper quartiles. Black points indicate possible outliers Figure S3. Boxplot representation for each protected area (Cabañeros National Park, Font Roja Natural Park, and Sierra Espadán Natural Park) of SES.Functional Richness (functional richness after removing the effect of the number of species) and functional redundancy, considering the species loss scenarios of functionally rare species (intervals going from 0 to 70% species loss). Different letters indicate significant differences based on the Kruskal-Wallis test between intervals (P < 0.05). The central rectangle indicates the second and third quartiles, and the black line the median value. Red points indicate mean values and whiskers lower and upper quartiles. Black points indicate possible outliers. Figure S4. Boxplot representation for each protected area (Cabañeros National Park, Font Roja Natural Park, and Sierra Espadán Natural Park) of SES.Functional Richness (functional richness after removing the effect of the number of species) and functional redundancy, considering the species loss scenarios of the less abundant species (intervals going from 0 to 70% species loss). Different letters indicate significant differences based on the Kruskal-Wallis test between intervals (P < 0.05). The central rectangle indicates the second and third quartiles, and the black line the median value. Red points indicate mean values and whiskers lower and upper quartiles. Black points indicate possible outliers.