What makes an island a suitable refuge for an endemic lizard?

To identify and evaluate the environmental properties of the Mediterranean islands that make them a suitable habitat for their endemic lizards. This would represent an important improvement in our knowledge and would assist the conservation of these lizards, which are threatened by both the expansion of mainland species and the increasing intensity of human disturbance on these islands.

Extinctions are also common on larger islands, but in some situations, endemic populations may remain depending on their interactions with the introduced mainland species or the existence of refuge habitats (Fox & Fox, 2000;Olson et al., 2006).
The Mediterranean Basin provides an interesting setting that can be used to evaluate patterns of coexistence among mainland and endemic island faunas. This basin comprises more than 3,000 islands that show great geophysical disparity and are spread over an area of 2,500,000 km 2 (Coll et al., 2010). The groups of islands include some that can be classified as "near-oceanic," such as the Balearic Islands or Crete, which have experienced long periods of isolation and endemic evolutionary radiation initiated after the Messinian salinity crisis Lalueza-Fox et al., 2005). However, most are continental islands, which remained connected to the mainland during glacial-eustatic regressions (Palombo, 2009). The isolation of all of these archipelagos ended with the first human settlers, who arrived during the Upper Epipalaeolithic (Cherry & Leppard, 2018). The arrival of humans was synchronous with significant changes in the fauna and botanical composition across the Mediterranean islands, which suggests that human activity led to significant environmental impacts (Gippoliti & Amori, 2006;Jalut et al., 2009). The changes involved the introduction of numerous continental species, a process that has been maintained until the present (Barun et al., 2011).
Despite the long history of disturbance, reptile extinctions in the Mediterranean islands have been scarce, at least when compared with other groups of vertebrates (Bover et al., 2008). These extinctions affected large species, such as giant tortoises, turtles and lizards on Malta, and isolated populations of extant species, such as Testudo marginata on Crete and Podarcis lilfordi on Mallorca and Menorca (Bachmayer et al., 1976;Böhme & Zammit-Maempel, 1982;Terrasa et al., 2009). The extinction of reptiles on the larger islands occurred in prehistoric times, possibly as a result of overharvesting by the first human settlers and the introduction of saurophagic snakes (Kotsakis, 1977;Mencía et al., 2017). On the smaller islands, extinctions followed the accidental translocation of congeneric lizards, such as P. siculus during the 20th century (Corti et al., 2014;D'Amico et al., 2018).
In this study, I evaluate the factors that influence the occurrence of endemic and mainland lizards on the Mediterranean islands. I expected greater habitat heterogeneity (e.g. island size, topography and vegetation) to be correlated with higher species richness within both groups of lizards. However, larger islands can also harbour more predators, which shifts the species balance in favour of the mainland lizards (Van Moorleghem et al., 2020). I also expected that distance from the continent would be correlated with a greater presence of endemic species (Kougioumoutzis & Tiniakou, 2015). Marine channels, even if they are relatively shallow, can constitute an almost impassable barrier for some small lizards (Hurston et al., 2009), although there are also opposite examples (Phillips et al., 2019;Raia et al., 2017). Consequently, I hypothesized that isolation relative to other islands in an archipelago would be negatively associated with both groups of lizards.
A total of 31 species of mainland lizards also occur in the same archipelagos, most of which are confined to the larger islands (Kotsakiozi et al., 2018;Speybroeck et al., 2016) (Figure 2). Data related to the presence or absence of lizards on 552 islands were obtained from biogeographic atlases (Delaugerre & Cheylan, 1992;Mayol, 1997;Sindaco et al., 2006;Valakos et al., 2008; Appendix S1). The environmental responses of the species were evaluated by biogeographic subregions as proposed by Coll et al. (2010) and Kougioumoutzis et al. (2017).

| Environmental data
Several ensembles of variables were selected based on their influence on the composition of the island biotic communities (Ding et al., 2006). The variables describe the geophysical, climatic and biotic characteristics of the islands, as well as their isolation relative to the mainland (i.e. continental territories) and surrounding islands (Table 1). The geophysical characteristics of the islands include the surface area (km 2 ) and elevation (which can be used as a proxy for habitat heterogeneity; Kougioumoutzis & Tiniakou, 2015), obtained from published data (Koster, 2005) and Google Earth Pro 7.3.2.5776 (Google LLC).
The isolation of the islands was assessed by calculating the shortest distance to the mainland and the average elevation around the 10 km perimeter of the island (Kreft et al., 2008). These data were generated from a raster-based geographical information system ( Figure 1) showing negative or positive elevational values relative to sea level (Becker et al., 2009). The island isolation relative to other islands within an archipelago was evaluated after constructing a connection network (Kougioumoutzis et al., 2017). This network was based on simple proximity relationships, linking the islands with their closest neighbour, the largest islands with each other and these in turn with the mainland. Network properties were used to estimate the connectivity, using simulated annealing (Carstensen et al., 2012). The connectivity indicates the relative importance of the nodes (islands) within the entire network, based on the estimation of the Zi parameter (Guimerà & Amaral, 2005). I also calculated the island radiality (i.e. the ease with which a node can reach other network nodes, which is determined using the geographical distance) (Borgatti & Everett, 2006). If the islands are close to each other, radiality values are expected to be high, whereas islands far from the main clusters will show lower values. Connectivity and radiality parameters were calculated with the rnetcarto (Doulcier & Stouffer, 2015) and igraph (Csardi & Nepusz, 2006) packages for the R environment (R Development Core Team, 2020).
The island environment was also characterized based on climate and vegetation. I used two variables to describe the climate, that is mean annual temperature and mean annual precipitation, which were obtained from the WorldClim 2 database (Fick & Hijmans, 2017).
Vegetation was characterized using the enhanced vegetation index (EVI), a proxy for primary productivity and food resource availability (Huete et al., 2002;Wei et al., 2019); EVI data were obtained for the mid-spring (April-May) period from the high-resolution Landsat 8 Collection Tier 1 32-Day EVI composite images (Chander et al., 2009). Ten islets had a surface area that fell below the raster resolution (250 m/pixel; Chander et al., 2009). Given that the EVI value approaches zero for those islets with a surface area of less than 500 m 2 (mean EVI = 0.016, n = 15), a zero value was assigned to these ten islets smaller than the pixel area. Human population density and snake presence were also included in the models. Human density because it is associated with the alteration of island ecosystems and greater penetration of mainland species (Escoriza, 2020a;Silva-Rocha et al., 2019) and snake presence because these reptiles are thought to play a role in moulding island lizard assemblages F I G U R E 1 Map of the study region, where the terrain elevation has been superimposed. Upper panel, Number of islands populated by endemic species, grouped within a radius of 50 km (red circles); lower panel, number of islands populated by species of mainland origin, grouped within a radius of 50 km (blue circles) (Brock et al., 2015;Pérez-Mellado et al., 2008). Human population density data were obtained from open-access censuses.

| Data analysis
Prior to analysis, variables with values of skew and kurtosis that suggested that they followed a non-normal distribution were logarithmically transformed (Lewis, 1977), generating a new distribution with skewness and kurtosis values near zero. All of the continuous variables were normalized. The response of the species was evaluated using binomial generalized linear models (BGLMs). The models were built using two dependent variables: the presence of endemic and mainland lizards. The relative contribution of the pooled variables to the models was evaluated using an automated procedure (Calcagno & de Mazancourt, 2010), ranking the models with the Bayesian information criterion (BIC; Schwarz, 1978). BIC was recommended over the Akaike information criterion when models are optimized for explanatory purposes (Aho et al., 2014). The importance of each variable was determined according to its model-averaged weighting based on the 100 best models. These analyses were conducted using the package glmulti (Calcagno & de Mazancourt, 2010) in R.
Species associations with the island environment were assessed using outlying mean index (OMI) analysis (Dolédec et al., 2000). This method describes the species responses by quantifying their ecological marginality, that is the distance between the species centroid and the average environmental conditions (OMI; Dolédec et al., 2000).
An OMI value approaching zero indicates a higher resemblance between the species position and the environmental average. The OMI analysis also provides an estimate of the niche width of the species (tolerance index). In this study, I used an extension of this OMI analysis, the within outlying mean index (WitOMI) approach, which was developed to evaluate the responses of species within spatial subsets (Karasiewicz et al., 2017). Following arcsine transformation, the distributions of the OMI indices were compared using a t test. These analyses have been carried out using the packages ade-4 (Dray & Dufour, 2007) and subniche (Karasiewicz et al., 2017) for the R environment.

| RE SULTS
The results obtained from the BGLMs show that the parameters that best explain the variations in the occurrence of the endemic lizards are mainland distance, annual temperature and mainland lizard species richness (Table 2). This association was positive and statistically significant with mainland distance (β = 0.673, p = 1.27e -8 ) and connectivity (β = 0.467, p = .0019), but negative and statistically significant with annual mean temperature (β = -0.578, p = 7.27e -7 ) and mainland lizard species richness (β = -0.384, p = 1.05e -6 ; The parameters that best explain the variations in the occurrence of the mainland species are the mainland distance and annual mean temperature (Table 2). This association was positive and statistically significant with island surface area (β = 1.198, p = 6.20e -6 ), annual mean temperature (β = 0.787, p = 1.12e -7 ) and EVI (β = 0.514, p = .0068), but negative and statistically significant with mainland distance (β = -1.292, p = 2.85e −11 ; Table 2). The analyses also showed a positive association with human population density, although it was not statistically significant (β = 10.213, p = .2164; Table 2).
The first two axes of the OMI explained 74.35% of the total inertia (axis 1 = 56.45%, axis 2 = 17.90%; Figure 3). The first axis described a gradient according to the island size and topography ( Figure 3). The second axis represented a mix of climatic conditions and geographical isolation, separating species from the northern Mediterranean from the southern, and those that occur on continental islands from those that occur on near-oceanic islands (Figure 3).
Overall marginality (WitOMI) was lower for the endemic island species (mean = 29.6 ± 30.6 SD), than for those from the mainland (mean = 51.8 ± 35.7 SD), resulting in significant statistical differences (t = 2.393, p = .0189; Table 3). The values for each species within subregions for niche parameters are shown in Appendix S2.
The dispersion of the WitOMI was also greater for the mainland species (Figure 4), indicating that they occupy a more variable niche than do the endemic species. Tolerance showed slightly higher values for the endemic species (mean endemics = 17.8 ± 15.3 SD; mean mainland = 17.5 ± 17.2 SD), but this difference was not statistically significant (t = 0.042, p = .9667; Table 3 and Figure 4).

| D ISCUSS I ON
The analysis, presented above, of the lizard populations on 552 Mediterranean islands reveals that the endemic lizards differ from the mainland species in their environmental preferences. In particular, two findings have implications that will improve our understanding of the spatial patterns of coexistence between these groups of TA B L E 2 General patterns of occurrence of lizards in the Mediterranean islands, assessed by automated GLM selection

F I G U R E 3
Plots showing the results of the OMI analysis. The left panel showed the position of the species and groups convex hulls, while the right panel showed the predictor variables, represented as vectors. In both panels, the origin represented the average environmental conditions. Red: endemics; blue: species of mainland origin that occur on islands lizards: (i) endemic species tend to occur on islands distant from the mainland, but well-connected within an archipelago, and (ii) the size of the island is more relevant to the presence of the mainland species than for the endemics. These findings imply that large islands are more easily colonized by mainland species than smaller islands (MacArthur & Wilson, 1967) and that island inter-connectivity is relevant only for those species that have been established in an archipelago for long periods of time. An alternative interpretation of these results would be that island endemic lizards have a greater capacity for transmarine dispersal, but this hypothesis is not supported by the colonization times estimated for peripheral islets (Rodríguez et al., 2013(Rodríguez et al., , 2014Terrasa et al., 2009). The larger islands are more easily colonized by mainland species for two reasons: first, the greater probability of being reached (due to the effect of the surface and also favoured by human transportation), and second, the greater capacity to support stable populations (Whittaker et al., 2008).
Most of the endemic species are well represented on the satellite islets, and in some cases, these islands constitute their last sanctuaries . With the possible exception of the larger species (Čerňanský et al., 2016), the low occurrence of mainland species on the islets is more likely to be a reflection of the difficulty of reaching these smaller islands in the first place, rather than of their inherent habitat limitations. The observation that accidentally translocated species become quickly established on the islets supports this hypothesis (Corti et al., 2014;Koren et al., 2011).
Mainland species are usually superior competitors and display Central Aegean Endemic 4 52.9 ± 39.8 6.0 ± 5.5

Mainland 11
53.2 ± 36.8 14.9 ± 12.9 Southern Aegean Endemic 1 2.8 29.9 Mainland 10 50.9 ± 32.3 9.9 ± 10.9 Note: Endemic: species with exclusive or mainly insular natural distribution; mainland: species with exclusive or mainly continental natural distribution. n, number of species. WitOMI, distance between the species' environmental centroid and the subregion environmental centroid. Tolerance, breadth of the species' niche. TA B L E 3 Niche indices per subregion (mean ± standard deviation)

F I G U R E 4
Violin-box plots for the niche indices, generated by WitOMI analysis. Red: endemics; blue: species of mainland origin that occur on islands aggressive behaviour and more efficient antipredatory responses (Corti et al., 2009;Nikolic et al., 2019). Therefore, the absence of mainland species in small islets is possibly a consequence of the lack of human activity (Corti et al., 2014;Grbac & Brnin, 2006).
Temperature is another environmental factor that is positively associated with the occurrence of mainland species. The most widespread mainland species (i.e. Chalcides ocellatus, Hemidactylus turcicus and Tarentola mauritanica) display thermophilous habits, as they evolved in the southern Mediterranean (Kornilios et al., 2010;Martínez Rica, 1974). These species were transported by humans in historical times and were able to colonize several Mediterranean archipelagos (Kornilios et al., 2010;Rato et al., 2011). However, their presence in the northern Mediterranean is more localized, possibly because of the effect of the lower temperatures (Delaugerre & Cheylan, 1992). niches. This means that these species are able to colonize many heterogeneous islands, highlighting their propensity for opportunistic implantation (Barreiros et al., 2010). In the case of at least one species, P. siculus, expansion has been correlated with the rapid decline of island endemic species, but adverse effects have also been postulated for the other species (D'Amico et al., 2018;Mayol-Serra, 1977).
Although islets are colonized by species of alien reptiles that are quite diverse phylogenetically (Escoriza, 2020b), those evolutionarily close to endemics, as using similar resources, are frequently associated with the collapse of native populations (Nikolic et al., 2019).

| CON CLUS IONS
The results of this study suggest that the populations of endemic island lizards could be threatened by the expansion of mainland lizards. Most of the endemic species occupy niche positions similar to the environmental average, but without using the entire range of available niches within the archipelago, particularly large islands with the presence of competitors and/or predators. Satellite islets play a key role as refuges for endemic species, and this role is supported by their inherent isolation. Large islands are more easily colonized by mainland species, thereby acting as reservoirs for secondary colonization.

ACK N OWLED G EM ENTS
The author thanks several anonymous reviews provided helpful comments that improved the quality of this manuscript.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13206.

DATA AVA I L A B I L I T Y S TAT E M E N T
The author confirms that the data supporting the findings of this study are available within the article and its supplementary materi-