Patterns and drivers of species richness and turnover of neo-endemic and palaeo-endemic vascular plants in a Mediterranean hotspot: the case of Crete, Greece
Journal of Biological Research-Thessaloniki volume 26, Article number: 12 (2019)
Exploring species richness and turnover patterns and their drivers can provide new insights into underlying mechanisms shaping community assembly, with significant implications for biodiversity conservation. Here, we explored diversity patterns of non-endemic, neo-endemic and palaeo-endemic vascular plants in Crete, Greece, a Mediterranean hotspot of plant richness and endemism. We evaluated the relationship between α-diversity and environmental (bioclimatic variables, topography), and anthropogenic variables by Generalized Additive Models, after accounting for spatial autocorrelation. Then, we quantified turnover using the novel concept of zeta diversity (the number of shared species by multiple sites), a framework which allows to explore the full spectrum of compositional turnover, the contribution of rare and widespread species to observed patterns and the underlying processes shaping them. Finally, we explored the abiotic and biotic effects, i.e. how well one category of species (non-endemics, palaeo-endemics, neo-endemics) predicts the patterns of the other categories, on zeta diversity by multi-site Generalized Dissimilarity Modelling.
We found a strong correlation between neo-endemic and palaeo-endemic α-diversity, with climate, topography, and human impact driving species richness. Zeta diversity analysis revealed a sharper decrease of shared palaeo-endemic species, followed by neo-endemics, and then by non-endemics with the number of sites considered to estimate compositional turnover. Perhaps, the narrow distributions of palaeo-endemics as relict species and often habitat specialists, thus persisting locally, and of neo-endemics that may have not reached yet their potential geographical range, resulted in the observed zeta diversity decline pattern. Deterministic processes controlled species turnover of rare non-endemic and neo-endemic species, while deterministic and stochastic processes contributed similarly to palaeo-endemic turnover. However, stochasticity dominates in the case of widespread species in all occasions. The environmental and anthropogenic variables were poor predictors of compositional turnover, especially of widespread species. However, the non-endemic species composition was correlated to rare palaeo-endemics and neo-endemics, highlighting the importance of biotic effects in driving turnover patterns.
It seems that centers of neo-endemism of vascular plants coincide with centers of palaeo-endemism in Crete, but species richness and species turnover are shaped by different drivers.
The understanding of the diversity patterns along spatial scales provides invaluable insights into species distribution and underlying assembly processes [1,2,3] with significant implications for biodiversity conservation . Whittaker [5, 6] proposed the partitioning of diversity into three components: α-, β-, and γ-diversity defined as the species richness at local scale, the variation in species composition, and species richness at regional scale, respectively. Patterns of species richness across spatial scales have been widely studied [7,8,9], whereas there is a growing research interest for β-diversity patterns. β-diversity reflects how communities respond to environmental gradients and changes, and climate change [4, 10,11,12,13]. Among the most often used metrics of β-diversity are pairwise (dis)similarity indices (e.g. Jaccard index) quantifying changes (or similarity) in species composition between a pair of sites [14, 15]. However, pairwise β-diversity metrics do not quantify compositional differences across more than two sites failing to fully describe turnover patterns, while they are sensitive to rare species with widespread species contributing less to turnover [16,17,18].
The novel concept of zeta (ζ) diversity defined as the number of species shared by multiple sites was proposed by Hui and McGeogh  to resolve these issues. Zeta diversity by quantifying the overlap of species distributions across multiple sites, overcomes the limitation of the pairwise comparisons of many widely used β-diversity metrics, reflecting the full spectrum of multi-site compositional turnover patterns . Therefore, it offers a framework that links diversity patterns across spatial scales, with this link representing a crucial desideratum in the biodiversity conservation . Furthermore, zeta diversity allows us to evaluate the contribution of rare, intermediate-ranging and widespread species to turnover patterns . This is an important property, since although species rarity and the associated extinction risk is commonly used for prioritization of species and conservation planning [20, 21], not only rare, but also common species, the ones “that shape the world around us”  contributing substantially to the ecosystem functioning , are calling for conservation . To sum up, zeta diversity as it is linked to all facets of diversity , although is comparatively still in its infancy, has the potential to provide in-depth insights into the turnover patterns and the underlying community assembly processes driving them [18, 26,27,28,29], the species-area relationship  and the scaling of endemism .
In the Mediterranean region, the geological and historical events created a unique biodiversity hotspot where approximately 10% of the world’s higher plants are found, with an astounding number of endemic species [31, 32]. Within the region, the island of Crete is considered a plant richness and endemism hotspot . The endemic flora of Crete is composed of relict species of a past flora (palaeo-endemics sensu Stebbins and Major ) and recently evolved species (neo-endemics sensu Stebbins and Major ) . Cretan relict flora consists of lowland species, while diversification processes occurred at mountains . Palaeo-endemics have suffered range contraction due to past climatic changes (e.g. during the Pleistocene) and are restricted to a fraction of their original distribution, often persisting in refugia (e.g. cliffs in the Mediterranean area). They are considered taxonomically isolated taxa including usually mono- or oligo-typic genera. On the other hand, neo-endemics are taxa that have evolved recently, include polytypic genera, and might have not reached yet their potential distribution due to their young age. In this context, palaeo-endemics and neo-endemics are linked to “museums” and “cradles” of biological diversity i.e. speciation centers [32, 36], and therefore understanding their diversity patterns and drivers shaping them will provide useful input for biodiversity conservation.
In the present study, we explored the patterns and drivers of species richness and turnover of vascular plant species in a Mediterranean hotspot considered as center of endemism [31, 37, 38], the island of Crete, Greece. We performed the investigation for three species categories: non-endemic, palaeo-endemic and neo-endemic species. We quantified turnover with zeta diversity, a suitable approach for comparing different species groups’ turnover within the same study area , to capture the full range of compositional turnover and the contribution of species ranging from rare to widespread of each category. First, we explore if the species of different categories differ in their elevational and geographical range. Next, we ask if centers of palaeo-endemism, neo-endemism and of non-endemic species richness coincide, and which climatic and anthropogenic factors affect the species richness of different categories. Finally, we explore the turnover and the co-occurrence patterns quantified by zeta diversity of non-endemic, palaeo-endemic and neo-endemic species, and the drivers and the underlying processes shaping the observed patterns. We expect that as palaeo-endemics are narrow-ranging species that evolved in past climatic conditions and are perhaps specialists persisting only locally will exhibit a sharper decline in multi-site compositional turnover possibly driven by deterministic processes. Contrarily, given the reported reduced turnover of island floras due to the wide distribution of non-endemic species driven by both stochastic and deterministic processes , a less steep decline is expected for the non-endemic group. Finally, we expect that neo-endemics will lie between two extremes.
Crete, the largest island of Greece, is located in the southern part of the Aegean Sea. The climate is Mediterranean with long hot and dry summers and mild winters, with mean annual temperature ranging from 9.94 to 19.13 °C and mean annual precipitation from 486.50 to 1035.32 mm (data from WorldClim ). The area is characterized by mountainous terrain with three mountain massifs: Lefka Ori (2452 m, west Crete), Psiloritis (2456 m, central Crete) and Dikti (2148 m, east Crete). Biogeographically, the island belongs to the floristic region of Kriti and Karpathos (including satellite islands = the Cretan area) which comprises 2079 species and 571 subspecies . The Cretan vascular flora includes 1647 species and subspecies (hereafter species) with approximately 10% of them being endemic to Crete  and 17.6% of them being endemic to Greece [32, 43].
We used floristic data obtained by the distribution maps (grid size: 8.25 × 8.25 km, 160 grid cells, hereafter sites) of individual plant species from Crete provided by the “Flora of the Cretan area: annotated checklist and atlas” of Turland et al. . Species were classified into three categories: non-endemics (NON-E), single-island neo-endemics (NE), and single-island palaeo-endemics (PE) using data provided in Kallimanis et al.  and any subsequent analysis was performed separately for each species category. Briefly, Kallimanis et al. classified species into different categories using published studies for specific taxa (e.g. Cellinese and Smith  for endemic Campanulaceae of Crete and Greuter  about the relict Cretan flora) and a criterion of systematic isolation, congruent with Stebbins and Major  definition of palaeo- and neo-endemism. Therefore, as palaeo-endemics were assigned isolated (reproductively and geographically) species with no close relatives and as neo-endemics all taxa at subspecific rank including vicariant species. We acknowledge that literature and systematic based formulation of endemic species categories—due to lack of phylogenetic information e.g. the public database TimeTree  included only 34 out of the 165 species [a short description of the performed analysis to generate timed phylogeny is presented in Supplementary along with the tree generated by TimeTree (Additional file 1: Figure S1)]—is certainly a limitation of our study.
For each site, we obtained topographical data (elevation, aspect, and slope) from a digital surface model produced in the framework of the Reference Data Access (RDA) Action of the EU GMES/Copernicus program (Copernicus land monitoring services 2018). The climate was quantified by the 19 bioclimatic variables of the WorldClim dataset . Specifically, we aggregated grid cell values from the high-resolution WorldClim dataset (30 arc-seconds, ~ 1 km) using zonal statistics in ArcGIS 10.3 to calculate mean variable values per site. Furthermore, to quantify human effect, we calculated the percentage of human land uses (i.e. artificial surfaces, arable land, permanent crops, pastures, and heterogeneous agricultural areas) using the CORINE Land Cover 2000 database (Copernicus land monitoring services 2018) and mean human population density per site. The calculations were performed in ArcGIS 10.3 (ArcGIS® software by ESRI).
Elevational range and species distribution
We explored if species belonging to different categories demonstrated different elevational and geographical ranges. To this end, we estimated (a) the mean elevation of each species’ occurrence, (b) the range i.e. minimum–maximum elevation of occurrence, and (c) the ratio of the number of sites each species occupies to the total number of sites. Then, we evaluated differences among categories by performing permutational one-way ANOVA with the function aovp of the R package lmPerm . In the case of significant differences, we implemented Tukey HSD post-hoc tests for pairwise comparisons. Furthermore, we estimated the C-score and the NODF index to quantify species co-occurrence and identify patterns of nestedness in species composition respectively with the bipartite package .
We estimated α-diversity i.e. species richness at site level for each species category, and we explored the relationship between species richness of different categories per two by Generalized Additive Models (GAM) with Poisson error distribution and log-link function (thin plate regression splines and 3 knots per spline). Furthermore, we evaluated the effect of bioclimatic variables, topography, and human effect on α-diversity of each species category with GAMs. To account for spatial autocorrelation we included the distance-based Moran’s eigenvectors with positive values that were estimated by the adespatial package  to the GAM. Prior to analysis, we performed Variance Inflation Factor (VIF) analysis to test for multicollinearity between bioclimatic variables, topographical variables and human effect related variables with the function vifstep of the usdm package , and setting VIF criterion < 10. The analysis indicated multicollinearity and the variables satisfying the criterion and used in the analysis were: isothermality, temperature seasonality, precipitation seasonality, precipitation of the wettest quarter of the year, precipitation of the warmest quarter of the year, aspect, slope, human population density, and the percentage of human land uses.
We estimated the zeta diversity decline and the ratio of zeta diversity decline for NON-E, NE and PE vascular plants with the package zetadiv . Briefly, zeta diversity is defined as the mean number of shared species in n number of sites ; the number of sites that is used to estimate zeta diversity and is referred as zeta order (hereafter ζi for i different number of sites). For zeta order of 1, the zeta diversity corresponds to the mean α-diversity i.e. mean species richness per site (ζ1), while for order 2 is the mean number of shared species in two sites (ζ2). Therefore, all the incidence-based pairwise β-diversity metrics can be estimated using ζ1 and ζ2. Zeta diversity for zeta orders from 3 to n is computed as the number of species in common in three to the total number of sites, with higher orders of zeta reflecting the contribution of more widespread species to compositional change and lower orders providing information on the rare species. Given that the shared species of n sites are necessarily shared in the n − 1 sites and the number of shared species declines as more sites are considered to calculate number of shared species, zeta diversity declines with zeta order. The zeta diversity decline is usually described by power-law or negative exponential function and it is informative about underlying processes driving differentiation of species composition, and the role of common and rare species in species turnover patterns . Specifically, a power-law decline, i.e. species have a unique chance to occur in a site, is indicative that niche-based processes drive species turnover, while an exponential decline i.e. species have an equal chance to occur in a site, implies that stochasticity overrules. The zeta ratio ζi/ζi+1 termed as the retention rate, quantifies the proportion of species which are retained with the addition of more sites to estimate zeta diversity. The relationship between zeta ratio and zeta order (zeta order here corresponds to the denominator i + 1) inform us about the probability of a species to be retained as more sites are considered in the computation, hence for the rarity and commonness of species . Here, following Latombe et al. , we estimated zeta decline and retention rate for each species category and estimated the parametric form of the former relationship piece-wisely to detect different patterns between widespread and rare species. Additionally, we estimated the β-diversity index N*  for comparison.
We applied the framework of Multi-Site Generalized Dissimilarity Modelling (MS-GDM)  to explore the effect of environmental (bioclimatic, topographic) and anthropogenic variables transformed with I-splines on zeta diversity (abiotic model) with the package zetadiv  and the functions provided by Latombe et al. . The importance of the predictor variable is evaluated by the maximum value of the spline, and the variation in slope across splines shows at which range of the predictor variable, the latter exerts more important effect on the differences of species composition . We quantified zeta diversity by Sørensen and Simpson version of zeta diversity. The indices are rescaled to 0–1, by dividing zeta diversity with mean α-diversity and the minimum richness of sites for the Sørensen and Simspon index, respectively. The scheme is analogous to Baselga’s partitioning of β-diversity into nestedness and turnover components , with the Sørensen version of zeta diversity reflecting nestedness component thus richness-dependent turnover and Simpson reflecting “true” turnover. Furthermore, we explored the effect of NON-E composition on PE and NE composition by adding the zeta diversity of ΝΟΝ-Ε as independent variable to the MS-GDM of NE and PE (biotic model I). Finally, we tested how the NE composition affected PE composition and vice versa (biotic model II). For each MS-GDM we estimated the variance explained as the Pearson R2 obtained by the relationship between the observed and predicted z values.
All analyses were performed in R version 3.5.2 (R Development Core Team, 2018).
Elevational range and species distribution
Palaeo-endemics and neo-endemics exhibited significant differences in mean elevation and minimum elevation of occurrence (Table 1), with non-endemics (NON-E) showing higher values, followed by neo-endemics (NE). NON-E had significantly wider distributions than palaeo-endemics (PE), and NE exhibited greater, but no significant, geographical range than PE. The lowest value of C-score and the lowest value of NODF were estimated for PE.
Among the 1647 species, 1482 were NON-E, 91 NE, and 74 PE. There were no significant differences in α-diversity between PE and NE (permutational ANOVA), but there were more NE than PE species in approximately 60% of the sites (Additional file 1: Figure S2). We observed significant positive relationships between NON-E and NE, and NE and PE (R2 = 0.24, deviance explained = 30.00%, and R2 = 0.80, deviance explained = 74.20% respectively, p < 0.05), but the relationship between NON-E and PE was hump-shaped (R2 = 0.13, deviance explained = 22.50%, p < 0.05) with declining PE species richness when NON-E species richness was high (Additional file 1: Figure S3). Furthermore, we found that all the bioclimatic, topographical, and human effect variables had a significant effect on NON-E α-diversity, after accounting for spatial autocorrelation (R2 = 0.18, deviance explained = 36.30%, p < 0.05, Additional file 1: Figure S4). NE α-diversity was significantly correlated to precipitation seasonality, slope and human population density (R2 = 0.26, deviance explained = 43.50%, p < 0.05, Additional file 1: Figure S5). Finally, temperature seasonality, precipitation of the wettest and warmest quarter and human population density had a significant effect on PE α-diversity (R2 = 0.48, deviance explained = 55.80%, p < 0.05, Additional file 1: Figure S6).
Palaeo-endemics exhibited the highest β-diversity according to N* index, followed by neo-endemic. The decline of zeta diversity (rescaled to 0–1) with zeta order was slightly sharper for PE, followed by NE (Fig. 1a). Therefore, the retention rate was lower for PE. This means that the number of shared PE species declined sharper as more sites were included for the zeta diversity estimation. The relationship of the zeta ratio with order showed that PE and NE retention rate increased sharply up to 6 and 8 zeta order respectively, and then decreased, with the decreasing part being sharper for PE (Fig. 1b). On the other hand, the retention rate of NON-E species increased with order initially and then reached a rough plateau to decrease at higher zeta orders. Based on the shape of the retention rate relationship with zeta order, we formulated two categories of commonness for each species category, differentiating species into rare (low zeta orders corresponding to the increasing part of the retention rate) and widespread (higher zeta orders corresponding to the decreasing part of the retention rate), and estimated the parametric form piece-wisely of the zeta decline for rare and widespread. The exploration revealed that the zeta decline was described by a combination of negative exponential and power-law functions for rare species, while the decline was exponential for the widespread species (Table 2). The PE exhibited similar exponential and power-law coefficient for rare species, whereas for rare NON-E and NE the power-law coefficient was greater (Table 2).
The abiotic Multi-Site Generalized Dissimilarity Model (MS-GDM) including bioclimatic variables, human population density, and the percentage of human land uses explained a small part of the variation in Sørensen and Simpson zeta diversity indices, with the explained variance decreasing with zeta order (Fig. 2a). In the case of Sørensen index, variance explained was 0.031 for NON-E, 0.013 for NE and, and 0.042 for PE for zeta order equal to one. We observed similar low values for the Simpson index (Fig. 2b). Due to the low variance explained, the abiotic MS-GDMs provided little information about environmental and anthropogenic drivers of differences in species composition (Additional file 1: Figures S7–S12).
The biotic model I MS-GDM i.e. the inclusion of zeta diversity of NON-E as predictor of differences in species turnover of NE and PE exhibited higher variance explained than the abiotic model at low zeta orders. The effect was slightly better for PE than NE, with variance explained for Sørensen index ranging between 0.001 and 0.15 (Simpson index: 0.008–0.11) for NE, and between 0.01 and 0.17 for PE (Simpson index: 0.003–0.12). The species composition of the NON-E had the highest effect on differences in NE and PE species composition. For NE Sørensen index, the biotic effect was followed in importance by topography (aspect, highest coefficient: first I-spline) and human population density (coefficient: third I-spline), while isothermality had an effect for zeta orders higher than 2 (highest coefficient: second I-spline) (Fig. 3). For PE Sørensen index, for zeta diversity of order 2 temperature seasonality (highest coefficient: third I-spline), topography (highest coefficient: third I-spline), and human population density (highest coefficient: second I-spline) contributed to differences in species composition (Fig. 3). With increasing order, we detected the effect of precipitation seasonality (highest coefficient: third I-spline) and precipitation of the warmest quarter (highest coefficient: second I-spline). Regarding Simpson index, apart from the prominent biotic effect, there was an effect of topography (aspect and slope, highest coefficient: third I-spline) and human population density (highest coefficient: third I-spline) for NE (Fig. 4). For PEs, species composition were also affected by temperature seasonality (highest coefficient: third I-spline), precipitation seasonality (highest coefficient: third I-spline), topography (aspect, highest coefficient: third and first I-spline for zeta order < 4 and zeta order = 4, respectively; slope, highest coefficient: first I-spline) and human population density (highest coefficient: second I-spline) (Fig. 4). The biotic model II for PE i.e. including NE zeta diversity as predictor of PE species turnover, exhibited higher variance explained than the abiotic model, but lower variance than the biotic model I. Biotic effect, distance and all environmental variables but temperature seasonality and human population density contributed to differences in PE Sørensen and Simpson index (Additional file 1: Figures S13, S14). However, for NE, the PE zeta diversity did not improve the performance of the MS-GDM (Additional file 1: Figures S13, S14).
Elevational range and species distribution
Cretan flora consists of a similar number of palaeo-endemic and neo-endemic species. Our results demonstrated that NE exhibited higher mean elevation and minimum elevation of occurrence than PE, in concordance with Trigas et al.  suggesting that Cretan relict flora consists mostly of lowland species, while diversification at higher (middle) elevations gave rise to neo-endemic species. Not surprisingly, endemic species were more range-restricted than NON-E, while PE exhibited narrower, but not significantly, distributions than NE. There are many possible explanations for PE and NE narrow distributions. Palaeo-endemics are relicts of a past flora occurring in a fraction of their original distribution, and niche-based processes acting over a long time may have limited them in marginal environments for which they are well-adapted. Neo-endemics may have not fully expanded their distribution due to their relatively recent differentiation. Furthermore, range-restricted endemics may be low competitors with low dispersal ability , investing in local persistence rather than in high dispersal to survive [56, 57].
There was a strong correlation between NE and PE α-diversity suggesting accumulation of endemics at similar sites across Crete i.e. centers of neo-endemism and palaeo-endemism may overlap, as it happens in the eastern part of the island. Regarding drivers of α-diversity, we found that climate, topography, and human effect were important for all categories, with the effect being stronger for PE. Local species richness of all species categories declined with human population density and showed a hump-shaped relationship with slope. Araújo  reported that human population density is positively correlated to plant species richness, but not related to species richness of narrow endemics in Europe. However, other studies found negative relationship between species richness and human population density [59, 60], in accordance with the observed pattern here. Kougioumoutzis and Tiniakou  found that human population density affected negatively the total number of endemic species of Cycladic islands in the Aegean area. According to Lavergne et al. , range-restricted endemic species in Mediterranean Basin are related to low human population density—often limited to inaccessible areas probably due to human pressure —and slope. Furthermore, topographic relief is considered to promote endemic species richness , favoring neo-endemism through spatial divergence, but also palaeo-endemism in Mediterranean hotspots [57, 64]. Regarding climatic variables, we found that precipitation (precipitation seasonality for NE, and precipitation of wettest and warmest quarter for PE) shaped species richness patterns. The role of climate as determinant of species richness patterns has been extensively studied [65,66,67,68], especially in the face of climate change [69,70,71]. Molina-Venegas et al.  argue that palaeo-endemics tend to occur in wetter conditions than neo-endemics that are related to less benign environmental conditions.
The species richness patterns and the geographical range size of different species categories were well reflected in their retention rate and the pattern of zeta diversity decline with zeta order. Specifically, PE with lowest species richness and narrower geographical range showed also a sharper decline of zeta diversity than NE and NON-E, and higher β-diversity as it was estimated by the N* index. Therefore, as zeta order increased, less PE than NE or NON-E species were shared between sites. This was more prominent in the low orders of zeta, reflecting the contribution of the rare species to the observed pattern. PE are ancient isolated species with small geographical range due to environmental change and habitat loss. Perhaps, the environmental change has rendered their traits less adaptive to the present-day environmental conditions , with PE being often habitat specialists restricted to specific environmental conditions. These species may have persisted in refugia  e.g. ecosystems at high elevations . Mountainous topography favors high endemism  e.g. cliffs in the Mediterranean are considered refugia for plant species with unique species composition . On the other hand, NE exhibited significantly lower species richness and narrower geographical range than NON-E. A possible explanation could be that NE may have not dispersed very far reaching their full potential in terms of geographical range  or they are of lower competitive ability . Finally, the zeta diversity decline for rare species was described by both power-law and negative exponential model, indicating that deterministic processes and stochasticity drive species turnover of rare species. The deterministic processes were stronger for NON-E and NE (greater power-law coefficient), while stochasticity and deterministic processes contributed similarly to rare PE turnover. The latter finding may be interpreted by the role of environmental stochasticity in shaping PE turnover (see below). Contrarily, the turnover pattern of widespread species was driven by stochasticity, as also noted by Latombe et al.  who found that species turnover of rare vascular plants with different residence time is driven by both deterministic and stochastic processes, whereas stochasticity prevails in the case of widespread species.
The environmental and anthropogenic variables were strong determinants of species richness i.e. for zeta diversity of order 1 (ζ1), but the variance explained of the abiotic MS-GDM predicting species turnover was low across zeta orders, with slightly better performance for PE. The declining explanatory power of predictors from ζ1 (species richness) to higher orders indicates that different drivers shape species richness and species turnover of vascular plants in Crete. Perhaps, environmental variables not considered here are better predictors of species turnover e.g. soil properties are considered strong drivers of β-diversity patterns , especially for the heterogeneous Mediterranean ecosystems .
Biotic effects in the Biotic model I (zeta diversity of NON-E as predictor of NE and PE turnover) increased the variance explained by the MS-GDM, especially for rare species. The variance explained increased more in the case of the richness-dependent Sørensen index. This was also indicated from the significant correlation of species richness (ζ1) of NON-E with both NE and PE species richness. The incorporation of biotic effects has shown that improve significantly the predictive power of species distribution models of plants across spatial scales [78,79,80], and our results suggest that this may apply to turnover also. These biotic effects can be negative plant–plant interactions e.g. competition affecting species distribution and composition [81,82,83] or positive e.g. facilitation of expansion of species distribution toward higher elevations . Although increased values of C-score provided evidence for disaggregation, the impact of biotic effects on turnover remains an open question and further research on the co-occurrence patterns is required. However, the predictive power of MS-GDM diminished for widespread species, suggesting that turnover of widespread NE and PE is independent of the species composition of NON-E. The environmental and anthropogenic variables exerting significant effect on species richness appeared to have a role on the richness-dependent turnover i.e. topography and human population for NE, and additionally, temperature seasonality and precipitation related variables for PE. Similar patterns were observed in the case of Simpson index. Human population density was more important when it was intermediate for NE or high for PE (i.e. urbanization effect). Furthermore, seasonality was significant for PE when it was high. The link between turnover and seasonality in Mediterranean landscapes has been highlighted for vertebrates by Martin and Ferrer , and perhaps this applies also to plants. It is highly probable that many palaeo-endemics may have emerged before the establishment of typical Mediterranean climate  and their distribution is now restricted at high elevations where the climate resembles the pre-Mediterranean one with low (or no) seasonality. On the other hand, biotic effects in model II (zeta diversity of NE as predictor of PE and vice versa) increased the variance explained by the MS-GDM only for PE, and once again especially for rare species. It seems that only species distribution of rare PE relates to the species distribution of NE. A possible explanation could be that hotspots of NE are in areas where the habitats are suitable also for PE. A potential limitation of our study is the classification of endemic species into palaeo-endemics and neo-endemics according to published studies and the systematic isolation criterion. Recent studies have attempted to classify neo-endemics and palaeo-endemics using time-calibrated phylogenetic trees. This approach might lead to differences in the observed patterns, thus future research incorporating phylogenetic information could shed more light on the diversity patterns of endemism.
The spatial patterns of neo-endemic and palaeo-endemic richness imply that hotspots of neo-endemism coincide with palaeo-endemism hotspots, perhaps due to the effect of climate, topography, and anthropogenic variables on α-diversity. More interesting were the species turnover patterns which highlighted the contribution of rare species to the observed patterns. The profile of zeta diversity decline with zeta order revealed that palaeo-endemics showed a sharper decline, followed by neo-endemics and then by non-endemics. A possible explanation could be that palaeo-endemics, as relict species less adapted to current environmental conditions and often habitat specialists, exhibited narrow distributions, and thus the number of shared species across sites decreased sharply. The narrow distribution of neo-endemics is restricted due to recent differentiation. Despite the significant effect of environmental and anthropogenic variables on α-diversity, their effect on zeta diversity was weak implying that different drivers affect the spatial patterns of species richness and species turnover of vascular plants in Crete. Finally, our analysis emphasized the role of biotic effects in driving turnover patterns.
Availability of data and materials
The biodiversity dataset used/analyzed during the current study is published in Turland et al. , the climate and the land use data were derived from WorldClim dataset  and CORINE Land Cover 2000 database (Copernicus land monitoring services 2018) respectively. All sources are publicly available.
Borcard D, Legendre P, Avois-Jacquet C, Tuomisto H. Dissecting the spatial structure of ecological data at multiple scales. Ecology. 2004;85:1826–32.
Holyoak M, Leibold MA, Holt RD. Metacommunities: spatial dynamics and ecological communities. Chicago: University of Chicago Press; 2005.
McKnight MW, White PS, McDonald RI, Lamoreux JF, Sechrest W, Ridgely RS, et al. Putting beta-diversity on the map: broad-scale congruence and coincidence in the extremes. PLoS Biol. 2007;5:2424–32.
Socolar JB, Gilroy JJ, Kunin WE, Edwards DP. How should beta-diversity inform biodiversity conservation? Trends Ecol Evol. 2016;31:67–80.
Whittaker RH. Vegetation of the Siskiyou Mountains, Oregon and California. Ecol Monogr. 1960;30:280–338.
Whittaker RH. Evolution and measurement of species diversity. Taxon. 1972;21:213–51.
Rahbek C. The role of spatial scale and the perception of large-scale species-richness patterns. Ecol Lett. 2005;8:224–39.
Kessler M, Kluge J, Hemp A, Ohlemuller R. A global comparative analysis of elevational species richness patterns of ferns. Glob Ecol Biogeogr. 2011;20:868–80.
Stein A, Gerstner K, Kreft H. Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecol Lett. 2014;17:866–80.
Anderson MJ, Crist TO, Chase JM, Vellend M, Inouye BD, Freestone AL, et al. Navigating the multiple meanings of beta diversity: a roadmap for the practicing ecologist. Ecol Lett. 2011;14:19–28.
Lazarina M, Sgardelis SP, Tscheulin T, Devalez J, Mizerakis V, Kallimanis AS, et al. The effect of fire history in shaping diversity patterns of flower-visiting insects in post-fire Mediterranean pine forests. Biodivers Conserv. 2017;26:115–31.
Lazarina M, Sgardelis SP, Tscheulin T, Kallimanis AS, Devalez J, Petanidou T. Bee response to fire regimes in Mediterranean pine forests: the role of nesting preference, trophic specialization, and body size. Basic Appl Ecol. 2016;17:308–20.
Hodapp D, Borer ET, Harpole WS, Lind EM, Seabloom EW, Adler PB, et al. Spatial heterogeneity in species composition constrains plant community responses to herbivory and fertilisation. Ecol Lett. 2018;21:1364–71.
Koleff P, Gaston KJ, Lennon JJ. Measuring beta diversity for presence-absence data. J Anim Ecol. 2003;72:367–82.
Chao A, Chazdon RL, Colwell RK, Shen TJ. A new statistical approach for assessing similarity of species composition with incidence and abundance data. Ecol Lett. 2005;8:148–59.
Chao A, Jost L, Chiang SC, Jiang YH, Chazdon RL. A two-stage probabilistic approach to multiple-community similarity indices. Biometrics. 2008;64:1178–86.
Latombe G, Hui C, McGeoch MA. Multi-site generalised dissimilarity modelling: using zeta diversity to differentiate drivers of turnover in rare and widespread species. Methods Ecol Evol. 2017;8:431–42.
Hui C, Vermeulen W, Durrheim G. Quantifying multiple-site compositional turnover in an Afrotemperate forest, using zeta diversity. For Ecosyst. 2018;5:15.
Hui C, McGeoch MA. Zeta diversity as a concept and metric that unifies incidence-based biodiversity patterns. Am Nat. 2014;184:684–94.
Brooks TM, Mittermeier RA, da Fonseca GAB, Gerlach J, Hoffmann M, Lamoreux JF, et al. Global biodiversity conservation priorities. Science. 2006;313:58–61.
Gauthier P, Debussche M, Thompson JD. Regional priority setting for rare species based on a method combining three criteria. Biol Conserv. 2010;143:1501–9.
Gaston KJ, Fuller RA. Biodiversity and extinction: losing the common and the widespread. Progr Phys Geogr. 2007;31:213–25.
Baker DJ, Garnett ST, O’Connor J, Ehmke G, Clarke RH, Woinarski JCZ, et al. Conserving the abundance of nonthreatened species. Conserv Biol. 2019;33:319–28.
Ellison AM. Foundation species, non-trophic interactions, and the value of being common. Iscience. 2019;13:254–68.
Mcgeoch MA, Latombe G, Andrew NR, Nakagawa S, Nipperess DA, Roigé M, et al. Measuring continuous compositional change using decline and decay in zeta diversity. Ecology. 2019. https://doi.org/10.1002/ecy.2832.
Britton AW, Day JJ, Doble CJ, Ngatunga BP, Kemp KM, Carbone C, et al. Terrestrial-focused protected areas are effective for conservation of freshwater fish diversity in Lake Tanganyika. Biol Conserv. 2017;212:120–9.
Latombe G, Richardson DM, Pysek P, Kucera T, Hui C. Drivers of species turnover vary with species commonness for native and alien plants with different residence times. Ecology. 2018;99:2763–75.
Leihy RI, Duffy GA, Chown SL. Species richness and turnover among indigenous and introduced plants and insects of the Southern Ocean Islands. Ecosphere. 2018;9:e02358.
Simons AL, Mazor R, Stein ED, Nuzhdin S. Using alpha, beta, and zeta diversity in describing the health of stream-based benthic macroinvertebrate communities. Ecol Appl. 2019;29:e01896.
Kunin WE, Harte J, He F, Hui C, Jobe RT, Ostling A, et al. Upscaling biodiversity: estimating the species-area relationship from small samples. Ecol Monogr. 2018;88:170–87.
Medail F, Quezel P. Biodiversity hotspots in the Mediterranean basin: setting global conservation priorities. Conserv Biol. 1999;13:1510–3.
Medail F. The specific vulnerability of plant biodiversity and vegetation on Mediterranean islands in the face of global change. Reg Environ Change. 2017;17:1775–90.
Stebbins GL, Major J. Endemism and speciation in the California flora. Ecol Monogr. 1965;35:1–35.
Georghiou K, Delipetrou P. Patterns and traits of the endemic plants of Greece. Bot J Linn Soc. 2010;162:130–53.
Trigas P, Panitsa M, Tsiftsis S. Elevational gradient of vascular plant species richness and endemism in Crete—the effect of post-isolation mountain uplift on a continental island system. PLoS ONE. 2013;8:e59425.
Lopez-Pujol J, Zhang FM, Sun HQ, Ying TS, Ge S. Centres of plant endemism in China: places for survival or for speciation? J Biogeogr. 2011;38:1267–80.
Myers N. Threatened biotas: “ hot spots” in tropical forests. Environmentalist. 1988;8:187–208.
Canadas EM, Fenu G, Penas J, Lorite J, Mattana E, Bacchetta G. Hotspots within hotspots: endemic plant richness, environmental drivers, and implications for conservation. Biol Conserv. 2014;170:282–91.
König C, Weigelt P, Kreft H. Dissecting global turnover in vascular plants. Glob Ecol Biogeogr. 2017;26:228–42.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Dimopoulos P, Raus T, Bergmeier E, Constanundis T, Iatrou G, Kokkini S, et al. Vascular plants of Greece: an annotated checklist. Supplement. Willdenowia. 2016;46:301–47.
Turland NJ, Chilton L, Press JR. Flora of the Cretan area: annotated checklist and atlas. London: HMSO; 1993.
Kallimanis AS, Bergmeier E, Panitsa M, Georghiou K, Delipetrou P, Dimopoulos P. Biogeographical determinants for total and endemic species richness in a continental archipelago. Biodivers Conserv. 2010;19:1225–35.
Kallimanis AS, Panitsa M, Bergmeier E, Dimopoulos P. Examining the relationship between total species richness and single island palaeo- and neo-endemics. Acta Oecol. 2011;37:65–70.
Cellinese N, Smith SA, Edwards EJ, Kim ST, Haberle RC, Avramakis M, et al. Historical biogeography of the endemic Campanulaceae of Crete. J Biogeogr. 2009;36:1253–69.
Greuter W. The relict element of the flora of Crete and its evolutionary significance. Taxonomy, phytogeography and evolution. London: Academic Press; 1971.
Kumar S, Stecher G, Suleski M, Hedges SB. TimeTree: a resource for timelines, TimeTrees, and divergence times. Mol Biol Evol. 2017;34:1812–9.
Wheeler B, Torchiano M. lmPerm: permutation tests for linear models. R package version. 2010;1(1.2). https://CRAN.R-project.org/package=lmPerm. Accessed 18 Jan 2019.
Dormann CF, Gruber B, Fründ J. Introducing the bipartite package: analysing ecological networks. R News. 2008;8:8–11.
Dray S, Blanchet G, Borcard D, Clappe S, Guenard G, Jombart T. Adespatial: multivariate multiscale spatial analysis. R package ver. 0.1–1. 2018. https://cran.r-project.org/package=adespatial. Accessed 12 Nov 2018.
Naimi B. usdm: uncertainty analysis for species distribution models. R package version. 2015;1:1–12. http://CRAN.R-project.org/package=usdm. Accessed 2 Feb 2019.
Latombe G, McGeoch M, Nipperess D, Hui C. Zetadiv: functions to compute compositional turnover using zeta diversity. R package version 1.0. 2017. https://cran.r-project.org/package=zetadiv. Accessed 10 Sept 2018.
Lazarina M, Sgardeli V, Kallimanis AS, Sgardelis SP. An effort-based index of beta diversity. Methods Ecol Evol. 2013;4:217–25.
Baselga A. Partitioning the turnover and nestedness components of beta diversity. Glob Ecol Biogeogr. 2010;19:134–43.
Lesica P, Yurkewycz R, Crone EE. Rare plants are common where you find them. Am J Bot. 2006;93:454–9.
Lavergne S, Thompson JD, Garnier E, Debussche M. The biology and ecology of narrow endemic and widespread plants: a comparative study of trait variation in 20 congeneric pairs. Oikos. 2004;107:505–18.
Molina-Venegas R, Roquet C. Directional biases in phylogenetic structure quantification: a Mediterranean case study. Ecography. 2014;37:572–80.
Araujo MB. The coincidence of people and biodiversity in Europe. Glob Ecol Biogeogr. 2003;12:5–12.
Pautasso M. Scale dependence of the correlation between human population presence and vertebrate and plant species richness. Ecol Lett. 2007;10:16–24.
Carboni M, Thuiller W, Izzi F, Acosta A. Disentangling the relative effects of environmental versus human factors on the abundance of native and alien plant species in Mediterranean sandy shores. Divers Distrib. 2010;16:537–46.
Kougioumoutzis K, Tiniakou A. Ecological factors driving plant species diversity in the South Aegean Volcanic Arc and other central Aegean islands. Plant Ecol Divers. 2015;8:173–86.
Lavergne S, Thuiller W, Molina J, Debussche M. Environmental and human factors influencing rare plant local occurrence, extinction and persistence: a 115-year study in the Mediterranean region. J Biogeogr. 2005;32:799–811.
Vetaas OR, Grytnes JA. Distribution of vascular plant species richness and endemic richness along the Himalayan elevation gradient in Nepal. Glob Ecol Biogeogr. 2002;11:291–301.
Molina-Venegas R, Aparicio A, Lavergne S, Arroyo J. Climatic and topographical correlates of plant palaeo- and neoendemism in a Mediterranean biodiversity hotspot. Ann Bot. 2017;119:229–38.
Hawkins BA, Field R, Cornell HV, Currie DJ, Guegan JF, Kaufman DM, et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology. 2003;84:3105–17.
Currie DJ, Mittelbach GG, Cornell HV, Field R, Guegan JF, Hawkins BA, et al. Predictions and tests of climate-based hypotheses of broad-scale variation in taxonomic richness. Ecol Lett. 2004;7:1121–34.
Field R, Hawkins BA, Cornell HV, Currie DJ, Diniz-Filho JAF, Guegan JF, et al. Spatial species-richness gradients across scales: a meta-analysis. J Biogeogr. 2009;36:132–47.
Petanidou T, Kallimanis AS, Lazarina M, Tscheulin T, Devalez J, Stefanaki A, et al. Climate drives plant-pollinator interactions even along small-scale climate gradients: the case of the Aegean. Plant Biol. 2018;20:176–83.
Midgley GF, Hannah L, Millar D, Rutherford MC, Powrie LW. Assessing the vulnerability of species richness to anthropogenic climate change in a biodiversity hotspot. Glob Ecol Biogeogr. 2002;11:445–51.
Thuiller W, Lavorel S, Araujo MB, Sykes MT, Prentice IC. Climate change threats to plant diversity in Europe. Proc Natl Acad Sci USA. 2005;102:8245–50.
Steinbauer MJ, Grytnes JA, Jurasinski G, Kulonen A, Lenoir J, Pauli H, et al. Accelerated increase in plant species richness on mountain summits is linked to warming. Nature. 2018;556:231–4.
Zliobaité I, Fortelius M, Stenseth NC. Reconciling taxon senescence with the Red Queen’s hypothesis. Nature. 2017;552:92–5.
Harrison S, Noss R. Endemism hotspots are linked to stable climatic refugia. Ann Bot. 2017;119:207–14.
Vogiatzakis I. Mediterranean mountain environments. Chichester: John Wiley & Sons; 2012.
Stebbins GL. The genetic approach to problems of rare and endemic species. Madroño. 1942;6:241–58.
Fitzpatrick MC, Sanders NJ, Normand S, Svenning JC, Ferrier S, Gove AD, et al. Environmental and historical imprints on beta diversity: insights from variation in rates of species turnover along gradients. Proc R Soc B. 2013;280:20131201.
Molina-Venegas R, Aparicio A, Pina FJ, Valdés B, Arroyo J. Disentangling environmental correlates of vascular plant biodiversity in a Mediterranean hotspot. Ecol Evol. 2013;3:3879–94.
Lortie CJ, Brooker RW, Choler P, Kikvidze Z, Michalet R, Pugnaire FI, et al. Rethinking plant community theory. Oikos. 2004;107:433–8.
Thuiller W, Munkemuller T, Lavergne S, Mouillot D, Mouquet N, Schiffers K, et al. A road map for integrating eco-evolutionary processes into biodiversity models. Ecol Lett. 2013;16:94–105.
Wisz MS, Pottier J, Kissling WD, Pellissier L, Lenoir J, Damgaard CF, et al. The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling. Biol Rev Camb Philos Soc. 2013;88:15–30.
Chesson P. General theory of competitive coexistence in spatially-varying environments. Theor Popul Biol. 2000;58:211–37.
Sexton JP, McIntyre PJ, Angert AL, Rice KJ. Evolution and ecology of species range limits. Annu Rev Ecol Evol Syst. 2009;40:415–36.
Aschehoug ET, Brooker R, Atwater DZ, Maron JL, Callaway RM. The mechanisms and consequences of interspecific competition among plants. Annu Rev Ecol Evol Syst. 2016;47:263–81.
le Roux PC, Virtanen R, Heikkinen RK, Luoto M. Biotic interactions affect the elevational ranges of high-latitude plant species. Ecography. 2012;35:1048–56.
Martin B, Ferrer M. Temporally variable environments maintain more beta-diversity in Mediterranean landscapes. Acta Oecol. 2015;68:1–10.
This research is co-financed by Greece and the European Union (European Social Fund- ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project “Reinforcement of Postdoctoral Researchers” (MIS-5001552), implemented by the State Scholarships Foundation (ΙΚΥ).
This research is co-financed by Greece and the European Union (European Social Fund- ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project “Reinforcement of Postdoctoral Researchers” (MIS-5001552), implemented by the State Scholarships Foundation (ΙΚΥ).
Ethics approval and consent to participate
The present study was based on analysis on already collected and published data, thus no ethical approval or consent was required.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
(i) the generated timed phylogeny analysis of the vascular plants of Crete used in our study, Figures S2–S6 (ii) the species richness patterns and their drivers of non-endemic, neo-endemic and palaeo-endemic species, and Figures S7–S14 (iii) the I-splines of abiotic (predictors: bioclimatic variables, human population density and the percentage of human land uses) and biotic II (predictors: predictors of abiotic model and neo-endemic zeta diversity as predictor of palaeo-endemic zeta diversity and vice versa) Multi-Site Generalized Dissimilarity Models showing the contribution of different predictors to explaining zeta diversity of different species categories
About this article
Cite this article
Lazarina, M., Kallimanis, A.S., Dimopoulos, P. et al. Patterns and drivers of species richness and turnover of neo-endemic and palaeo-endemic vascular plants in a Mediterranean hotspot: the case of Crete, Greece. J of Biol Res-Thessaloniki 26, 12 (2019). https://doi.org/10.1186/s40709-019-0106-x