Iberian Protected Areas Capture Regional Functional, Phylogenetic and Taxonomic Diversity of Most Tetrapod Groups

Llorente-Culebras, Sonia; Molina-Venegas, Rafael; Barbosa, A. Márcia; Carvalho, Silvia B.; Rodríguez, Miguel Á.; Santos, Ana M. C.

doi:10.3389/fevo.2021.634653

ORIGINAL RESEARCH article

Front. Ecol. Evol., 26 July 2021

Sec. Biogeography and Macroecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.634653

This article is part of the Research Topic Roles and Implications of Functional Traits and Phylogenies to Characterize Refugia Under Increasing Climate Variability View all 8 articles

Iberian Protected Areas Capture Regional Functional, Phylogenetic and Taxonomic Diversity of Most Tetrapod Groups

$\r\nSonia Llorente-Culebras,$ Sonia Llorente-Culebras^1,2

Rafael Molina-Venegas^2,3

A. Márcia Barbosa⁴

Silvia B. Carvalho⁵

Miguel Á. Rodríguez^2,3

Ana M. C. Santos^*1,3,6

¹Terrestrial Ecology Group (TEG-UAM), Departamento de Ecologia, Universidad Autónoma de Madrid, Madrid, Spain
²Universidad de Alcalá, Department of Life Sciences, Alcalá de Henares, Spain
³Universidad de Alcalá, GLOCEE - Global Change Ecology and Evolution Group, Department of Life Sciences, Alcalá de Henares, Spain
⁴Spatial Biology Lab, CICGE, Universidade do Porto, Vila Nova de Gaia, Portugal
⁵CIBIO/InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Oporto, Portugal
⁶Centro de Investigación en Biodiversidad y Cambio Global (CIBC-UAM), Universidad Autónoma de Madrid, Madrid, Spain

Protected areas (PAs) have been created with the purpose of preserving biodiversity, acting as refuges from anthropogenic pressures. Traditionally, PAs have been designed and managed to represent mainly taxonomic diversity, ignoring other diversity facets such as its functional and phylogenetic components. Yet, functional and phylogenetic diversity are, respectively, connected with species’ roles on ecosystems and evolutionary history held within communities. Here, we focused on the amphibian, reptile, resident breeding bird, and non-flying mammal faunas of the national and natural parks of the Iberian Peninsula, to evaluate whether these PAs are adequately representing regional functional, phylogenetic, and taxonomic diversity of each group. Specifically, we computed functional and phylogenetic diversity within each PA, and then compared those values to the ones obtained from a random assembly of species from the regional pool, that was defined as the region encompassing the PA and a neighboring area of 50 km beyond its boundary. We also calculated the proportion of species in each regional pool that were present within the PAs. In general, the functional and phylogenetic diversity of amphibians, reptiles and non-flying mammals found within PAs did not differ significantly from random expectations generated from the species pertaining to the regional pool, although a few PAs showed a higher diversity. In contrast, resident breeding birds presented lower functional and phylogenetic diversity than expected by chance in many of the PAs, which could relate to climatic variables and the habitat specificity of some species. The proportion of species from the regional pools that are present in the PAs was high for amphibians, reptiles and mammals, and slightly lower for birds. These results suggest that the Iberian natural and national parks are effectively capturing the functional, phylogenetic and taxonomic diversity of most tetrapod assemblages present at the regional level. Future studies should identify priority areas to expand the representation of these biodiversity components, and assess potential effects of climate and land-use changes on current patterns.

Introduction

Human activity is causing serious ecological disturbances on Earth, which has led researchers to propose a new geological epoch, the Anthropocene (Waters et al., 2016). Land-use change, habitat loss and fragmentation, alien species, and climate change are prime drivers of ecosystem alterations (MEA - Millennium Ecosystem Assessment, 2005; Allan et al., 2019). At the global scale, the impact of these actions is significant and is inducing waves of species extinctions and population declines (Pimm et al., 1995; Hoffmann et al., 2010; Barnosky et al., 2011; Davidson et al., 2017; Ceballos et al., 2020), as well as disruption of ecosystem services (Isbell et al., 2017; Pecl et al., 2017; IPBES, 2019).

Protected areas (PAs) are legally protected geographic spaces that have been created to preserve biodiversity and mitigate the impact of human actions. These areas can potentially act as refuges, giving protection from ongoing temporal and/or spatial perturbations (Keppel et al., 2012). In the long-term, PAs may even act as refugia (sensu Keppel et al., 2012, 2018), if they are capable to maintain species occurrence when conditions become adverse in the surrounding areas (e.g., due to pollution or climate change). The design of PAs has been traditionally achieved through information on species distributions, being particularly based on the concepts of species representativeness and species protection from human impacts (Margules and Pressey, 2000), while accounting for economic and opportunity costs (Moilanen et al., 2009), but focusing mostly on taxonomic diversity (Myers et al., 2000; Brooks et al., 2006; Jenkins et al., 2013). By contrast, functional diversity (defined as the range, values, and relative abundance of functional features or traits of a given community; Díaz and Cabido, 2001; Díaz et al., 2007; Harrington et al., 2010) and phylogenetic diversity (i.e., the amount of evolutionary history encompassed by a set of taxa; Faith, 1992) have generally been left aside when designing PAs (Winter et al., 2013; Guilhaumon et al., 2015; Pollock et al., 2017). Yet, the emphasis on taxonomic diversity implies assuming that all species contribute to biodiversity equally, which ignores: (1) that species differ in terms of evolutionary distinctiveness (Faith, 1992; Rodrigues and Gaston, 2002); (2) that some ecosystem functions are performed by particular species (Flynn et al., 2011); and (3) that extinctions are often not random among species [e.g., across vertebrate groups, higher extinction probability has been reported for species from certain lineages such as carnivores or primates (e.g., Purvis et al., 2000), or with certain traits, such as species with a large or small body-mass, slow pace of life and/or a small distribution range size (e.g., Fritz et al., 2009; Ripple et al., 2017; Carmona et al., 2021)]. Therefore, it is essential to develop projects and proposals that integrate more facets of diversity and sources of data into conservation efforts (e.g., Devictor et al., 2010; Barak et al., 2016). One example comes from the EDGE – “Evolutionary Distinctive and Globally Endangered” project, which aims at prioritizing endangered and evolutionary distinct species in conservation planning (Isaac et al., 2007).

The effectiveness of PAs in preserving biodiversity has been traditionally evaluated with a focus on taxonomic diversity (e.g., Long et al., 2004; Araújo et al., 2007; González-Maya et al., 2015; Kukkala et al., 2016; Rosso et al., 2018; Hanson et al., 2020), while attention to functional and phylogenetic diversities is only beginning to arise (Devictor et al., 2010; Villamor and Becerro, 2012; Mazel et al., 2014; Guilhaumon et al., 2015; Thuiller et al., 2015; Brum et al., 2017; Pollock et al., 2017). Even so, since most of these recent studies have evaluated the representation of the gamma diversity hosted within the PAs network at broad geographical extents, how well each PA preserves all facets of biodiversity relative to their neighborhoods remains an open question (see Gaston et al., 2008). Critically, by assuring the maintenance and conservation of functionally and phylogenetically diverse communities–presumably better prepared to adapt to global change (Sgro et al., 2011; Mori et al., 2013)–we can, in general, also preserve the resilience of the ecosystems and the continuity of services they provide, particularly at larger spatial scales (Loreau et al., 2001; Díaz et al., 2013).

Focusing on the Iberian Peninsula (Iberia herein), comprising mainland Portugal and Spain, this study aims to evaluate whether each national and natural park is adequately representing the regional functional, phylogenetic, and taxonomic diversity across tetrapod assemblages (sensu Fauth et al., 1996; i.e., amphibians, reptiles, resident breeding birds, and non-flying mammals). With an extent area of 583,254 km², which corresponds to 5.7% of Europe’s area, Iberia stands out in this continent for hosting approximately half of Europe’s animal and plant species (Williams et al., 2000; Araújo et al., 2007). The effectiveness of Iberian PAs as a conservation tool has been evaluated before. Araújo et al. (2007) concluded that Iberian PAs reasonably represent the taxonomic diversity of the vertebrates and plants that inhabit this region, and Rosso et al. (2018) suggested that Natura 2000 Network is an effective representation of the Iberian endemic fauna. There are also some studies focusing on other biodiversity components. For instance, Villamor and Becerro (2012) evaluated fish functional diversity in Spanish coastal PAs, and Carvalho et al. (2017) identified priority areas for conservation of Iberian herpetofauna phylogenetic diversity.

Yet, no study to date has examined whether each Iberian PAs is safeguarding regional phylogenetic and functional diversity across the taxonomic groups considered in this study. Here, we used a null model approach to evaluate whether each Iberian PA holds non-random fractions of both functional and phylogenetic diversity of the region where they stand. Protected areas showing functional or phylogenetic diversity values that are not different, or are higher, than expected by chance could be deemed as adequate representations of the regional faunas, and those PAs showing lower values as poor representations. Additionally, we evaluated the proportions to which the species present in the regional pools were also present in their corresponding PAs. Finally, we also investigated whether PA’s climate and habitat characteristics were related with these proportions and/or with deviations from randomness in functional and phylogenetic diversity, to explore potential aspects influencing PA’s representativeness of regional diversity.

Materials and Methods

Study Area

We focused on the two most important protection categories of Iberian PAs, i.e., National and Natural parks, of which there are 137 in mainland Spain (9 and 114, respectively) and Portugal (1 and 13, respectively) (Supplementary Table 1 in Data Sheet 1). Polygon maps for these PAs were downloaded from the websites of the former Spanish Ministry of Agriculture and Fisheries, Food and Environment (MAPAMA - Ministerio de Agricultura, Pesca, Alimentación y Medio Ambiente, 2017b), and the Portuguese “Instituto da Conservação da Natureza e das Florestas” (ICNF - Instituto da Conservação da Natureza e das Florestas, 2017). Afterward, they were combined, merging adjacent parks from both countries to form single analysis units (this occurred for 40 parks; see Supplementary Table 1 in Data Sheet 1), thus rendering a total of 113 different PAs to be considered in this study (the PAs’ main characteristics can be found in Supplementary Table 1 in Data Sheet 1).

Species Data

Presences/absences of terrestrial tetrapod species native to Iberia at the PAs were obtained by overlapping PAs’ polygon maps with several species occurrence data sources available at the original Universal Transverse Mercator (UTM) 10 × 10 km grid. Thus, amphibian and reptile data were extracted from Loureiro et al. (2008) and the Spanish Herpetological Society database SIARE (SIARE (Servidor de Información de Anfibios y Reptiles de España) - Asociación Herpetológica Española, 2017), while bird and mammal data were obtained from the atlas of mammals (Bencatel et al., 2017) and breeding birds of Portugal (ICNB – Instituto da Conservação da Natureza y da Biodiversidade, 2008), and the Spanish inventory of terrestrial species database (MAPAMA - Ministerio de Agricultura, Pesca, Alimentación y Medio Ambiente, 2017a). Bats, non-resident and non-breeding birds were excluded, thus rendering 383 native species in our database, of which 26 were amphibians, 50 reptiles, 247 birds, and 60 mammals.

Trait data were compiled from the literature and were related to: trophic behaviour, habitat, body size, daily activity, reproduction, longevity, and displacement mode (Table 1). Due to the physiological and anatomical differences between tetrapod groups, some of the traits differed across groups, yet they represented common biological functions (see Table 1).

TABLE 1

Table 1. Functional traits compiled for amphibians, reptiles, birds, and mammals using different sources: 1-García-París et al., 2004; 2-Trochet et al., 2014; 3-Carvalho et al., 2017; 4-Grimm et al., 2014; 5-Salvador, 2014; 6-Cramp et al., 2004; 7-Dunning, 2008; 8-Wilman et al., 2014; 9-Hume, 2016; 10-SEO - Sociedad Española de Ornitología, 2018; 11-Palomo et al., 2007; 12-Aulagnier, 2009; 13-Jones et al., 2009; 14-Safi et al., 2011.

Regarding the phylogenetic data, for amphibian and reptile species, we used the species-level ultrametric phylogenies published in Carvalho et al. (2017). The trees were inferred using Bayesian methods (BEAST v.1.8.0) based on three mitochondrial markers – 12S rRNA (12S), 16S rRNA (16S), and the nuclear gene cytochrome b (CYTB). The amphibian tree was inferred using an additional nuclear marker – the recombination-activating gene 1 (RAG1) – and the reptile tree with three additional nuclear markers – the oocyte maturation factor Mos (CMOS), the melanocortin receptor 1 (MC1R), and the recombination-activating gene 1 (RAG1) (see Carvalho et al., 2017 for full details). For birds and mammals, the phylogenies were built de novo using CYTB, COI, 16S, and 12S mitochondrial markers and a maximum likelihood approach (Roquet et al., 2013). The resultant trees were time-calibrated using the TreePL software (Smith and O’Meara, 2012) and node age estimations from Meredith et al. (2011) and Prum et al. (2015; see Supplementary Appendix 1 in Data Sheet 1 for full details on the phylogenetic procedures). Since the phylogenetic inference approaches used for amphibians and reptiles are different from those used for mammals and birds, there could be slightly different topologies toward the shallower branches (relationships among species within genera). However, accumulating evidence suggests that genus-level phylogenies are appropriate for phylogenetic diversity analyses (e.g., Allen et al., 2019; Qian and Jin, in press), and therefore, the putative differences in among-species evolutionary relationships that may result from the use of different inference methods would have a negligible impact in the analyses. It is also important to notice that the height of the four chronograms was scaled to unit to make them comparable.

Diversity Metrics

To evaluate whether PAs capture the functional, phylogenetic and taxonomic diversity of the region where they stand, we first defined species regional pools for each of them. Four concentric areas were designed placing their limits 15, 25, 50, or 100 km beyond the boundaries of the PA. Then, we defined four nested regional pools of increasing size, each corresponding to both the list of species occurring within the PA and those occurring in each of the four buffer zones. However, as diversity metrics [i.e., proportion of species and Standard Effect Size (SES) values; see below] strongly correlated across regional pools (Supplementary Table 2 in Data Sheet 1), we only present results pertaining to the 50 km buffer area. We chose the latter because larger buffer areas overlapped between some PAs and smaller areas resulted in depauperated regional pools.

For each tetrapod group, we obtained a functional dendrogram including the species observed across all the Iberian PAs and buffers (hereafter “Iberian functional dendrogram”). To do so, the complete group’s trait matrix was converted into a distance matrix using a modified version of Gower’s distance (Pavoine et al., 2009; “dist.ktab” from package “ade4”; Dray and Dufour, 2007), which was in turn used to obtain the Iberian functional dendrograms (“hclust” from package “ade4”, using UPGMA as the agglomerative method; Dray and Dufour, 2007). This methodology has been widely used in the literature in studies including similar goals to ours (e.g., Roa-Fuentes et al., 2019; Santos et al., 2020). Then, the Iberian functional dendrogram of each tetrapod group was pruned to each regional pool (composed by the species observed in each PA and respective buffer zone). These PA-specific functional dendrograms served to compute functional diversity (FD) of the PAs and the corresponding null assemblages of species (see below) as the minimum spanning path connecting the set of species in each case (Petchey and Gaston, 2002). Although the computation of Gower’s distance allows missing data, trait completeness was overall high across all groups (97.9%, 97.8%, 99.9%, and 100% completed for amphibians, reptiles, birds, and mammals, respectively). Likewise, Faith’s (1992) Phylogenetic Diversity (PD) index was used to calculate the phylogenetic diversity of each PA and the null assemblages of the respective regional pool, following the same procedure described above for FD but using the molecular phylogenies instead of the functional dendrogram (see Santos et al., 2020 for a similar approach).

We assessed whether the set of species of each taxonomic group inhabiting each PA represents its regional functional and phylogenetic diversity by comparing the FD and PD observed in the PAs with a null distribution generated by the random assembly of species from their corresponding regional pool (i.e., PA plus buffer zone). For each PA, we created 1,000 random assemblages from its regional pool, with the same number of species as in the PA, by shuffling species labels across the tips of the corresponding functional and phylogenetic dendrograms (Supplementary Figure 1 in Data Sheet 1). We performed 1,000 randomizations as these provide very similar FD and PD values as those obtained when using a higher number of randomizations, allowing to compromise between the robustness of the results and optimization of computational time (see Supplementary Figures 2–5 in Data Sheet 1). This was achieved using the function “ses.pd” of the “picante” R package (Kembel et al., 2010) and applying the null model “taxa shuffle” to calculate standardized effect sizes (SES; Gotelli and Graves, 1996; Kembel, 2009) for FD (herein SES.FD) and PD (herein SES.PD) as:

S E S . F D = \frac{F D_{o b s} - m e a n F D_{n u l l}}{s d F D_{n u l l}};

and S E S . P D = \frac{P D_{o b s} - m e a n P D_{n u l l}}{s d P D_{n u l l}}

where the subscripts obs and null refer, respectively, to the observed diversity and the randomly generated null distribution of diversity values. For each SES, a p-value was calculated by placing the observed diversity in the ranked null distribution of randomized values and dividing its position by the number of randomizations +1. Thus, using a 5% nominal alpha, a positive SES.FD (or SES.PD) with p-value > 0.975 corresponds to a PA with significantly higher FD (or PD) than expected by a random assembly from the regional pool, while a negative SES.FD (or SES.PD) with p-value < 0.025 indicates the opposite. All values obtained through this methodology (observed FD/PD, SES.FD/PD, and p-values) can be found in Supplementary Table 1.

Regarding taxonomic diversity, we evaluated how much of the regional species diversity is present within each PA, by calculating the proportion of species from the regional pool that are also present in the PA (number of species present in the PA/number of species present in the total regional pool; herein species proportion). The values vary from 1 when all the species of the regional pool are present in the PA, and decrease as the number of species present in the PA is lower in comparison with its associated region.

Habitat and Climatic Data

Habitat diversity is one of the most important determinants of not only species richness (e.g., Tews et al., 2004; Moreno-Rueda and Pizarro, 2008), but also of functional diversity (García-Llamas et al., 2019) and phylogenetic diversity (Franke et al., 2020). In order to explore if there is any relationship between the habitat types and the regional diversity represented in the Iberian PAs, we compiled Iberian land uses from the CORINE Land Cover inventory, using them as proxies for habitats (for simplification, the different land covers will be referred to as habitats). This inventory classifies European land uses into three nested levels, from which we selected the second level with 15 categories (see details in Supplementary Table 3 in Data Sheet 1). We downloaded the 2012 update of these data (CORINE Land Cover, 2019) and processed them with Quantum GIS (version 3.6.0; QGIS Development Team, 2018) to generate a set of habitat variables for each PA: area of each CORINE land-use type, number of habitats, and diversity of habitats (calculated using the Shannon-Wiener index). We performed a preliminary analysis based on a paired Student t-test, to assess if the habitat configuration within the PAs (i.e., habitat richness, habitat diversity and the proportion of area occupied by each habitat) differs from what occurs in their surrounding region (see more details in Supplementary Appendix 2 in Data Sheet 1). Overall, this comparison indicated significantly different habitat configurations between inside and outside the PAs (see Supplementary Table A1 in Supplementary Appendix 2 in Data Sheet 1).

Climatic factors, as temperature and precipitation, influence vertebrates’ distribution and diversity (Moreno-Rueda and Pizarro, 2008), affecting particularly the environmental suitability for physiologically limited species (i.e., ectotherms; Clarke and Gaston, 2006). To explore the effect of climate on diversity metrics, we obtained climatic data for the 1970 to 2000 range from the WorldClim 2.1 database (WorldClim, 2020) at a resolution of ∼5 km × 5 km. We selected several climatic variables: Annual Mean Temperature, Temperature Seasonality, Minimum Temperature of Coldest Month, Annual Precipitation, Precipitation of Driest Month, and Precipitation Seasonality (see details in Supplementary Table 1 in Data Sheet 1). We processed these variables with Quantum GIS (version 3.6.0; QGIS Development Team, 2018) and then we obtained mean values of each climatic variable for each PA, by calculating the average values of the cells included within its boundaries. A preliminary analysis based on a paired Student t-test indicated there are no significant differences in the climatic conditions found within the PAs and their surrounding regions (see more details in Supplementary Appendix 2 in Data Sheet 1). However, as climate is known to have important effects on vertebrates’ functional and/or phylogenetic diversity, both at continental (Santos et al., 2020) and smaller scales (García-Llamas et al., 2019), we considered important including these variables in the analyses (see below).

Model Selection

To explore which PA’s characteristics might be related to deviation from randomness in functional and phylogenetic diversity, and also with species proportions of each taxonomic group, we used Generalized Linear Models (GLMs) in a multi-model inference procedure. This approach provides the best models (among all possible models) that explain the variation of the different diversity metrics between PAs.

We evaluated 12 response variables, resulting from the combination of taxonomic groups (amphibians, reptiles, birds, and mammals) and diversity metrics (functional diversity, phylogenetic diversity, and species proportion). The response variables for each vertebrate group were SES.FD, SES.PD, and the proportions of species from the regional pools that are present in PAs. The initial candidate predictor variables were related to three categories: habitat (area of each of the fifteen land use categories of CORINE Land Cover, number of habitat types, and habitat diversity of each PA), climate (mean of six bioclimatic variables; explained above), and the area of the PAs. All predictors were standardized to mean = 0 and standard deviation = 1. We then identified pairs of candidate predictors that were strongly correlated (|Spearman’s ρ| ≥ 0.8) and used the “corSelect” function of the “fuzzySim” R package (Barbosa, 2015) to exclude the variable with the weakest relationship with each particular response variable.

The remaining predictors were included in a multi-model inference procedure based on Akaike’s Information Criterion (AIC; Akaike, 1973) as implemented in the “glmulti” R package (Calcagno and de Mazancourt, 2010), using Gaussian family with “identity” link when the SES values were the response variables, and binomial family with “logit” link in the case of the proportion of species from the regional pool. We first selected those models with lower AICc values, i.e., AIC values corrected for small sample sizes (see Symonds and Moussalli, 2010), identifying the best models as those with ΔAICc < 2. Afterward, we applied a model averaging procedure to these models using the function “coef.glmulti” of the “glmulti” R package (Calcagno and de Mazancourt, 2010). This function provides the averaged estimates of all the variables which appear at least in one of the best models (presented in Supplementary Table 2), as well as the importance of each predictor variable. The importance corresponds to the sum of the AICc weights of the models in which each predictor variable appears. The importance varies from 1, when a variable appears in all the best models, decreasing to 0, when a variable appears in a reduced number of models or none. For simplification, we only discuss those variables that have an importance value of 1 for each response variable, but the importance values of the other variables can be found in Table 2. The amount of deviance that the best model for each response variable accounted for (D²; the equivalent to R² of least square models, for general linear models; Guisan and Zimmermann, 2000) was calculated with the “Dsquared” function of the “modEvA” package (Barbosa et al., 2013) and is presented in Supplementary Table 2, along with the equation of the best models.

TABLE 2

Table 2. Importance value of each predictor after the model averaging of the most substantially supported GLM models (i.e., with ΔAICc < 2) quantifying how well protected areas (PAs) represent the functional (FD) and phylogenetic (PD) regional diversities of four tetrapod groups (see Materials and Methods).

All statistical analyses were performed in R Version 1.0.136 (R Core Team, 2017). The R code used can be found in Supplementary Data Sheet 2).

Results

Amphibians’ functional and phylogenetic diversity of most PAs did not differ from the diversity of the regional pool (Figures 1A, 2A; see also Supplementary Figures 6A, 7A in Data Sheet 1). The exceptions are seven PAs from northern Iberia that presented significantly higher functional diversity than expected from the regional pool (i.e., with significantly positive SES.FD values), and one PA with significantly high phylogenetic diversity (i.e., with a significantly positive SES.PD value; see Figures 1A, 2A). Environmental GLM models for both functional and phylogenetic diversities were generally weak for this group (D² = 0.17 and 0.15, respectively), but included a few environmental predictors with high importance values in each case. Related to this, our model selection procedure detected that SES.FD values tended to decrease in PAs presenting more seasonal climate (see Table 2). Moreover, SES.PD values tended to decrease in PAs with open spaces and scarce vegetation, and areas with marine waters, and to increase along with the cover of arable land. Regarding amphibian species proportions, more than half of the PAs (64 out of 113) contained at least 80% of the species present in the corresponding regional pools, of which 10 presented the same number of species as in the region where they stand (Supplementary Figure 8A in Data Sheet 1). Additionally, environmental predictors modeled did not explain the variability of amphibians species proportion (D²≈0 in the best model; see Supplementary Table 2).

FIGURE 1

Figure 1. Representation of the protected areas according to the level of significance of the standardized effect sizes of functional diversity for different tetrapod groups: (A) amphibians, (B) reptiles, (C) birds, and (D) mammals. Colours indicate whether each protected area holds significantly higher (red) or lower (blue) functional diversity than it is expected from the region where it stands (an area designed placing its limits 50 km beyond the boundaries of the PA, including also the PA); parks in grey correspond to those where functional diversity does not differ from a random draw of species from the region where they stands (see main text for more details).

FIGURE 2

Figure 2. Representation of the protected areas according to the level of significance of the standardized effect sizes of phylogenetic diversity for different tetrapod groups: (A) amphibians, (B) reptiles, (C) birds, and (D) mammals. Colours indicate whether each protected area holds significantly higher (red) or lower (blue) phylogenetic diversity than it is expected from the region where it stands (an area designed placing its limits 50 km beyond the boundaries of the PA, including also the PA); parks in grey correspond to those where phylogenetic diversity does not differ from a random draw of species from the region where they stands (see main text for more details).

In the case of reptiles, again neither functional nor phylogenetic diversity of most PAs differed significantly from the regional pool (Figures 1B, 2B; see also Supplementary Figures 6B, 7B in Data Sheet 1). Only two PAs presented significantly high SES.FD values, while five PAs showed the opposite trend, having significantly lower SES.FD values than expected by chance from the regional pool (Figure 1B). On the other hand, eight PAs presented significantly high SES.PD values, while three PAs showed the inverse trend (Figure 2B). Environmental GLM models for reptiles explained more variation in SES.PD values (D² = 0.29) than in SES.FD values (D² = 0.16), a pattern also observed for birds and mammals (see Table 2). According to these models, SES.FD values tended to increase in PAs presenting high mean annual temperature, while SES.PD values were higher in habitat-rich PAs with reduced precipitation and scarce forest cover (variables which showed the highest importance values; see Table 2). Regarding reptile species proportions, only two PAs included all of the species present in the regional pool, and more than half of the PAs (66 out of 113) contained at least 75% of the regional taxonomic diversity (Supplementary Figure 8B in Data Sheet 1). As in the case of amphibians, we did not find a relationship between species proportion and the environmental predictors modeled (D²≈0 in the best model; Supplementary Table 2).

Contrary to other groups, birds exhibited lower functional and/or phylogenetic diversity than expected in their regional pool in a large number of PAs, or, in other words, they present significantly low SES.FD and/or SES.PD values (Figures 1C, 2C; see also Supplementary Figures 6C, 7C in Data Sheet 1). Specifically, 53 and 26 PAs showed significantly low SES.FD and SES.PD values, respectively. Among these PAs, 18 presented significantly low SES values for both facets of diversity (Figures 1C, 2C). Environmental GLM models for birds detected a trend for SES.FD values (D² = 0.19) to increase in PAs with higher habitat diversity and winter temperatures, and reduced precipitation seasonality (variables showing the highest importance values; see Table 2), and for SES.PD values (D² = 0.37) to increase in warmer habitat-rich PAs, with inland wetlands (Table 2). It is remarkable that the species proportion of birds in the PAs was the poorest of all groups investigated, with about half of the PAs studied (55) encompassing less than 65% of the regional bird diversity, and none having more that 85% of the species present in the regional pool (Supplementary Figure 8C in Data Sheet 1). Species proportion was also unrelated with the environmental variables considered (D²≈0 in all cases; Supplementary Table 2).

For mammals, we found once again a majority of PAs not differing from their regional pools in terms of functional and phylogenetic diversity (Figures 1D, 2D; see also Supplementary Figures 6D, 7D in Data Sheet 1). The exceptions in this case corresponded to five PAs that showed significantly low values of functional diversity and seven PAs which showed significantly high values of phylogenetic diversity (Figures 1D, 2D). Mammal models were the most explicative ones for each diversity facet (cf. D² values in Table 2), which suggests a stronger influence of environmental factors on this group. In these models, both SES.FD and SES.PD values tended to increase at PAs with lower annual mean temperature and higher annual precipitation. Also, SES.FD values increased at PAs with more arable land, whereas SES.PD values increased at PAs with more temperature seasonality and decreased with the area occupied by inland wetlands and by mine, dump and construction sites within the PAs (Table 2). On the other hand, 65 PAs comprised over 75% of mammal species present in the regional pool, of which only two encompassed 100% of all regional species (Supplementary Figure 8D in Data Sheet 1). As for the rest of the considered groups, the species proportion of mammals was also unrelated with the environmental factors considered (D²< 0.05 in all cases; Supplementary Table 2).

Finally, it has to be noted that, across all tetrapod groups and GLMs evaluated, PAs area only appeared as a potentially relevant factor for reptile SES.FD values, and in this case with a very low importance value (0.06) (see Table 2). This reinforces the notion that the patterns we found were not artifacts due to area effects.

Discussion

The effectiveness of the Spanish and Portuguese networks of PAs for safeguarding biodiversity has been evaluated several times through different approaches (e.g., Araújo, 1999; Carrascal et al., 2003; Araújo et al., 2007; Hernández-Manrique et al., 2012; Villamor and Becerro, 2012; Guilhaumon et al., 2015; Hermoso et al., 2015; Rosso et al., 2018; Estrada and Real, 2018; Rodríguez-Rodríguez and Martínez-Vega, 2018). Yet, the ability of each area to represent functional and phylogenetic diversity components of the regional pool remained unknown. Here, we found that, except for birds, Iberian national and natural parks are generally adequate, and in some exceptional cases even optimal (i.e., better than expected by chance) representations of the regional functional and phylogenetic diversities. Iberian PAs also generally preserve a high proportion of the regional species richness. These results complement previous findings indicating that, overall, the existent network of PAs in this area adequately represents diversity for some taxonomic groups (e.g., Araújo et al., 2007; Rosso et al., 2018), although this does not necessarily hold for all the groups evaluated so far (e.g., Araújo et al., 2007; Sánchez-Fernández et al., 2008; Hernández-Manrique et al., 2012; Hermoso et al., 2015).

Previous studies focusing on European PAs showed that these achieve a significant representation of functional distinctive amphibian diversity (Thuiller et al., 2015). These results are in line with our findings, as the functional diversity of most PAs did not differ from that of the regional pool, and seven PAs located in northern Iberia even presented higher functional diversity than expected from the regional pool. Interestingly, some of these PAs match the main hotspots of amphibian richness of this region (Rey-Benayas and de la Montaña, 2003; Martins et al., 2014). As most amphibians are highly dependent on humid habitats to complete their reproductive cycle (Arnold and Burton, 1982; Vences and Köhler, 2007; da Silva et al., 2012), it is likely that the generally wetter climate of northern Iberia, which is maintained throughout the year, favors a functionally diverse amphibian fauna (Rey-Benayas and de la Montaña, 2003; Martins et al., 2014), a notion also supported by the decrease we found in amphibian’s functional diversity within the PAs presenting higher climatic seasonality in comparison with PAs with lower climatic variation. Thus, PAs located in the north of the peninsula could be acting as essential refuges for Iberian amphibian species, where they have found the ideal conditions for maintaining high levels of functional diversity. On the other hand, amphibian diversity seems to be affected by the kind of habitats present in the PA, which usually differ from those occurring outside the PA’s boundaries (see Supplementary Appendix 2 in Data Sheet 1). Protected areas dominated by open, scarcely vegetated surfaces showed lower amphibian phylogenetic diversity, possibly because desiccation and temperature are typically higher in bare soils (Kleidon et al., 2000), which may deem these areas inhabitable for many lineages. A similar negative association occurred between phylogenetic diversity and marine waters, likely because amphibians rely on freshwater for survival. By contrast, PAs with larger arable lands showed higher values of phylogenetic diversity, which was unsurprising as agricultural lands are regularly used by amphibians, some even in areas with high-intensity crop cover (Koumaris and Fahrig, 2016; Collins and Fahrig, 2017), particularly when artificial agricultural ponds are present and become alternative breeding sites for these species (Knutson et al., 2004). Finally, regarding species proportions, previous studies showed that amphibians are particularly underrepresented in PAs, both in Europe (Lobo and Araújo, 2003; Abellán and Sánchez-Fernández, 2015; see also Araújo et al., 2007) and at a global scale (Nori et al., 2015). This contrasts with our findings, as the Iberian PAs are safeguarding a large proportion of the species found in the region where they stand.

Existing evidence suggests that reptiles are underrepresented in protected areas of the Iberian Peninsula (Araújo et al., 2007; Carvalho et al., 2011, 2017), and also of Europe (Abellán and Sánchez-Fernández, 2015; Thuiller et al., 2015). Still, we found that the functional and phylogenetic diversity of reptiles in most PAs does not differ from the adjacent region, indicating that regional diversity is well captured in the current network of natural and national parks. Note, however, that previous studies have focused on evaluating the effectiveness of the overall Iberian PAs network, and not the regional representativeness of each PA individually, so our results do not necessarily contradict previous findings, as regional diversity could be simply low, and the PAs could just be mirroring such low diversity. Protected areas hosting high habitat richness tend to be associated with reptiles’ phylogenetic overrepresentation. Also, PAs presenting lower annual temperatures and higher precipitation levels showed a tendency toward an underrepresentation of functional and phylogenetic diversities, respectively. This could be related with reptiles’ ectothermy and their need for sunlight and heat to regulate body temperature and increase metabolic rates (Clarke and Gaston, 2006; Salvador, 2014), which could also explain the trend of reptile phylogenetic diversity to be underrepresented in PAs with larger shaded forest habitats. Finally, our results show that the species composition harbored by most of the PAs reflects a large proportion of the regional reptile species richness, although the results found by Araújo et al. (2007), and by Lobo and Araújo (2003) suggest that the current locations of Iberian PAs do not capture the hotspots of Iberian reptile’s species richness.

In the case of mammals, we have found that Iberian PAs are generally adequate representations of their regional functional and phylogenetic diversity. Interestingly, although Brum et al. (2017) concluded that the current PAs global network does not adequately represent the different facets of mammalian biodiversity and advocated for its expansion, these authors did not included Iberia among the high priority regions, which is in line with our findings. However, although the climatic characteristics within the PAs did not differ significantly from the regions where they are located (see Supplementary Appendix 2 in Data Sheet 1), our results showed that both diversity facets were influenced by the climatic conditions found within the PA’s boundaries. These results agree with the findings of Santos et al. (2020), who concluded that current climate strongly affects both functional and phylogenetic diversity gradients of European mammals. Particularly, we found a tendency for higher functional and phylogenetic diversity than expected from the regional pool in PAs with high annual precipitations and low annual mean temperature. Small mammals and rodents are the most common groups of this taxonomic group in Iberia, and they usually prefer temperate, moist and rainy habitats since their basic diet components (insects and herbaceous plants) proliferate throughout the year in these areas (Vickery and Rivest, 1992; Milstead et al., 2007; Hsu et al., 2012); therefore, PAs with very warm and arid climates could be harboring species with similar traits or lineages which are adapted to survive under these climatic conditions. Conversely, PAs with larger areas covered with inland wetlands and mine, dump and construction sites showed a tendency toward an underrepresentation of phylogenetic diversity, which is not unexpected because these habitats are extremely hostile for most terrestrial mammals where only a reduced number of species can survive (Ardente et al., 2016; Lawer et al., 2019). We found higher functional diversity of mammals than expected from the regional pool in PAs with larger areas covered by arable lands, which could be due to the presence of generalist species that are known to thrive in agricultural landscapes (Dotta and Verdade, 2011; Magioli et al., 2016). Also, taxonomic-based assessment of mammal diversity indicates that this group is fairly represented in Iberian PAs (Araújo et al., 2007) and in Natura 2000 Network (Rosso et al., 2018), which is much in line with our findings (more than half of PAs host at least the 75% of mammal species from their region).

Different from other groups, bird functional and phylogenetic diversities were underrepresented in a large number of the Iberian PAs, which agrees with the findings of other studies conducted at different scales. Pollock et al. (2017) concluded that both functional and phylogenetic diversity of birds are poorly represented in protected areas worldwide because the design of these areas was based on optimizing taxonomic diversity, which is an inefficient conservation solution for other diversity facets. Zupan et al. (2014) also reported incongruences between hotspots of birds’ phylogenetic diversity and European PAs, which despite their effectiveness in preserving high species numbers might hold low phylogenetic diversity due to low extinction rates and recent events of higher diversification affecting only some clades. At smaller scales, Devictor et al. (2010) found a mismatch between the location of French PAs and the distribution of functional diversity, which appears underrepresented in comparison with taxonomic diversity. In the Iberian Peninsula, birds are the most diverse tetrapod group, comprising many habitat specialist species (Cramp et al., 2004) (e.g., Ardea purpurea or Podiceps cristatus in wetlands, and Tetrax tetrax in crop steppes). Indeed, our results indicate that PAs with high habitat diversity and richness favor the presence of functional and phylogenetic diverse bird communities, respectively. Therefore, PAs highly dominated by a particular habitat type, that might in turn be missing outside the PA’s limits (see Supplementary Appendix 2 in Data Sheet 1), are expected to harbor a reduced portion of the regional functional and phylogenetic spectrum. That could also explain the slightly lower proportion of regional bird species in the PAs when compared with other taxonomic groups. It is noteworthy that PAs including larger areas of inland wetlands correspond to those with higher representation of regional phylogenetic diversity, a result that complements the findings of Brazner and MacKinnon (2020), who identified northern Europe wetlands as hotspots of bird diversity. Our results also indicate that climate conditions influence functional and phylogenetic diversity (note that climatic characteristics within the PAs did not differ from those of the regions where they are located; see Supplementary Appendix 2 in Data Sheet 1), a result that is in line with Ferger et al. (2014), who concluded that temperature and precipitation strongly affect species richness of birds, having temperature the strongest overall effect on bird richness. However, our results do not necessarily indicate that PAs are inefficiently protecting species that are important targets in conservation programs. Rather, they indicate that PAs are safeguarding species that are functional and/or phylogenetically close (e.g., habitat specialist species) or that belong to a particular guild (Duckworth and Altwegg, 2018), thus diverging from presumably generalist species that may occur outside the PAs. Therefore, if we are to achieve a high representation in the PAs of the regional functional and phylogenetic diversity, then it might be necessary to expand the current PAs in order to increase the amount of habitat heterogeneity under protection.

On the other hand, as commented above, the proportion of regional species hosted in each PA is generally not related to any of the predictors considered, i.e., habitat diversity, habitat type, climate, or the area of the PA itself. This can be indicating that other predictors different from those considered in this study could be causing such patterns. Future studies should explore the determinants of these surprising results.

In summary, while the current network of Iberian national and natural parks generally reflects regional functional, phylogenetic and taxonomic diversity of amphibians, reptiles and mammals, birds’ regional diversity is underrepresented in several PAs. These results are somewhat unexpected, considering that the definition of these PAs has mainly been done ad hoc (Araújo et al., 2007), and that previous studies in general reported a misrepresentation of species diversity in Iberian PAs. Also, we found a very limited overrepresentation of the regional pool diversity within the parks, indicating that PAs are not hotspots of the overall regional diversity, although they might be acting as refuges for particular species, providing more suitable conditions than the surrounding region and safeguarding population declines driven by human impacts. In the long-term, the persistence of these PAs could allow them to act as refugia (sensu Keppel et al., 2012, 2018), as they could serve as retreats for biotas when conditions are unfavourable in the surrounding areas, and then act as sources of recolonization once suitable conditions are restored. Future studies should focus on other types of protected areas (e.g., Natura 2000), including those established specifically to preserve birds (i.e., Special Protection Areas for Birds, ZEPAs acronym in Spain), which were not considered in this study. Such future developments should also be complemented with the identification of hotspots of functional and phylogenetic diversity, and an evaluation of whether these can be preserved by the current PAs network. Finally, considering the ongoing and future effects of climate change, future research should also focus on the importance of PAs for biodiversity preservation in the long-term. Therefore, it is necessary to identify which geographical areas are in fact acting as biodiversity refugia, and where are the regions that will host suitable climatic conditions in the future. This knowledge will allow us to evaluate whether there is a spatial match between such areas and the current network of PAs, and eventually propose the design of new PAs that will assure the maintenance of biodiversity under the current scenario of global change.

Data Availability Statement

The original contributions presented in the study, as well as datasets and R code generated, are included in the Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

SL-C, MÁR, and AMCS conceived the ideas. SL-C and AMB compiled the distribution database. SL-C and SBC compiled the trait data. RM-V and SBC built the phylogenies. SL-C, AMCS, and AMB analyzed the data. SL-C and AMCS wrote the manuscript with the help from all the authors. All authors contributed to the article and approved the submitted version.

Funding

This research was partially supported by the MINECO (Ministerio de Economía y Competitividad, Spain; project “Bioregionalization, graph theory and simulated worlds: revisiting foundational goals of biogeography with 21st century tools”- CGL2017-86926-P), by the Universidad de Alcalá (project CCG2016/BIO-018), and by the Fundação para a Ciência e Tecnologia (project PTDC/BIA-BIC/3545/2014), supported by the Norte Portugal Regional Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF).

SL-C was supported by the Universidad de Alcalá grant “Initiation to Research Activity”. RM-V was supported by the TALENTO program of the Regional Government of the Community of Madrid (2018-T2/AMB-10332). SBC was supported by National funds by FCT (Foundation for Science and Technology), through the individual scientific employment program-contract (CEECIND/01464/2017). AMCS was supported by a “Juan de la Cierva” Fellowship (IJCI-2014-19502) funded by the current Spanish ‘Ministerio de Ciencia y Innovación”.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Alberto Jiménez-Valverde and Salvador Carranza for providing feedback on the statistical analyses and Ana Paula Portela for assisting with data collection.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.634653/full#supplementary-material

References

Abellán, P., and Sánchez-Fernández, D. (2015). A gap analysis comparing the effectiveness of Natura 2000 and national protected area networks in representing European amphibians and reptiles. Biodivers. Conserv. 24, 1377–1390. doi: 10.1007/s10531-015-0862-3