Distribution and relative age of endemism across islands worldwide

Islands have remarkable levels of endemism and contribute greatly to global biodiversity. Establishing the age of island endemics is important to gain insights into the processes that have shaped the biodiversity patterns of island biota. We investigated the relative age of monocots across islands worldwide, using different measures of phylogenetic endemism tested against null models. We compiled a species occurrence dataset across 4,306 islands, and identified 142 sites with neo-, paleo-, mixed and super-endemism. These sites were distributed across the world, although they tended to be more common at low latitudes. The most frequent types of endemism were mixed and super-endemism, which suggests that present-day island biodiversity has frequently been shaped by processes that took place at different points in times. We also identified the environmental factors that contributed most to different types of endemism; we found that latitude, habitat availability and climate stability had a significant impact on the persistence of ancient taxa and on recent diversification events. The islands identified here are irreplaceable both for the uniqueness and the evolutionary history of their flora, and because they are a source of “option values” and evolutionary potential. Therefore, our findings will help guide biodiversity conservation on a global scale.


DATA FILTER
Step 1 Step 2 Step 3 Step 4 Step 5 Supplementary Figure S2: Selection procedure of GBIF data.

A. Corrected PE, PEalt and RPE values -addendum
A limitation in the GBIF data is the difference in geographic coverage between islands 1 . We took this uncertainty into account using a two-step procedure. First, we calculated an index of geographic coverage for each island. To do so, we used modelling to predict species richness on each island. We estimated how species number on an island varied according to a We performed calculations to correct both PEE and PER. However, the identification of areas of neo-, paleo-, mixed and super-endemism already relied on null models that removed dependency on species richness. We therefore did not correct, with our index of geographic coverage, the calculations that determined which category of endemism an area belongs to.
Doing so would have included circularity when looking for predictors of phylogenetic endemism. Species diversity theoretically increases with proximity to the continent, mainly because of colonization 4 . The distance from the mainland is one of the strongest predictors of plant species richness on islands 5 . However, despite having low species richness, remote islands may harbor high rates of endemism due to in situ speciation 6 . SLMP, which accounts for the size and shape of the coastline of the surrounding landmass, has been found to be a valuable isolation metric to explain island plant diversity on a global scale 7 . Latitude and longitude may highlight differences related to geographic position, especially along the Latitudinal Diversity Gradient 8 .
Latitude, longitude, and minimum distance from the continent were obtained from the UNEP-WCMC database 9 . SLMP came from Weigelt et al. 10 . ii.
Habitat availability: area (km²), elevation (meters), and number of ecoregions on an island. Species richness and speciation rates theoretically increase with area 4 , which is among the strongest predictors of plant species richness 5 . Maximal elevation is related to topography and environmental heterogeneity as, for instance, evident from the temperature decrease with altitude 10 11 . Elevation has also been shown to account to a large extent for insular plant species richness 5 .
The number of ecoregions per island is defined as the number of areas with distinct environmental conditions and distinct assemblages of natural communities sharing a large majority of species 12 .
Island area was obtained from the UNEP-WCMC database 9 . Elevation was extracted from Weigelt et al. 10 and the number of ecoregions per island was obtained from Olson et al. 12 .
iii. Bioclimatic factors: mean annual temperature (C°), mean annual rainfall (mm), temperature seasonality (C°), rainfall seasonality (mm), mean annual solar radiation (kJ m -2 day -1 ) and its standard deviation (sd), mean annual wind speed (m s -1 ) and its sd, mean annual water vapor pressure (kPa) and its sd, and isotherm (C°). Bioclimatic factors may select against lineages adapted to certain environmental conditions, a process called environmental filtering 13 . Temperature and rainfall were shown to be among the main predictors of plant species richness 5 14 15 . Wind speed may act on plant diversity by favoring long-distance dispersal 16 , influencing plant growth 17 , and selecting for lineages adapted to cool or harsh wind conditions. Evapotranspiration, measured from water vapor pressure, and solar radiation are key components related to energy availability. They were both shown to be predictors of plant richness on islands as well as on continents 14 15 16 .
Temperature, rainfall, wind speed, solar radiation and vapor pressure data were extracted from Worlclim version 2 19 , and missing data were completed with the  20 .
v. Sampling effort: the ICEr metric described in Supplementary Method S3.C.
By testing collinearity between variables we excluded from our analysis variables with a Spearman coefficient higher than 0.5, i.e. temperature seasonality, isotherm, mean annual solar radiation, standard deviation in solar radiation, mean annual vapor pressure and standard deviation in wind speed.

C. Estimating sampling effort
To estimate sampling effort we chose the Incidence-based Coverage Estimator (ICE 21 ) which is calculated from the number of rare species in a sample and from species accumulation curves. Compared with other estimators of species richness, ICE may best satisfy the requirements of an ideal species-richness estimator, including the capacity to reach a stable value independently of sample size, and low sensitivity to sampling patchiness and density.
We calculated ICE by defining a sub-sample within an island as a set of observations obtained at a given date, and using the R function spp.est from the 'Fossil' package 22 . We then calculated the ICEr index as the ratio of the observed number of species in an area over the expected number of species estimated with ICE.  and easy implementation of complex and/or multi-way interactions. We first looked for the best combination of learning rate, tree complexity and bag fraction in order to select for the minimum number of trees that achieved minimum prediction error 2 . Learning rate is used to shrink the contribution of each tree as it is added to the model, tree complexity is the number of nodes in a tree and the bag fraction is the proportion of randomly selected data at each step for model fitting. Ideally, low learning rates and a tree complexity that reflects the true interaction order of the response being modelled (which is almost always unknown) should be used 2 . The bag fraction should be within the 0.5 -0.75 range. To allow for tree convergence a minimum of 1,000 trees is recommended by Elith et al. 2 . To satisfy these criteria we generally chose a learning rate of 0.005, a tree complexity of 5 and a bag fraction of 0.5.

D. Correlation between variables
Boosted Regression Trees algorithms were run with the 'dismo' package in 23 .
To test for the significance and the direction of the relationship between category of endemism and our set of variables, we used generalized linear models assuming a binomial data distribution. We ran one model for each possible combination of the tested variables using the dredge function of the MuMIn R package version 1.9.13 24 . However, we had to remove missing values to run multi-model selection and thus tested these models on a subset of 2,184 islands. The most complex model included the fixed additive effects of localization, habitat availability, bioclimatic, historical and sampling efforts. We then generated a set of best models, selected on the basis of Akaike's Information Criteria (AIC; adjusted for small sample size, i.e. AICc). The lower the AICc, the better the model. The 'best model set' was defined as the model set for which cumulative AICc weight (w, i.e. a measure of relative statistical support) reached 95% of the total AICc weights. Parameter estimates were averaged across the selected models using the model averaging function (full average 24 ). This procedure enabled us to account for model selection uncertainty. All statistical analyses were performed under R version 3.4.0 3 .

F. Sensitivity analysis based on variation in geographic coverage
Specifically, we ran our analysis on data-sets where islands with an index of predicted species richness inferior to 0.25, 0.50 and 0.75 were excluded. With these data-sets, nearly all islands identified previously as priority areas of endemism were again selected as such, i.e. in these 3 data-sets only 5,6 and 7 priority islands were not selected. Moreover, changes in category of endemism was minor (Table A1). The table shows that changes are more frequent in the "super-endemism" category. This was expected because the p-value required to be classified in the "super" category is very high: as our method is relative, even little change in the pvalue caused by changes in species occurrence data may influence the selection of "superendemism" areas. However, among these islands not selected as "Super" in the sensitivity analysis, a high proportion of them were actually identified as "Mixed" which is in accordance with the fact that the "super" category is part of the "mixed" category.  Table: proportion of category change when islands were excluded depending on an index of coverage