Current and future geographic patterns of bird diversity dimensions of the Yucatan Peninsula and their representativeness in natural protected areas

ABSTRACT

The deviation of expected relationships between taxonomic (TD), phylogenetic (PD), and functional (FD) diversity may inform about some processes (speciation, extinction, competition, and migration) responsible for current biodiversity patterns. We studied the relationship between different dimensions of bird diversity (TD, PD, and FD) of the Yucatan Peninsula (YP) in a context of future climate change. We used ecological niche models to predict the current and future potential distribution of 257 bird species of the YP and estimate their TD, PD, and FD. We calculated a standardized effect size of PD (SES.PD) and FD (SES.FD) to provide an estimation of phylogenetic and functional diversity, independent from species richness. Finally, we evaluated the effectiveness of the system of PAs comparing the observed diversity values in each PA versus what is expected for a null model. We found a positive correlation between PD and TD, and a negative correlation between FD and TD in the current and future scenarios. Finally, we found that none of the PAs protect more diversity of birds than expected by a null model. Our results suggest that macroevolutionary processes have played an important role in the composition of the current Yucatan Peninsula avian assemblages. Our assessment of the effectiveness of the PAs suggests the need to adopt an integrative approach to biodiversity conservation in the YP.

KEYWORDS:

• Birds diversity
• climate change
• conservation
• Yucatan Peninsula

Introduction

Large-scale geographical biodiversity patterns and the exploration of the underlying mechanisms that structure them are central topics of ecology and biogeography [1]. Among the different dimensions of biodiversity, patterns of taxonomic diversity (TD) are the most recurrently studied [2]. Environmental variations and habitat heterogeneity are important factors in the spatial configuration of species richness patterns for different taxonomic groups. Places with greater complexity offer more resources, shelters, and opportunities for isolation where divergent adaptation improves the coexistence, persistence, and diversification of species [3]. Frequently, conservation efforts focus on TD as a measure in prioritizing conservation areas [4,5]. However, since TD ignores differences between species in terms of other key features for representing the diversity of the living world, such as functional and evolutionary characteristics [6], prioritize for conservation-based only on this measure fails to capture other biodiversity dimensions [7]. Hence, incorporating other dimensions of biodiversity when establishing conservation goals implies properly identifying the differences among them [8]. For instance, phylogenetic diversity (PD) incorporates the amount of evolutionary history among species [9,10], so it can be used to quantify the uniqueness of species in terms of evolutionary history [11]. Previously, it has been postulated that areas such as ecotones or biogeographic contact can present high PD values, on the one hand because they tend to present high speciation rates [12] or alternatively, these zones can host assemblages from both biogeographic regions, which implies an important role for dispersion [13]. In contrast, functional diversity (FD) denotes morphological, physiological, and ecological traits within communities and is useful to understand current dynamics of ecosystems, their resilience, and environmental services [14–16]. It is known that the FD of taxonomic groups such as birds increases with the diversity of the vegetation cover present in a place [17,18].

Although TD is expected to be positively related to PD and FD, the precise form of these relationships may vary spatially [19]. Communities with a similar TD value may differ in the phylogenetic relationships of species as a result of differences in their evolutionary histories (PD) as well as in their functional characteristics (FD) due to the existence of contrasting environmental conditions [20]. The deviation of expected relationships between TD, PD, and FD may inform about some processes (speciation, extinction, competition, and migration) responsible for current biodiversity patterns [21]. For instance, regions with high TD but low PD may indicate areas where recent adaptive radiation has occurred [22]. Conversely, areas with low TD associated with high PD could be the outcome of the presence of ancient lineages that may result from a process of speciation with low radiation, or from the extinction of other species in the same clades [23]. On the other hand, FD may be higher than expected by TD when resources are limited, presenting an overdispersion in the functional traits; in the opposite scenario, FD would be lower than expected by TD in the presence of strong environmental filters that could imply functional redundancy [23].

Adopting an integrated vision of biodiversity dimensions is also a challenge for conservation planning since global changes can have a specific effect on different diversity dimensions [24,25]. It has been found that Climate Change (CC) has altered demography, reproductive success, and the geographic distribution of many bird species [26,27]. Ecological niche models (ENMs) are increasingly being used to predict the potential effects of future CC on species geographic ranges [28]. ENMs correlate the presence of species to environmental conditions, to estimate their physiological tolerances [29]. Species could move through the landscape as the climate changes, thus modifying their ranges of distribution [30]. In this case, the integration of these temporal dynamics, together with the different dimensions of biodiversity, could maximize conservation efforts in the face of current climate change [31]. This is since through ENMs it is possible to make estimates of how species distributions may change in the future, and this can inform conservation strategies.

Understanding the spatial patterns of diversity is crucial in the study of the factors that modulate them, as well as the aspects that must be considered in their conservation, in countries with high biodiversity. Within Mexico, the Yucatan Peninsula (YP) is one of the richest areas harboring ~ 38% of the birds of the country [32]. The YP has been defined as a biogeographic province [33,34], with four distinctive biotic zones regarding the avifauna. The first is located in the northwest over the Yucatan state, and it is defined by tropical dry forest habitats. The second ranges from the center of the state of Yucatan and the northern part of the state of Quintana Roo, and the predominant habitats are deciduous and sub-deciduous tropical forests. The third comprises the state of Campeche and the south of the state of Quintana Roo. This zone is distinguished by presenting higher humidity and precipitation, where the evergreen forest predominates. And finally, the fourth corresponds to the island of Cozumel that contains elements of the second biotic zone, but also of the Antilles [34].

Here, we used ENMs to study the relationship between different dimensions of bird diversity of the Yucatan Peninsula in a context of future climate change. We selected the YP since this region will face different potential threats due to climate change, such as increases in temperature, changes in the frequency and intensity of cyclones, and sea-level rise [35]. We expect that i) TD, PD, and FD correlate positively, with the relationship between PD and FD being the highest under the assumption that closely related species could share characteristics such as functional traits [9]. ii) TD will be higher in areas of high precipitation and with topographic complexity [3,36]. iii) FD is associated with the areas with greater environmental variability and diversity of vegetation types (i.e grasslands, forests) [18]. iv) the contact zones between the biotic areas within the YP present the highest PD [37]. Finally, we evaluated the current and future effectiveness of the protected areas (PAs) system of the Yucatan Peninsula for conserving bird diversity dimensions (we expected that only TD is well represented).

Methods

Study area

We conducted this study in the YP, which is located in the southeast of Mexico. Predominant climate is warm and subhumid, however, across the region there is a gradient that goes from warm-dry to cold-wet from the northwest to the southeast [38]. The predominant vegetation of the YP includes humid and tropical dry forests, as well as flood grasslands [39]. Within the YP, there are 10 PAs that are included in IUCN categories IA and IV (IA – Strict Natural Reserve, IB – Wild Area, II – National Park, III – Monument or Natural Feature, IV – Habitat Management Area/Species), which exhibit clear conservation objectives (Figure 1 and Table S1).

Figure 1. Study area. Mexican portion of the Yucatan Peninsula and Protected Areas of categories Ia to IV according to the IUCN.

Presence data

We obtained geo-referenced presence records for bird species of the YP from the Global Biodiversity Information Facility database [40] and the VertNet portal [41] Only birds that are resident throughout the year in the YP were included in this study since we consider that they are the ones that compose the regional pool and therefore need to be prioritized in the conservation schemes of the region [42]. Each presence record was verified regarding the distribution ranges reported in online databases [43]. We also filtered in environmental space presence records of each species to reduce sampling bias and model overfitting [44]. This filtering was performed using the first two main axes of a principal component analysis (PCA) carried out using as inputs the environmental predictors described in the following section. In this step, the presence records were first filtered with the first axis where one presence record was retained for each value associated with the first PCA value. The presence records selected in the first step were filtered again, in this case, using the second PCA. We used the first two components in the filtering since together they include 80% of the global variance. Filtering was carried out using the “gridSample” function from Dismo package [45]. We used the “chekerboard1” function of ENMeval package [46] to divide the resulting databases into training and testing presence records.

Environmental data

We downloaded 15 bioclimatic variables from WorldClim v1.4 [47] at a spatial resolution of 2.5’ (~5 km²). These represent annual trends and deviations of temperature and precipitation for the period 1950–2000. We obtained future analogous variables from projections of global circulation models (GCM) for the 2050 and 2070 periods (CMIP5) [48]. Two GCM were selected (CCSM4 and MIROC5) to include variation and uncertainty between predictions [49]. In addition, for each GCM we choose two representative concentration pathways (RCP); one optimistic (RCP45) and the other pessimistic (RCP85). In RCP45, it is expected that the mean temperature of the planet increase ~1°C, while in the RCP85 it is expected that global mean temperature increase ~3°C [48].

We performed a principal component analysis (PCA) to reduce multicollinearity and dimensions of environmental variables using the PCARaster function of the ENMGadgets package [50]. The components were projected to the future variables using the same package. We retained the first five components that explained 97% of the total variance, as predictors to build niche models. This analysis was applied to the environmental layers, which were previously cropped using as a mask a polygon that represents a hypothetical area of historical accessibility of each species (area M of the BAM diagram [51]). Models calibrated on M areas are less prone to including regions inaccessible to species and, also prevent the evaluation metrics from inflating [52,53]. For each species, the M was generated by joining all WWF terrestrial ecoregions [54] in which at least one presence record intersects.

Ecological niche modeling

We built ecological niche models (ENMs) for each species using Maxent [55] in the “kuenm” package [56]. This package was designed to calibrate, evaluate, and select models in Maxent varying the algorithm settings via an automated process. We defined as arguments while calibrating models: Four regularization multipliers (1, 2, 3, and 4) and the “basic” combination of features (“l”, “lq”, “lqp”, “lqpt”, “lqpth”). The candidate models were evaluated and selected in a three-step protocol. First, the significance of each model was estimated via the partial receiver operating characteristic (ROC) technique [57] by analyzing the proportion of bootstrap replicates with area under the curve (AUC) ratios > 1. Second, among the significant models, we chose those models that presented an allowed omission rate (E) of less than 5%. Thus, we are considering those that exclude some presence records that may be errors not identified in the cleaning phase (e.g. those from sink populations). Third, based on the Akaike information criterion for small samples, we selected those models with least complexity [58]. Models that met the evaluation criteria were built and transferred to future scenarios by applying two transfer procedures: clamping and truncation. With the clamping procedure, the response in areas with novel environments is fixed to the suitability value levels presented at the periphery of the calibration region in environmental space. While, under the truncation procedure, the response is set to zero if the environments in transfer areas are more extreme than those in areas across which the models were calibrated [55].

For each future period (years 2050 and 2070) and RCP scenarios we used the medians of the medians of the replicates obtained by each GCM and transfer procedure. Finally, we projected the ENM to geographic space as environmental suitability maps. Then, in order to estimate the potential distribution per species, we reclassified suitability maps into two categories (species potentially absent/species potentially present) by applying a threshold of 5% of allowed omission in the training presences. This is, training presences were ranked according to their respective suitability value. Then, we considered as threshold the least suitability value associated with the 95% occurrences with the highest values. Thus, cells with environmental values above this threshold were defined as potential distribution.

Diversity measures

We computed presence-absence matrices (PAM) using potential distribution models with the “lets.presab” function [59]; with PAMs we estimated the number and identity of species per cell (i.e. a raster pixel with the same spatial resolution of predictors ~5 km²). The above was what we considered as taxonomic diversity (TD). To estimate PD, we used cophenetic distances of 100 trees derived from the phylogeny of birds published by Jetz et al. [60], while accounting for phylogenetic uncertainty. In the case of FD, we used the Gower distance to estimate the paired functional dissimilarity matrix between all species, based on the following traits: body mass, trophic guild, forage stratum, habit (diurnal or nocturnal species), clutch size, and egg size. These functional traits include ecological, life history, and morphological characteristics, and were downloaded from the EltonTraits database [61].

We used Rao’s quadratic entropy index (QE [62];) to measure PD and FD per cell. This index indicates the average of all dissimilarities in the dissimilarity matrix, which represents the dissimilarity between a random pair of species [63], which were estimated with the “melodic” function [64]. In order to ensure that all indices were comparable in a biological sense, we transformed PD and FD values to number equivalents. We applied the corrections PDj = 1/(1-PD) and FDj = 1/(1-FD) resulting in PD and FD in terms of the effective number of species needed to obtain the given value of the diversity index per cell [65].

Identifying spatial patterns of diversity dimensions

To analyze spatial patterns of congruence and discrepancy between TD, FD, and PD we computed correlation coefficients, where significance levels were calculated using Dutilleul’s degrees of freedom (considering the spatial autocorrelation). Then, we calculated the standardized effect size of PD (SES.PD) and FD (SES.FD) to provide an estimation of phylogenetic and functional diversity, independent from species richness [66,67]. For each metric, we generated a null distribution of 10 000 randomizations by shuffling the species names while maintaining the species richness in each cell. We categorized assemblages as phylogenetically or functionally overdispersed (O.PD/O.FD) when values were significantly higher than the randomized ones, and phylogenetically or functionally clustered (C.PD/C.FD) when observed values were significantly lower than randomized values at the 5% confidence level [68,69].

Effectiveness of Protected Areas

We estimate the cumulative proportion of total taxonomic, phylogenetic, and functional diversity included in protected areas. For each diversity dimension, we ranked all cells from the most to the least diverse. Along with this decreasing diversity gradient, we calculated the cumulative percentage of protected cells. For any given value of the diversity gradient, if the protected area system is unbiased, one expects the proportion of protected cells to match the overall proportion of protected cells among all cells of the YP (i.e.13.6%). Therefore, proportions of protected cells located over and under this value reveal the over- and underrepresentation, respectively, of cells located in protected areas [8].

We compared the estimated values of TD, significantly SES.PD or SES.FD in each PA of the YP versus what is expected for a null model. We fix the values of diversity and randomize the position of the PAs within the YP based on the protocol described by Ferro [70] and Ribeiro [71] in which it is possible to maintain the size, shape, and orientation of the PAs. We carried out 1000 randomizations and a PA was considered effective when the observed TD, SES.PD or SES.FD was higher than expected in at least 95% of the randomizations (p-value ≤ 0.05).

Results

The regional pool of the YP birds consisted of 257 species that represented 56 families and 21 orders. In the current scenario, the highest values of TD were found in the south and the eastern region of the peninsula between the states of Campeche and Quintana Roo (maximum values range between 225 and 252 species). While in the center and west coast of the YP we found the lowest number of species (less than 100 species). The highest PD was located towards the south-central part, extending eastward, also including the island of Cozumel. The areas with the highest FD corresponded to the center of the YP (Figure 2).

Figure 2. Change in time of the spatial configuration of the bird diversity of Yucatan Peninsula. Small maps represent diverse values for each future scenario.

Abbreviations: TD: Taxonomic diversity, PD: Phylogenetic diversity, FD: Functional diversity.

In future scenarios, we found that the areas with the highest TD would be restricted to the north and south of the state of Quintana Roo. The maximum PD values would be located on the basis and in the center of the YP (in the areas of lower TD), while the areas of maximum FD would be dispersed in the center, northern-western, and southwest of the YP. When comparing the spatial distribution of diversities in future scenarios regarding the current, there would be a potential increase of species up to 8% towards the coast of the YP and reductions in TD towards the center of the YP and particularly on the swamps and mangroves of the north coast of the state of Yucatan, ranging between 29% and 98%. The increase of PD is expected to be higher in the south-western and center of the YP, conversely, there would be areas of reduction of up to 41% of PD in the north and west of the state of Quintana Roo. The areas of FD potential increase would be dispersed concentrically in the YP, with an increase of up to 45% near the Sian Ka’an natural reserve. Rather, our models show no potential losses of FD (Figure 2).

We found a positive correlation between PD and TD, and a negative correlation between FD and TD under current climatic conditions (Table 1), in which about 57% and 76% of the variation of PD and FD, respectively, is not explained by TD. For the future scenarios, the dynamics of these relationships were variable. The relationship between TD and PD decreased compared to the current distributions, while the relationship between TD and FD increased. On the other hand, PD and FD in all cases showed a positive correlation, being higher in the present (Table 1). We only identified areas of functional over-dispersion in the center and north of the YP, at current. These areas could increase to the northeast and southwest of the YP in future scenarios (Figure 3).

Table 1. Pearson correlation coefficients with Dutilleul correction between taxonomic (TD), phylogenetic (PD), and functional (FD) diversities for the different scenarios

	TD-PD Correlation	TD-FD Correlation	PD-FD Correlation
Present	0.432***	−0.238*	0.482***
2050 rcp45	0.247*	−0.306***	0.367***
2050 rcp85	0.117	−0.335***	0.371***
2070 rcp45	0.054	−0.446***	0.424***
2070 rcp85	0.074	−0.384***	0.437***

p < 0.001***, p < 0.01**, p < 0.05*

Figure 3. Spatial distribution of bird phylogenetic and functional diversity across the Yucatan Peninsula. First column SES.PD, standardized effect size of PD. Second column SES.FD, standardized effect size of FD. Third column O.FD, functional over-dispersion.

When evaluating whether the cells with the highest or lowest values of the diversity dimensions were included in the PAs, we found that approximately 10% of the highest PD values are overrepresented in the PAs, while these same values for TD and FD are underrepresented. Further, it is noted that FD is always found underrepresented. In the year 2050, we forecast a scenario very similar to the current one, while by 2070 only 25% of the highest TD values would be overrepresented in less than 20% of the protected cells (Figure 4). Finally, in both periods (current and future) we found that none of the PAs protect more diversity (in any of the evaluated dimensions) of birds than expected by a null model (Table S2).

Figure 4. Cumulative proportion of taxonomic, functional, and phylogenetic diversity included in Protected Areas.

Abbreviations: TD: Taxonomic diversity, PD: Phylogenetic diversity, FD: Functional diversity.

Discussion

This study is the most complete analysis of current and future geographical patterns of bird diversity in the YP. Higher TD in the south of the YP is consistent with the peninsular effect previously observed among other taxonomic groups in this region where species richness decreases from base to tip [72]. Moreover, the high TD east of the peninsula coincides with the region with the highest rainfall, where dominant vegetation corresponds to evergreen tropical forest, thus fulfilling the spatial prediction of the distribution of TD. The structural complexity, both horizontal and vertical, of this type of vegetation, allows the coexistence of many birds, as it offers them shelter and food resources [73,74].

According to Arita [75], current fauna of the north of the YP is of recent origin. Also, it is suggested that these bird assemblages would have originated largely by the dispersal of species from the base of the peninsula. This idea supports the high PD found at the south of this region, indicating that bird communities in these areas have a high singularity in terms of evolutionary history [9]. Our results further highlight that some areas with a higher PD agree with the contact areas of the second and third biotic zones proposed for the YP [34], suggesting possible dispersal events of different groups of species between these zones [37]. In the areas with the highest FD of birds (center of the YP over the state of Yucatán) predominates sub-deciduous forests. Structural variation of these vegetation communities could facilitate the presence of different food guilds and foraging strategies and could also differentially influence the success of birds nesting, promoting the existence of species with different reproductive strategies [18,76].

The areas of TD loss coincide with the areas that will potentially have the higher increases in temperature and sea level. In these areas, the predominant habitats are mangroves and wetlands. These natural systems are considered especially vulnerable to climate change due to their greater exposure to droughts (in the case of wetlands) [73], and changes in hydrological cycles due to variations in the rate of elevation of sediments and the sea level (in the case of mangroves) [77]. On the other hand, our future projections estimate potential PD changes ranging from 8% gains to 41% losses in the YP. Although these percentage changes are lower than expected for TD, they suggest a future spatial restructuring of species distributions resulting for instance in the confluence of phylogenetically distant species on the southwestern and center of the YP coinciding with the areas of greatest TD loss. On the other hand, the potential changes in FD would be given by the functional differences or similarities of the species that are added or lost in an assembly (cell) with respect to the species for which their favorable conditions remained there [78]. For example, near the Sian Kan’an natural reserve, the maximum potential increase in FD is estimated and it is in these same regions where an increase in TD is expected, therefore, it is likely that the species that are added there possess unique functional traits or different from the species found in the current scenario.

The positive relationship between TD and PD observed here coincides with previous studies at continental and global scales [21,37]. However, the negative relationship between TD and FD differs from what has been previously reported in the literature [23,37]. This could indicate that in the areas with high taxonomic richness inhabit species with similar functional traits. This has several implications. First, it could increase competition between species, as they are so functionally similar. Second, competition could lead to small population sizes in the species, due to its negative effect on the population growth rate of the species [76]. Third, functional redundancy in communities with a high number of species could promote the functional stability of an ecosystem, in the face of change scenarios (climate change, disturbance) since, in functional terms, one species can replace another [14]. The coexistence of functionally similar species could be favored by the low topographic heterogeneity of the YP, due to the lack of physical barriers to the dispersal of individuals. Additionally, previous studies have found that the functional redundancy of some bird communities in the Neotropics is associated with low seasonality and high resource availability [76], suggesting that finer niche subdivision may be a greater driver in allowing coexistence of species with high functional similarity [77,7980–85]. For instance, the largest functional group in this study are the insectivorous birds which may be exploiting different resources (i.e. different heights, seasons). On the other hand, a variation between species has also been found in functionally similar ones, as in the case of ardeid species where they present different heights of their nesting sites, even when they present similar clutch sizes [78].

Positive values of functional overdispersion were observed in the north-central part of the YP indicating signs of high dispersion in the functional traits of the communities that occur in that region. Among the factors that could generate high FD values are competition, high environmental heterogeneity, and rapid evolution of traits [23]. At current in the northeast of the YP there is a higher FD than expected by the TD, but in future conditions, they will not be the same areas.

Our results have implications for the conservation of birds in the YP. Protected areas are the basis of current conservation strategies and therefore have a vital role in the protection of species against climate change. Even when potential losses in FD are not estimated, it is interesting that the highest values of this measure of diversity were always underrepresented in the PAs. This indicates that if the conservation purpose is centered on the highest values of diversity, the TD and PD do not correspond to a good surrogate of the FD. This suggests the need to avoid implementing strategies where only one measure of diversity is used [7,11] Although the highest TD and PD values were represented in the PAs in some scenarios, we found that they are not effective in protecting the diversity dimensions of the birds evaluated when compared to what is expected by a system of PAs established randomly in the YP. Therefore, it suggests the need for complementary analysis to evaluate the effectiveness of PAs in the protection of bird diversity. That is, on the one hand, it could be estimated if the beta diversity of the three dimensions is well represented [8] or it could also be identified which could be potential priority areas, so that the different diversities could be well represented in the present and in the future, using methodological strategies as Zonation [86].

This study represents a first approach in the evaluation of the species vulnerability and it also could help as a basis to inform decisions regarding the large infrastructure projects that will take place in this region (e.g. the Mayan train, and various solar power and wind energy stations), to minimize the negative effects of climate change on the diversity of birds present in the YP.

Author contributions

C.Y.A and F.V designed the study. J.E.L collected the data. J.E.L and S.C.Q did the data analysis. J.E.L wrote the manuscript with revisions of all authors.

Acknowledgments

This article is part of the requirements for obtaining the Master’s degree in Biological Sciences, in the field of ecology of Posgrado en Ciencias Biologicas of the Universidad Nacional Autonóma de México. We thank the Posgrado en Ciencias Biologicas UNAM for support in the academic training, also Bruno R. Ribeiro for sharing the protocol of randomization of PAs, PAPIIT UNAM IA205817 and Secretaria de Investigación, Innovación y Educación Superior del estado de Yucatán (Yucatan Initiative Project – Research Coordination Network: Biodiversity, Genomics, and Niche Modeling Across Multiple Scales in Yucatan) for funding.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Consejo Nacional de Ciencia y Tecnología [CVU: 889404]; Secretaria de Investigación, Innovación y Educación Superior del estado de Yucatán [Yucatan Initiative Project – Research Coordination Network: Biodiversity, Genomics, and Niche Modeling Across Multiple Scales in Yucatan]; Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (PAPIIT [IA205817].