Spatial patterns of genus-level phylogenetic endemism in the tree flora of Mediterranean Europe

Diversity and Distributions. 2021;27:913–928. | 913wileyonlinelibrary.com/journal/ddi Received: 1 October 2020 | Revised: 14 January 2021 | Accepted: 17 January 2021 DOI: 10.1111/ddi.13241 BIODIVERSITY RESEARCH Spatial patterns of genuslevel phylogenetic endemism in the tree flora of Mediterranean Europe Marwan Cheikh Albassatneh1,2,3 | Marcial Escudero4 | Anne- Christine Monnet1,2 | Juan Arroyo4 | Gianluigi Bacchetta5 | Francesca Bagnoli6 | Panayotis Dimopoulos7 | Arndt Hampe8 | Agathe Leriche2 | Frédéric Médail2 | Toni Nikolic9 | Loïc Ponger10 | Giovanni Giuseppe Vendramin6 | Bruno Fady11 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Diversity and Distributions published by John Wiley & Sons Ltd. 1Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum national d’Histoire naturelle (MNHN), CNRS, Sorbonne Université, EPHE, Université des Antilles, Paris, France 2Aix Marseille Université, Avignon Université, CNRS, IRD, IMBE, Aix- en- Provence, France 3Institute of Ecology and Environmental Sciences, Sorbonne University, Paris, France 4Department of Plant Biology and Ecology, University of Seville, Seville, Spain 5Department of Life and Environmental Sciences, University of Cagliari, Cagliari, Italy 6CNR, IBBR, Florence, Italy 7Department of Biology, Laboratory of Botany, University of Patras, Rio, Greece 8BIOGECO, INRAE, Bordeaux University, Cestas, France 9Department of Botany, Faculty of Science, University of Zagreb, Zagreb, Croatia 10Structure et Instabilité des Génomes, "Muséum National d'Histoire Naturelle" (MNHN), CNRS, Paris, France 11Ecology of Mediterranean Forests (URFM), INRAE, Avignon, France Correspondence Bruno Fady, INRAE, URFM, Ecology of Mediterranean Forests, Avignon, France. Email: br[email protected] Funding information Fondation de France; Labex OT- Med; Fondation pour la Recherche sur la Biodiversite Editor: Chris Burridge Abstract Aim: The Mediterranean Basin is a major hotspot of plant biodiversity, including forest trees. Over the past centuries, Mediterranean forests have been fragmented and overexploited, to which the threats of climate change are now added. Our aim is to better understand patterns and processes of tree biodiversity in the Mediterranean and to provide indicators complementing the traditional approaches to biodiversity conservation based on species counts and occurrences, using georeferenced phylogenetic diversity and endemism analyses in a spatial ecological context. Location: Mediterranean Europe. Methods: Using a dated phylogeny of the 64 Euro- Mediterranean tree genera, we calculated phylogenetic diversity for all 50 × 50 km2 grid cells spanning Mediterranean Europe (n = 643) and compared values with those obtained for genuslevel taxonomic diversity. Then, we tested the relative influence of geography, past and present climate, and soil on tree diversity (phylogenetic or taxonomic) and its geographical turnover. Geographical patterns of phylogenetic endemism were inferred using the Categorical Analyses of Neo- and Paleo- Endemism (CANAPE) methodology. Results: We showed that phylogenetic and taxonomic diversity within and among cells are correlated and influenced by soil parameters as well as current, Holocene and Late Glacial Maximum climate. Southern Spain, Cyprus and some Aegean islands contained areas of disproportionately high phylogenetic diversity and a concentration of phylogenetic paleoendemics, while phylogenetic neoendemism was high in eastern Sicily. Mixed phylogenetic endemism regions were detected in southern Spain and Portugal, in the Balkans and in Crete. Main conclusions: Our phylogenetic approach provides relevant indicators for better protecting forests of the Mediterranean, encompassing past and present evolutionary processes and factors. We consider areas that show a concentration of evolutionary history manifested by high phylogenetic endemism as high priority targets for the conservation of the European tree flora. 914 | CHEIKH ALBASSATNEH ET AL. 1 | INTRODUCTION Forests and woodlands harbour immense terrestrial and aquatic biodiversity and represent one of the most speciesrich habitat type worldwide (Gibson et al., 2011; Lindenmayer, 2009). Forests provide a wide range of critically important ecosystem services such as climate regulation, biomass production, water supply and purification, pollination, and provision of habitats for forest species (Brockerhoff et al., 2013; Decocq et al., 2016; Liang et al., 2016; Mori et al., 2017; Thompson et al., 2011). Over the last centuries, forest habitats have been destroyed at rates much higher than ever observed in human history (Gaston, 2000; Johnson et al., 2017; Turvey & Crees, 2019). As a result, a high number of species have become extinct and/ or suffered severe population declines (Mace et al., 2005), with many advancing at high speed to higher categories of threat every year, notably in the biodiversity hotspots of the world (Hoffmann et al., 2010; Le Roux et al., 2019). The forests of the Mediterranean Basin, a recognized hotspot of biodiversity, follow the same trend (Gauquelin et al., 2018; Médail & Quézel, 1997). Where and how to prioritize biodiversity conservation is a key political, societal and scientific issue for modern societies. Traditionally, biodiversity assessments have been based on species counts, valuing areas in terms of species richness and number of endemics and/or threatened species, mostly per spatial units (Kier et al., 2009; Myers et al., 2000; Norman, 2003). However, as species richness, endemism and threat can be influenced by many factors, such as the species concepts and taxonomic standards used, the spatial scale and the intensity of the sampling. Using such data can thus make biodiversity assessments and comparisons across taxonomic groups and sites limited in scope or misleading (Hillebrand et al., 2018; McKerrow et al., 2018; Van Jaarsveld et al., 1998). Limitations may also arise when considering all species equal entities and ignoring their particular functional role in the ecosystem, their associated communities and/or their evolutionary history (e.g., Doxa et al., 2020). The contribution of phylogenetic systematics to the biodiversity conservation debate stems from the very idea that species are not equivalent entities from an evolutionary point of view (Mishler, 2009). As tools to characterize species in terms of evolutionary history using DNA sequences became increasingly available, phylogenetic diversity (PD), the sum of branch lengths in a phylogeny connecting a set of taxa, was proposed as a measure to quantify the loss of evolutionary history (Faith, 1992; Nipperess et al., 2012). A major advantage of PD is that it is relatively resistant to changes in taxonomic status (not that we expect any here, as the genuslevel taxonomy of Mediterranean trees is well known and accepted, see Médail et al., 2019) and, thus, allows direct comparisons among different taxonomic groups. Additionally, PD can be also compared with other diversity measurements, such as functional diversity, including the role of evolutionary history in assembly rules in communities. For instance, facilitationmediated coexistence has been demonstrated to be driven by evolutionary history of interacting species (Valiente- Banuet & Verdú, 2007), as the functional syndromes of Mediterranean woody plants (Herrera, 1992; Verdú & Pausas, 2013) or the different niche response to fire in communities (Ojeda et al., 2010). Phylogenybased methods for measuring biodiversity have developed rapidly during the last decades, with one important area of development relating to endemism (Rosauer et al., 2009). Endemism has always been a major consideration in biogeography and one of the most important criteria when defining conservation priorities (Bacchetta et al., 2012; Linder, 2001). The concept of endemism is implicitly and closely related to the idea of irreplaceability (Margules & Pressey, 2000), and therefore of conservation value. However, the traditional definition of endemism has been narrowly interpreted as the complete restriction of a taxon to a certain area (i.e. absolute endemism). Relative endemism, which refers to the degree of restriction of taxonrange on a quantitative scale ranging from 1 (absolute endemism) to 0 (ubiquitous), broadens the concept (Crisp et al., 2001). In order to integrate evolutionary information, the concept of relative endemism has been extended to the geographic restriction of clades at any taxonomic level, a metric called phylogenetic endemism (PE, Mishler et al., 2014; Rosauer et al., 2009). Phylogenetic diversity and endemism are most often used to characterize alpha diversity (PD or PE) of given local sites. When entire geographical areas are considered, these estimators can be complemented by the quantification of phylogenetic beta diversity, that is, the turnover in PD among local sites (Graham & Fine, 2008). This estimator is similar to traditional measures of species turnover among sites with the exception that it refers to the shared and unshared branches of the phylogeny that are measured instead of the shared and unshared species. Turnover can also be measured for PE, by weighting rangerestricted branches more heavily than common branches, and is termed “phylogenetic range weighted turnover” (Laffan et al., 2016). All these measurements can enhance our understanding of biodiversity patterns and can be used for conservation assessment and planning (Thornhill et al., 2016). The Mediterranean biome worldwide contains circa 20% of the world's floristic richness on only 2% of its terrestrial surface (Médail & Quézel, 1997). Many biogeographical and phylogenetic studies have been conducted in the five Mediterranean climate regions of the world to elucidate the composition and structure of their floras (e.g., Rundel et al., 2016). Comprehensive biomewide phylogenetic surveys are also available for South Africa (Forest et al., 2007), KEYWORDS conservation, Euro- Mediterranean vascular flora, Nonmetric Multidimensional Scaling, phylogenetic diversity, phylogenetic endemism, relative environmental turnover | 915 CHEIKH ALBASSATNEH ET AL. California (Kling et al., 2019; Thornhill et al., 2017), Chile (Morlon et al., 2011; Scherson et al., 2017) and Australia (Mishler et al., 2014). Although many studies of the biodiversity of the Mediterranean Basin exist (Médail et al., 2019; Thompson, 2020), surprisingly, a comprehensive spatial phylogeny is still lacking for this region, by far the largest Mediterranean climate region in the world and the secondlargest terrestrial biodiversity hotspot of the world (Myers et al., 2000). Only spatially restricted phylogenies are available, notably in southern Spain (Molina- Venegas et al., 2015; Simon- Porcar et al., 2018). The most comprehensive phylogeny available to date for the region is the dated phylogeny at genus level for all tree taxa of the Euro- Mediterranean region of Cheikh Albassatneh et al. (2020). Here, we combine the genuslevel phylogeny of Cheikh Albassatneh et al. (2020) with the most exhaustive compilation of occurrence data for these tree genera across the European part of the Mediterranean Basin (Médail et al., 2019; Monnet et al., 2020). Our specific goal was to single out regions of high phylogenetic neo- and paleoendemism, thereby increasing our knowledge of biodiversity spatial patterns and of the ecological factors that can explain their emergence, as well as contributing to prioritizing areas of high conservation value. For this aim, we (a) calculated a series of taxonomic and phylogenetic diversity indices for each of the 50 × 50 km grid cells (n = 643) covering the entire European part of the Mediterranean Basin; (b) quantified amongcell variation and identified the main environmental factors explaining genus turnover and (c) calculated grid cell level phylogenetic endemism indicators using the Categorical Analyses of Neo- and Paleo- Endemism (CANAPE) methodology. 2 | MATERIAL AND METHODS 2.1 | Study area Our study area is the Euro- Mediterranean- region defined by Médail et al. (2019), the North Mediterranean terrestrial ecoregion in the biogeographic classification scheme of Olson et al. (2001). These biogeographical limits of the Euro- Mediterranean region are similar to those proposed by Médail and Quézel (1997), except for some mountain ridges in the Italian peninsula and the Balkans. The total area covered, including islands, is 1,610,200 km2, with a coastline approximately 45,200 km in length. The orographic variation in the area is large and elevations range from 0 metres above sea level (m.a.s.l) up to 3,482 m a.s.l. (Mulhacén, Sierra Nevada, Spain). Following the biogeographic delineation of Europe of Rivas- Martínez et al. (2004), our study area comprises fifteen biogeographical provinces (Figure 1), covering in parts or in total the following countries, from west to east: Portugal, Spain, France, Italy, Malta, Slovenia, Croatia, Montenegro, Albania, Macedonia, Greece and Cyprus. 2.2 | Data sources for tree occurrences, environmental variables and the phylogeny The woody vegetation of the Mediterranean is made of diverse types of shrublands (maquis, garrigue, phrygana) and forests, where tree morphology varies widely with environmental FIGURE 1 Distribution of the 15 biogeographical provinces (sensu Rivas- Martínez et al., 2004) of the Euro- Mediterranean region 916 | CHEIKH ALBASSATNEH ET AL. constraints (climate, geology, geomorphology, soil type), habitat types and anthropogenic activities. We adopted the criteria elaborated by Médail et al. (2019) concerning tree definition and used their checklist of 64 native tree genera. For genus occurrence and their spatial distribution, we used the data compiled by Monnet et al. (2020). We used a total of 643 grid cells of 50 × 50 km spatial resolution to completely cover the European Mediterranean Basin. The database of tree genus occurrence within each cell can be found in Appendix S1. The spatial resolution of environmental variables was originally 344 km2 for current climate and paleoclimate data (worldclim, https://world clim.org/) and 250 km2 for soil data (soilgrids, https://soilg rids.org/). To match our genus occurrence data, environmental data were aggregated to a 2,500 km2 spatial resolution (50 × 50 km) using the package raster (Hijmans et al., 2015) so as to obtain one average value per cell. For our analysis, we used 30 different environmental variables (Appendix S2) which are further described below. We used the dated genuslevel phylogenetic tree of Cheikh Albassatneh et al. (2020). This phylogeny was built using sequence data of three chloroplast DNA regions commonly used for phylogenetic and taxonomic barcoding purposes: the proteincoding rbcL and matK genes, and the noncoding intergenic spacer trnH- psbA (Hollingsworth et al., 2009; Kress & Erickson, 2007). The number of concatenated sites added up to 1,759 base pairs. The chronogram of Cheikh Albassatneh et al. (2020) is timed in million years and its topology and branch lengths are constrained using multiple calibration points. This phylogeny is largely congruent with the most recent Angiosperm Phylogeny Group classification (APGIV, Chase et al., 2016). Using a chronogram instead of a nondated phylogeny makes it possible to account for the slow mutation rate of plastid DNA and reduces the risk of not detecting the most recent phylogenetic endemism centres. 2.3 | Relative environmental turnover— diversity among cells and correlation with the environment To quantify amongcell variation and identify the main environmental factors explaining genus turnover, we estimated beta pairwise diversity using the betapart package in R (Baselga, 2012) and the following decomposition based on Sørensen's metric: 𝛽sor = 𝛽sim + 𝛽sne , where 𝛽sor is the total beta diversity, 𝛽sim is beta turnover and 𝛽sne is beta nestedness. This decomposition can be adapted to the calculation of beta taxonomic and phylogenetic diversities. We first computed incidencebased pairwise dissimilarities for the beta taxonomic (composition) diversity. The beta. pair function estimates three distance matrices: (a) turnover (replacement), (b) nestedness and (c) total dissimilarity (the sum of both components). We then computed pairwise phylogenetic dissimilarities for beta phylogenetic diversity, using the phylo.beta. pair function which similarly estimates three distance matrices: (a) phylogenetic turnover, (b) phylogenetic nestedness and (c) the sum of both values. Phylogenetic dissimilarities are based on Faith's phylogenetic diversity. For our analyses, we retained Simpson's beta turnover dissimilarity index ( 𝛽sim ) as the most relevant. Contrary to nestedness, spatial turnover implies the replacement of some taxa by others as a consequence of environmental sorting or spatial and historical constraints (Qian et al., 2005), our focus in this study. 𝛽sim values were in a matrix format (pairwise distance matrix of dissimilarities) for the 643 grid cells covering the Euro Mediterranean area. The relationships between 𝛽sim taxonomic and phylogenetic matrices were explored using Mantel tests and 999 permutations (Mantel, 1967). The diversity indices of the few cells that had a surface area less than 2,500 km2 (border cells) were calculated in the same way as those of the other cells. The spatial coordinates of a cell were those of its centroid. Relative environmental turnover (RET) was applied to examine the relationship between environmental variables and the tree genera of the Euro- Mediterranean area, using phylogenetic and taxonomic (composition) turnover (González- Orozco et al., 2016). As with previous studies (Buckley & Jetz, 2008), we used the term environmental turnover to explore rates of change of dissimilarity by taxonomic replacement in Mediterranean tree genera and their relationship to the environment depending on geographical distances. Then, we employed the βsim matrix to compute Nonmetric Multidimensional Scaling using the metaMDS method implemented in the R vegan packages (Oksanen et al., 2012). The metaMDS method performs Nonmetric Multidimensional Scaling (NMDS), a reduced representation rankorder ordination, and tries to find a stable solution using several random starts. In addition, it standardizes the scaling in the result, so that the configurations are easier to interpret, and adds taxa scores to the grid cell ordination. We reduced the number of dimensions to two axes, as recommended. We then fitted the environmental vectors onto the βsim ordination using the vector fitting envfit function of the vegan package in R (Oksanen et al., 2012). The environmental variables that best explained the patterns of turnover were then displayed as vectors only for the cases with high predictability (p < .001), assessed by 999 permutations. The “stress” values of the envfit results were used to estimate the efficiency of the NMDS ordinations. Stress values represent the difference between distances in the reduced, twoaxis dimension compared with the original complete multidimensional space. Stress values >0.05 provide an excellent reduced dimension representation of turnover, values >0.1 a great representation, values >0.2 a good representation, while values >0.3 provide a poor representation. The values of each grid cell along the first and second axes of the NMDS ordination of βsim distances were extracted and mapped into five arbitrary classes (that kept the number of cells per class as balanced as possible while emphasizing differences for easy visualization) over the study area using Quantum GIS (QGIS Development Team, 2014). NMDS values that are very different indicate strong differences in cell diversity. | 917 CHEIKH ALBASSATNEH ET AL. All environmental variables used are described in Appendix S3. All computations and analyses were performed using the R statistical environment (R Core Team, 2018). 2.4 | Phylogenetic diversity and endemism analyses Rao's phylogenetic diversity coefficient was calculated using the ade4 packages in R (Thioulouse et al., 1997). Biodiverse v.2.99 (Laffan et al., 2010) was used to calculate a set of six other diversity indices and associated randomizations for our 643 equalarea square grid cells (50 × 50 km) covering the Euro- Mediterranean area. These were as follows: Taxonomic Richness (TR), Weighted Endemism (WE; Crisp et al., 2001), Phylogenetic Diversity (PD; Faith, 1992), Phylogenetic Endemism (PE; Rosauer et al., 2009), Relative Phylogenetic Diversity (RPD; Mishler et al., 2014), Relative Phylogenetic Endemism (RPE; Mishler et al., 2014). RPE was then employed to calculate a set of four phylogenetic endemism indices using the Categorical analysis of neo- and paleoendemism (CANAPE; Mishler et al., 2014) method. Taxonomic richness (TR) is defined as the number of taxa observed in each cell of the grid. Weighted Endemism (WE; Crisp et al., 2001; Laffan & Crisp, 2003; Laffan et al., 2013) is a rangeweighted richness score, where the contribution of each taxon is weighted by its respective full geographic range in the European Mediterranean. Phylogenetic Diversity (PD; Faith, 1992) and Phylogenetic Endemism (PE; Rosauer et al., 2009) are the phylogenetic equivalents of TR and WE, respectively. PD is calculated as the sum of the branch lengths (from root to tips) of all the taxa observed in each cell whereas PE is the sum of branch lengths weighted by the full geographic range of the taxon. Branch lengths are proportional to time as the phylogeny used is timecalibrated. Relative Phylogenetic Diversity (RPD) and Relative Phylogenetic Endemism (RPE) are calculated as the ratios between the PD and PE obtained from the original tree and a reference tree with the same topology but with all branches of equal length (Mishler et al., 2014). Grid cells were assessed for their phylogenetic endemism using CANAPE (Mishler et al., 2014). CANAPE is a twostep process that assesses the contribution to PE from branches that are longer or shorter than expected, for locations that are first shown to be significantly high or low in PE. The process then assesses the significance of the RPE. All cells significant for one of these tests are classified into four nonoverlapping categories (González- Orozco et al., 2016; Mishler et al., 2014; Thornhill et al., 2016): (a) cells corresponding to centres of paleoendemism (few taxa with long branches, significantly high RPE), (b) cells corresponding to centres of neoendemism (few taxa with short branches, significantly low RPE), (c) cells corresponding to mixedendemism (rare long and rare short branches, not significant for RPE but significant denominator and significant numerator, pvalue < .05) and iv) cells corresponding to centres of superendemism (rare long and rare short branches, not significant for RPE but significant denominator and significant numerator, pvalue < .01). In all cases, endemism is meant at genus level. The significance of the observed PD, PE, RPD and RPE values was assessed using nonparametric tests based on a random reassignment of all the taxa into the grid cells. The distribution of the expected values of these indices under the null hypothesis was calculated from 999 trials of the randomization procedure. Indices in the highest 2.5% or the lowest 2.5% of the distribution were considered significant (twotailed test). For calculating the relationship between diversity indices and environmental factors, we used the factorial design of the PCA using the FactoMineR packages in R (Lê et al., 2008) and matrix of Spearman's Rank correlation coefficient using the Hmisc packages in R (Harrell & Harrell, 2019). We also calculated the relationship between the five CANAPE categories (neoendemism, paleoendemism, mixedendemism, superendemism, not significant) and the geographical and environmental variables which showed significant influence on indices of alpha diversity: longitude (Coords_X), latitude (Coords_Y), Last Glacial Maximum (LGM) maximum temperature (Tmax_LGM), Last Inter- Glacial precipitation (P_LIG), Mid- Holocene maximum temperature (Tmax_mid_H), solar radiation (SRAD) and current temperature variables (Tavg, Tmax and Tmin) (Appendix S3) overlapped with the five categories of the CANAPE results. As normality of the residuals (one of the fundamental assumptions for ANOVAs) was not fulfilled for the three comparisons (results not shown), we used a nonparametric Kruskal– Wallis test. If at least one of the distributions was significantly different from the others (p < .05), then Wilcoxon pairwise comparisons (twotailed tests with Holm's correction) were performed to disentangle which categories were significantly different from each other (p < .05). The comparison was plotted using the method of Kassambara (2017). Groups of cells showing significant PE using CANAPE were identified using an agglomerative UPGMA cluster analysis in Biodiverse v.2.99 (Laffan et al., 2010), and the rangeweighted phylogenetic turnover metric (Laffan et al., 2016). By focusing on the shared rangerestricted branches, this analysis highlights geographic regions within which the evolutionary makeup of the endemic flora is relatively homogeneous (Link- Pérez & Laffan, 2018; Thornhill et al., 2017). 3 | RESULTS 3.1 | Taxonomic and phylogenetic beta turnover— biodiversity among cells and correlation with environmental factors There was a highly significant relationship between the beta turnover taxonomic diversity and beta turnover Faith's phylogenetic diversity (Mantel test, observation = 0.8920115, simulated p value = .01). Figure 2 shows the results of the beta turnover phylogenetic diversity NMDS analysis plotted over the map of our study area. The coordinates of the NMSDS axes were aggregated into five 918 | CHEIKH ALBASSATNEH ET AL. classes (see Appendix S4 for the number of cells per NMDS class). Cells with the highest and lowest NMDS values indicate regions with high genus turnover. Particularly, Figure 2a shows that cells with high NMDS scores (light yellow cells) are located mostly in the Northern Mediterranean region, while cells with low NMDS scores (black cells) are mostly in the Southern Mediterranean and in islands. In Figure 2b, featuring NMDS axis 2, highvalue NMDS scores (yellow cells) are mostly in the Western Mediterranean region while lowvalue NMDS scores (black cells) are mostly in the Eastern Mediterranean. The scores of the two NMDS axes that characterize phylogenetic turnover could be explained with a moderately good fit (stress value of 0.2749) by different environmental variables. The NMDS plot (Appendix S5) and the table of environment vectors (Appendix S6) show the relationship between the beta turnover phylogenetic diversity and the environmental variables. Areas of either high- or low- NMDS values indicate regions with high genus turnover, whereas NMDS values close to zero indicate regions with the lowest genus turnover. The largest positive values of the NMDS axis 1 (light yellow squares in Figure 2a) were associated with Last Glacial Maximum precipitation (P_LGM), texture class (Texture_0.30), derived available soil water capacity (DASWC_30) and depth to bedrock (DBR_200), which corresponded to cells from the northern parts of the two Iberian biogeographic provinces, from the northern Cevenno- Pyrenean and Alpine provinces and the Apennino- Balkan province (Figures 1 and 2a; Appendix S5). The largest negative values of the NMDS axis 1 (black squares in Figure 2a) were associated with Mid- Holocene FIGURE 2 Phylogenetic turnover of Euro- Mediterranean tree genera measured using the nonmetric multidimensional scaling (NMDS) method of (a) axis 1 (NMDS1) and (b) axis 2 (NMDS2). Colours indicate ranges of coordinates on the NMDS axes and were kept different for axes 1 and 2 as they cover different ranges on the two different axes | 919 CHEIKH ALBASSATNEH ET AL. minimum temperature and precipitations (Tmin_mid_H and P_ mid_H), current average, minimum and maximum temperature (Tavg, Tmin and Tmax) and clay content (ClayCMFD_0.30), which corresponded to cells from the Betican and Murcian- Almerian provinces, the southern Graeco- Aegean province and most of the large islands (Figures 1 and 2a; Appendix S5). The largest positive values of the NMDS axis 2 (bright yellow squares in Figure 2b) were associated with sand content (SandCMF_0.30) and texture class (Texture_0.30), which corresponded to cells from most of the provinces of the Iberian Peninsula. The largest negative values of the NMDS axis 2 (dark purple squares in Figure 2b) were associated with Last Inter- Glacial minimum and maximum temperature (Tmin_LIG and Tmax_LIG), Last Inter- Glacial precipitation (P_LIG), water vapour pressure (WVP), cation exchange capacity of soil (CECSoil_0.30) and wind speed (WS), which corresponded to cells from most of the Central and Eastern Mediterranean provinces (Central Italy, Balkans and Cyprus) (Figures 1 and 2b; Appendix S5). 3.2 | Phylogenetic diversity and environmental drivers The Spearman's Rank correlation coefficient (Appendix S7) and factorial design of the PCA (Appendix S8) showed a positive relationship between PD, TR, and Rao phylogenetic diversity and the following variables (Appendix S3): annual precipitation (P), soil organic carbon (SoilOCST_0.15– 30), coarse fragments volume at depth 0.30 m (CFV_0.30), silt content at depth 0.30 m (SiltMF_0.30) and Last Glacial Maximum precipitation (mm/Years) (P_LGM). The extreme scattered points associated with these vectors are situated in the north and centre of Greece. On the other hand, the relationship between PD, TR and Rao phylogenetic diversity was negative with the solar radiation variables (kJ m−2 day−1; SRAD), Last Glacial Maximum temperature variables (Tmin_LGM and Tmax_LGM), Mid- Holocene temperature variables (Tmin_mid_H and Tmax_mid_H) and current maximum temperature (Tmax). The contribution of PE to the fist axes of the PCA is relatively low (Appendix S8). Relationships between PE and environmental variables appear on the 1– 4 PCA plane (Appendix S8). PE correlates positively with the Mid- Holocene precipitation (P_mid_H), and negatively with silt content at depth 0.30 m (SiltMF_0.30). All biodiversity indices were significantly correlated (Appendix S7). However, while TR was highly significantly related to PD (r = .97; Appendix S9a) without scatter (as in our beta diversity analysis), PE was less significantly related to PD (r = .45; Appendix S9b), and with high scatter, indicating that PE adds information to an analysis only focusing on TR (also see the maps of TR, PD, PE, PWE and Rao in Appendix S10). 3.3 | Categorical analysis of neo- and paleoendemism (CANAPE) Relative phylogenetic diversity divided the study area into three broad parts (Figure 3). In its western part, significantly high RPD indicating a concentration of long phylogenetic branches was found in the south of Iberian Peninsula while significantly low RPD, indicating a concentration of short phylogenetic branches, was found in the north of Iberian Peninsula (northeastern Portugal). In its central part, significantly high RPD was in the northeastern Italian Peninsula while significantly low RPD was in the islands (Sicily, Sardinia and Corsica) and the edges of the Italian Peninsula. In its eastern part, significantly high RPD was scattered in the north of the Balkan Peninsula and the islands of Crete and Cyprus. The CANAPE analysis identified 76 cells of high phylogenetic endemism (22 paleoendemism, 10 neoendemism, 40 mixedendemism and 4 superendemism cells) in Mediterranean Europe. Areas dominated by paleoendemism were in southeastern Spain as well as in the islands of Psara and Rhodes in Greece and in FIGURE 3 Map of the relative phylogenetic diversity (RPD) significance in Mediterranean Europe. Areas in blue have a concentration of significantly longer branches than expected; whereas, areas in red have a concentration of significantly shorter branches than expected 920 | CHEIKH ALBASSATNEH ET AL. Cyprus (Figure 4). The areas dominated by neoendemism were mostly located in eastern Sicily while the areas dominated by superendemism were scattered within few cells at the margins of the study area, in southern Portugal, in southeastern Italy, in southwestern Croatia and in the Greek Islands of Samothraki and the Cyclades (Naxos, Keros, Amorgòs, Santorini and Anafi). The areas dominated by mixedendemism were located at the eastern and western edges of our study area, in the southern and western regions of the Iberian Peninsula, in the Belasica Mountains of Bulgaria and Greece, and in Crete (more details in Appendix S10). Contributions from genera such as Juniperus, Chamaerops and Tetraclinis in the Iberian Peninsula, Liquidambar, Cupressus, Phoenix, Pyrus and Ceratonia in Greece, Acer, Cedrus, Cupressus, Quercus and Sorbus in Cyprus and Zelkova, Genista, Cytisus and Chamaerops in Sicily may explain at least in part the phylogenetic endemism structures observed. Several environmental variables could explain the distribution of the CANAPE categories of endemism (Figure 5; Appendices S11 and S12). Latitude (Coords_Y) and longitude (Coords_X) were correlated to paleo and mixed phylogenetic endemism, with many cells in either category located in southern islands and in either the Iberian or the Balkan Peninsula. Paleoendemism was more frequent in cells with relatively high LGM and mid- Holocene maximum temperature (Tmax_LGM and Tmax_Mid_H), while mixed and neoendemism cells were grouped in the areas of low Tmax_LGM. Current minimum temperature (Tmin) positively affected the distribution of Paleo and neoendemism. Paleoendemism was also more frequent in cells with relatively high Late Inter- Glacial precipitation (P_LIG). Paleoendemism cells were in the areas of the highest solar radiation (SRAD), in contrast to mixed and neoendemism cells which were grouped in the areas of low SRAD. The rangeweighed phylogenetic turnover analysis showed that the greatest dissimilarity was observed between the coast of the southern and western Iberian Peninsula, northeastern Sicily and Crete, and the rest of the Euro- Mediterranean region which formed a discrete cluster (Figure 6). 4 | DISCUSSION 4.1 | Taxonomic diversity, phylogenetic diversity and environmental factors Our results demonstrate a strong overall congruence between phylogenetic diversity (sensu Faith, 1992) and taxonomic richness of the genera of European Mediterranean trees. Although there is an expected congruence between phylogenetic and taxonomic diversity overall, strong spatial congruence is often rare because the shape of phylogenetic trees results from evolutionary processes that are not accounted for at taxonomic levels and because locally, taxonomic assemblages result from nonrandom ecological processes acting on regional pools (Cadotte & Tucker, 2018). Our analysis also highlighted areas of high beta turnover taxonomic diversity (and thus of high beta turnover phylogenetic diversity) within biogeographic provinces and pointed out the environmental factors that could explain their spatial structure. While the high diversity of the southernmost and northernmost cells of our study area were correlated with many relatively recent (Holocene FIGURE 4 Map of the major centres of phylogenetic endemism in Mediterranean Europe. Green grid cells do not show any signal of phylogenetic endemism. Red grid cells contain significantly lower Relative Phylogenetic Endemism (RPE) than expected given random sampling of the same number of taxa from a null tree, termed “centers of neoendemism.” Blue grid cells contain significantly higher RPE than expected, termed “centers of paleoendemism.” Purple grid cells are a mixture of neoendemism and paleoendemism, the most highly significant of which (darker purple) are termed “centers of superendemism.” Note that due to the scale of the map, not all significantly high phylogenetic endemism cells are visible | 921 CHEIKH ALBASSATNEH ET AL. FIGURE 5 Boxplot of the distribution of the environmental variables depending on CANAPE category of each of the 643 cells analysed. The six environmental variables displayed (Coords_X, Coords_Y, Tavg, Tmin, LIG_prec and SRAD) are significantly different (Kruskal– Wallis test, p < .05) among CANAPE categories (see Table S3: p < .05; pairwise comparison using Wilcoxon test with HB correction). For each box, the bold horizontal line corresponds to the median; the lower and upper bounds of the box correspond to first and third quartiles, respectively; the upper vertical line extends from the upper bound of the box to the highest value of the distribution, no further than 1.59 interquartile range (IQR, or the distance between the first and third quartiles); the lower vertical line extends from the lower bound of the box to the lowest value of the distribution, no further than 1.59 IQR; black dots are values beyond IQR (“outlier” values) -10 0 10 20 30 Mixed Neo Not_Sig Palaeo Super CANAPE X_sdrooC 36 38 40 42 44 46 Mixed NeoNot_Sig Palaeo Super CANAPE Y_sdrooC 30 40 50 60 70 MixedNeo Not_Sig Palaeo Super CANAPE GIL_P 5 10 15 Mixed NeoNot_Sig Palaeo Super CANAPE nimT 5 10 15 20 MixedNeo Not_Sig Palaeo Super CANAPE Tavg 12000 14000 16000 Mixed NeoNot_Sig Palaeo Super CANAPE DARS 928 | CHEIKH ALBASSATNEH ET AL. diversity. Molecular Ecology, 27, 1296– 1308. https://doi.org/10.1111/ mec.14520 QGIS Development Team (2014). QGIS Geographic Information System. Open Source Geospatial Foundation. Qian, H., Ricklefs, R. E., & White, P. S. (2005). Beta diversity of angiosperms in temperate floras of eastern Asia and eastern North America. Ecology Letters, 8, 15– 22. https://doi. org/10.1111/j.1461- 0248.2004.00682.x R Core Team (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Raven, P. H., & Axelrod, D. I. (1978). Origin and relationships of the California flora. University of California Press. Rivas- Martínez, S., Penas, A., & Díaz, T. (2004). Biogeographic Map of Europe. Sociedad Española de Geobotánica. Rodríguez- Sánchez, F., & Arroyo, J. (2008). Reconstructing the demise of Tethyan plants: Climatedriven range dynamics of Laurus since the Pliocene. Global Ecology and Biogeography, 17(6), 685– 695. Rosauer, D., Laffan, S. W., Crisp, M. D., Stearns, S. C., & Cook, L. G. (2009). Phylogenetic endemism: A new approach for identifying geographical concentrations of evolutionary history. Molecular Ecology, 18, 4061– 4072. https://doi.org/10.1111/j.1365- 294X.2009.04311.x Rundel, P. W., Arroyo, M. T. K., Cowling, R. M., Keeley, J. E., Lamont, B. B., & Vargas, P. (2016). Mediterranean biomes: Evolution of their vegetation, floras, and climate. Annual Review of Ecology, Evolution, and Systematics, 47, 3 8 3 – 4 0 7 . h t t p s : / / d o i . o r g / 1 0 . 1 1 4 6 / a n n u r e v - e c o l s y s - 1 2 1 4 1 5 - 0 3 2 3 3 0 Sandel, B., Arge, L., Dalsgaard, B., Davies, R. G., Gaston, K. J., Sutherland, W. J., & Svenning, J. C. (2011). The influence of Late Quaternary climatechange velocity on species endemism. Science, 334(6056), 660– 664. Scherson, R. A., Thornhill, A. H., Urbina- Casanova, R., Freyman, W. A., Pliscoff, P. A., & Mishler, B. D. (2017). Spatial phylogenetics of the vascular flora of Chile. Molecular Phylogenetics and Evolution, 112, 88– 95. Simon- Porcar, V., Escudero, M., Navarro, L., de Castro, A., Lorite, J., Molina- Venegas, R., & Arroyo, J. (2018). Using floristics, modern systematics and phylogenetics for disentangling biodiversity hotspots across scales: A Mediterranean case study. Plant Biosystems, 152, 1293– 1310. https://doi.org/10.1080/11263 504.2018.1445131 Stebbins, G. L., & Major, J. (1965). Endemism and Speciation in the California Flora. Ecological Monographs, 35, 1– 35. https://doi. org/10.2307/1942216 Takhtajan, A. (1986). Floristic regions of the world. University of California Press. Thioulouse, J., Chessel, D., Dole, S., & Olivier, J. M. (1997). ADE- 4: A multivariate analysis and graphical display software. Statistics and Computing, 7(1), 75– 83. Thompson, I., Okabe, K., Tylianakis, J., Kumar, P., Brockerhoff, E., Schellhorn, N., Parrotta, J., & Nasi, R. (2011). Forest biodiversity and the delivery of ecosystem goods and services: Translating science into policy. BioScience, 61(12), 972– 981. https://doi.org/10.1525/ bio.2011.61.12.7 Thompson, J. D. (2020). Plant evolution in the Mediterranean. Insights for conservation (2nd ed.). Oxford Univ. Press. Thornhill, A. H., Baldwin, B. G., Freyman, W. A., Nosratinia, S., Kling, M. M., Morueta- Holme, N., Madsen, T. P., Ackerly, D. D., & Mishler, B. D. (2017). Spatial phylogenetics of the native California flora. BMC Biology, 15, 1 – 1 8 . h t t p s : / / d o i . o r g / 1 0 . 1 1 8 6 / s 1 2 9 1 5 - 0 1 7 - 0 4 3 5 - x Thornhill, A. H., Mishler, B. D., Knerr, N. J., González- Orozco, C. E., Costion, C. M., Crayn, D. M., Laffan, S. W., & Miller, J. T. (2016). Continentalscale spatial phylogenetics of Australian angiosperms provides insights into ecology, evolution and conservation. Journal of Biogeography, 43, 2085– 2098. https://doi.org/10.1111/jbi.12797 Turvey, S. T., & Crees, J. J. (2019). Extinction in the Anthropocene. Current Biology, 29(19), R982– R986. Tzedakis, P. C., Lawson, I. T., Frogley, M. R., Hewitt, G. M., & Preece, R. C. (2002). Buffered tree population changes in a quaternary refugium: Evolutionary implications. Science, 297, 2044. https://doi. org/10.1126/scien ce.1073083 Valiente- Banuet, A., & Verdú, M. (2007). Facilitation can increase the phylogenetic diversity of plant communities. Ecology Letters, 10, 1029– 1036. Van Jaarsveld, A. S., Freitag, S., Chown, S. L., Muller, C., Koch, S., Hull, H., Bellamy, C., Kruger, M., Endrody- Younga, S., Mansell, M. W., & Scholtz, C. H. (1998). Biodiversity assessment and conservation strategies. Science, 279(5359), 2106– 2108. Vela, E., & Benhouhou, S. (2007). Assessment of a new hotspot for plant biodiversity in the Mediterranean basin (north Africa). Comptes Rendus Biologies, 330, 589– 605. Verdú, M., & Pausas, J. G. (2013). Syndromedriven diversification in a Mediterranean ecosystem. Evolution, 67, 1756– 1766. Vetaas, O. R., & Grytnes, J. (2002). Distribution of vascular plant species richness and endemic richness along the Himalayan elevation gradient in Nepal. Global Ecology and Biogeography, 11, 291– 301. BIOSKETCHES Marwan Cheikh Albassatneh is assistant professor in Syria, specialized in forest ecology and management and forest tree phylogenetics, currently with the Museum National d’Histoire Naturelle in Paris (France). Bruno Fady is a senior researcher at INRAE in Avignon (France) with wideranging interests in the population genetics and conservation of Mediterranean forest trees. All other authors, brought together in the Woodiv working group (https://www. fonda tionb iodiv ersite.fr/en/thefrb- inactio n/progr ams- and- proje cts/lecesab/ woodi v/), have interest in macroecology, community and population ecology, populations genetics, phylogeny and biodiversity of Mediterranean plants. Author contributions: M.C.A: methodology, data curation, data analysis, original draft preparation, writing, M.E.: methodology, data analysis, original draft preparation, writing, A.- C.M.: data curation, original draft preparation, reviewing, J.A.: reviewing, G.B.: reviewing, F.B.: reviewing, P.D.: reviewing, A.H.: reviewing, editing, A.L.: reviewing, F.M.: reviewing, T.N.: reviewing, L.P.: data analysis, G.G.V.: reviewing, B.F.: original draft preparation, writing, reviewing, editing and supervision. All authors took part in the conceptualization of the work and have read and approved the final version of the manuscript. SUPPORTING INFORMATION Additional supporting information may be found online in the Supporting Information section. How to cite this article: Cheikh Albassatneh M, Escudero M, Monnet A- C, et al. Spatial patterns of genuslevel phylogenetic endemism in the tree flora of Mediterranean Europe. Divers Distrib. 2021;27:913– 928. https://doi.org/10.1111/ddi.13241