Frontiers of Biogeography Species–area relationships of the Aegean, a comparative approach between six taxa

a small area on they Biogeography, herptiles, snails, isopods, tenebrionids and chilopods) and 20 different models. The aim was to determine which model(s) perform better for each taxon and also to compare the z and C parameters of the power model between animal groups, which are the only model parameters to date that have been linked with biological processes. We compared the relationship across different taxa for the entire archipelago and for the exact same islands, in two subgroups with similar paleogeographic history and environmental conditions in the central and eastern Aegean. For the taxonomic groups that were examined a strong correlation between the number of species and area was found, except for chilopods and herptiles. Although there is no universal best model for the SAR of the Aegean, the power model performed better for invertebrates, whereas concerning vertebrates there was more ambiguity in the shape of the relationship.


Introduction
Biological research on islands has been playing a crucial role in the deeper understanding of fundamental characteristics in ecology, evolution and biogeography.
Islands are considered to be extremely important in biogeography as they have easily delineated geographical borders, and are isolated and simple systems, making them ideal "natural laboratories" where biological phenomena, models and theories can be observed and tested Lomolino 1998, Whittaker andFernández-Palacios 2007).
One of the oldest, most ubiquitous patterns that has been recognised in ecology is the increase in species richness (S) with increasing sampling area (A): the species-area relationship (SAR). Although the SAR has been troubling scientists for more than 150 years (de Candolle 1855, MacArthur and Wilson 1967, Connor and McCoy 1979, He and Legendre 1996, it has broad applications in ecology and conservation and has been used as the basis for protected area design, prediction of species extinctions resulting from the loss of native habitat and estimating regional diversity from smaller-scale sample data (Guilhaumon et al. 2010). Hence, conservation biologists frequently rely on the species-area relationship to predict changes in species diversity resulting from habitat loss, and to develop strategies for conserving biological diversity in spatially limited reserves and fragmented ecosystems (Brooks et al. 1997, Smith 2010, Gerstner et al. 2014, Halley et al. 2014, de la Sancha and Boyle 2019. Generally, the usage and development of mathematical models in biology have as main objectives the understanding of biological phenomena and the prediction of patterns and characteristics of complex biological systems. Consequently, because of the great complexity of biological systems there is no perfect model that can include every parameter that affects complex ecological systems (Komineas and Harmandaris 2016). So, it is desirable to work with manageable models which maximise generality, realism, and precision toward the overlapping but not identical goals of understanding, predicting, and modifying nature (Levins 1966). Specifically for biogeography and in species-area relationship studies models are fitted to species-area plots for descriptive, explicative and predictive purposes. The identification of curve shape and the comparison of areas and studies constitute descriptive purposes. Explicative purposes -which are predominately the primary goals of this study-refer to the understanding of how observed or expected patterns in nature affect curve shape and the biological interpretation of parameters and their values. Finally, predictive purposes include extrapolation of total species numbers, prediction of extinction caused by fragmentation and habitat loss, and the identification of species hotspots (Tjørve 2009).
The SAR is commonly described by the power model, which in its logarithmic form is given by logS = logC + zlogA (where S= island species richness, A= island area, and z and logC are fitted parameters representing the slope and intercept of the model, respectively). It has been noted that fitting data to the log-log transformation yields only approximate estimates of the parameters of the power function, and may in fact produce significantly different estimates of z to the nonlinear model version (Conor and McCoy 1979, He and Legendre 1996, Williams el al. 2009, Tjørve and Tjørve 2017. During the past decade a series of studies have examined the form and shape of the SAR (Tjørve 2003, Dengler 2009, Williams et al. 2009, and more than 20 functions have been proposed to describe it, of which, however, the power function remains the most commonly used. On the contrary, analyses have often demonstrated substantial uncertainty in selecting the best SAR model for a given dataset (Stiles andScheiner 2007, Guilhaumon et al. 2008). Even though there is no universally best model, metaanalyses from true and habitat island datasets have shown the power model to provide the best general model from the 20 SAR models tested (Triantis et al. 2012, Matthews et al. 2016a. In the present work we consider (true) islands as bodies of land within a matrix of water, i.e. oceanic islands, continental-shelf islands and inland waterbody islands (Whittaker and Fernández-Palacios 2007) and thus we focus exclusively on island species-area relationships. This distinction is important because the processes producing SARs in true and habitat islands could be influenced in different ways (Rosenzweig 1995, Matthews et al. 2016a. Furthermore, fauna and flora in different types of true islands have been derived and shaped by different mechanisms (e.g. vicariance and dispersal). For the most part, on oceanic islands organisms arrive through dispersal and the majority of the species evolve from few ancestral species. On the contrary, on continental-shelf islands biodiversity and speciation is the product of vicariance but also the constant influence (dispersal) from the adjacent continental areas (Rosenzweig 1995).
The location of our study, the Aegean archipelago is of great biogeographical and evolutionary interest, as it accumulates numerous factors that make it exceptional for testing biogeographical and evolutionary theories and concepts (Sfenthourakis and Triantis 2017). It comprises a large number of islands and islets with a variety of geological and paleogeographical relationships with the neighbouring islands and the adjacent land regions. Also, it is geographically placed between three different continents, Europe, Asia and Africa. As a result, it has continuously been impacted and influenced by these three species sources and the assembly of life on its islands is dominated by species of European, Asian and African origin, which combined with the endemic species that have emerged on the islands, mainly due to their isolation, make the Aegean one of the biologically richest areas in the Mediterranean (Panitsa et al. 2018). The fauna (and flora) of the Aegean has emerged and diversified due to a complex mix of vicarianistic and dispersal events, as well as the long-lasting human presence, which has continuously interacted and impacted with the Aegean environment for more than 10,000 years (Poulakakis et al. 2014). Recently, Hammoud et al. (2021) highlighted that in terms of processes affecting species richness patterns, continental archipelagos differ fundamentally from oceanic systems because episodic connections with the mainland have profound effects on the biota of land-bridge islands.
A series of studies have examined SARs for multiple taxa in the Aegean (e.g., Watson 1964, Mylonas 1982, Sfenthourakis 1996, Welter-Schultes and Williams 1999, Dennis et al. 2001, Triantis et al. 2005 Frontiers of Biogeography 2021, 13.4, e52929 © the authors, CC-BY 4.0 license 3 Panitsa et al. 2006, Iliadou et al. 2014), but few have considered different models aside from the classic power model (Fattorini 2002, Simaiakis et al. 2012a, Sfenthourakis and Panitsa 2012, Kaloveloni et al. 2018. Although Kagiampaki (2011) studied the relationship focusing on plants and approached differences between them and several animal groups, and Fattorini et al. (2017) compared among five animal taxa in the Aegean archipelago, no study has compared the relationship between different animal organisms in the exact same subset of islands, i.e. where each taxon SAR comprised exactly the same islands.
In this framework, the aim of the present paper is to approach the species-area relationship for six different animal groups in the Aegean archipelago, namely breeding land birds, herptiles (amphibians and reptiles), land snails (including slugs), land isopods, tenebrionid beetles, and chilopods (henceforth birds, herptiles, snails, isopods, tenebrionids and chilopods). We fitted 20 models that have been proposed for SAR and calculated the best among them that described the relationship for each taxon using an information theoretic framework. Considering that most of the models we tested have merely statistical meaning, in other words, their parameters have not been linked to any biological processes, we also focused on the z and C parameters of the log transformed power model for each taxon, which remain the only ones that have been attributed biological explanation and significance.
Furthermore, in an effort to predominantly focus on the biological characteristics of each taxon that potentially influence and determine its SAR (e.g. the taxon's dispersal ability or overall ecology), we investigated the species-area relationship for different taxa on the exact same islands, specifically on two island subgroups that have common paleogeographical history and environmental conditions, the central Aegean for all the taxa, and the eastern Aegean for birds, herptiles and snails. These comparisons provide valuable insights since they allow us to concentrate on how the shape of the SAR for different taxa vary for the same islands of a specific complex, while reducing the noise that occurs from SARs that derive from different islands, which inevitably have different intrinsic characteristics. So, by comparing the SARs for different organisms between the same islands, aspects that are known to influence both the model selection and the z and C parameters of the power model, such as distance from species pool, the ratio of maximum to minimum island area (A max /A min ), number of islands (Triantis et al. 2012) are minimised or eliminated. Thus, this approach allows us to compare z and C values of SARs for different taxa not only in the entire archipelago, but also for SARs that have derived from the exact same islands.

Data Extraction
The numerous biogeographical, ecological or even taxonomical studies that have been conducted on the islands of the Aegean archipelago have resulted in a fairly notable data production and recording of the fauna for a considerable number of taxa on the Aegean islands. Hence, we collected data for species richness for 11 different taxonomical groups in the Aegean from published papers, books and doctoral theses. Regarding snails, data were also extracted from the collection of the Natural History Museum of Crete (NHMC). Moreover, reptiles and amphibians were grouped together as herptiles, because of their relatively low number of species. Taxa with data for less than 30 islands (e.g. butterflies, grasshoppers, bees and mammals) were not included in our analysis, as well as plants since they constitute an entire kingdom and hence have a disproportionately larger species number comparing to all taxa included in the study, which is documented to affect numerous aspects of the SAR (Triantis et al. 2012, Matthews et al. 2016a. It is evident that not all islands are equally studied for all taxa, therefore, the number of islands we considered in the analyses for each animal group was different. Island areas were always measured in km 2 , because C values change according to the unit used to express island surface (see Fattorini et al. 2017). In total, we analysed data from 150 islands, with area ranging from 0.05 km 2 to 1636 km 2 , for six animal groups, birds, herptiles, snails, isopods, tenebrionids and chilopods (Table S1).
In order to examine and compare the SAR for different taxa between the exact same islands that belong to Aegean island subgroups with fairly common paleogeographical history and environmental conditions, we took the intersection of the datasets for the central Aegean islands (Cyclades and Astypalaia) for all six taxa, and for birds, herptiles and snails for the eastern islands (Dodecanese, Samos and Ikaria) ( Fig. 1, Tables S2a and S2b).
Crete and its satellite islands and islets were not included in the analyses, although they geographically belong to the Aegean archipelago. Firstly, as Crete has been a discrete and isolated island for around 5 million years (Dermitzakis 1990), it could easily considered to be a continental fragment Fernández-Palacios 2007, Sfenthourakis andTriantis 2017) and along with its satellite islands constitutes a discrete island complex, with Crete serving as the main species pool for the surrounding islets. Secondly, its area is more than 5 times bigger than the second largest island in our dataset (Rhodes), which heavily influences the analysis and it is common practice these islands to be excluded from datasets and analyses (see Gao et al. 2019). The island of Evvoia was also not included in the dataset because, even though there are scarce data for some taxa, it is a severely undersampled island for most of the animal groups. We also opted not to include in the analyses, islets smaller than 0.05 km 2 .

Statistical Analyses
We compared 20 SAR models ( Table 1)  approach that can account for uncertainties in inferring the SAR, allowing the investigator to perform inferences while incorporating variability in both model selection and parameter estimation (Guilhaumon et al. 2008). Model residuals were evaluated for normality using the Lilliefors extension of the Kolmogorov normality test and for the homogeneity using a correlation of the residuals with the model fitted values (Matthews et al. 2019a). The fit of a model was considered to be satisfactory if both of these assumptions were met. Model performance was compared using the Akaike information criterion (AIC) and when necessary (n<40), its correction for small sample size (AIC c ; Burnham and Anderson, 2002). The smallest AIC value represents the best model for a given dataset, however, all models with AIC differences (∆AIC) between zero and two have equal support, whereas models that differ from the best model by between four and seven have limited support and those that differ from the best model by a value of ten or more, practically have no support (Burnham and Anderson, 2002). For the multimodel comparative analyses we only used the nonlinear implementation of the power model, but we also fitted the logarithmic form of the power model to each dataset using standard linear regression in log-log space to be able to make comparisons with previous studies (following Matthews et al. 2016). As the logarithmic form of the power function (Arrhenius, 1921) is the most frequently applied form for fitting SARs, it remains one of the few functions for which biological significance has been assigned to model parameters, and has a proposed, even if debated, theoretical basis (Preston 1962, Connor and McCoy 1979, Rosenzweig 1995, Martin and Goldenfeld 2006. Statistical analyses were performed in R (R core team, 2019), using the "sars" package (Matthews et al. 2019a), which allows the computation and comparison of 20 different mathematical models that have been proposed to describe the species-area relationship. Furthermore, to avoid statistical artefacts in the comparison of the power model parameters, we also performed unpaired t-tests for the fitted z and C values of the power regression lines of each taxon pair that was compared (see Trichas et al. 2008, Simaiakis et al. 2012b).

Aegean archipelago
Considering the entire Aegean archipelago and the log transformed power model, area explained 27 to 87% of the variance in species richness, concerning chilopods and isopods, respectively (Table 2a). Specifically, in chilopods species richness can only weakly be linked with area increase as even the best models according to AIC do not explain the data adequately ( 2 adj R <0.1 and p>0.05). For all the islands of the archipelago we found that z values for birds, herptiles and tenebrionids (0.26, 0.27 and 0.27) were statistically significantly higher than those of snails, isopods and chilopods (0.18, 0.18 and 0.12). On the other hand, C values ranged from 2.3 to 13.2 for herptiles and snails respectively, with herptiles having significantly lower C values and snails significantly higher from every other taxon (Tables 2a and 2b).
In addition to the classic log-log power model, we also plotted the multimodel averaged model, which is constructed weighting the fitted values of each of the 20 models by a model's AIC c weight (wAIC c ), a procedure that accounts for model selection uncertainty (Fig. 2).

Comparison of different taxa in two Aegean island subgroups
In the subset of central Aegean (Fig. 3), z values ranged from 0.16 to 0.4 (Table 3a) but they did not have statistically significant differences (Table 3b), due to the high standard error of the slope value and the relatively low number of islands (n=15). On the other hand, C values ranged from 1.5 for herptiles to 11.8 for isopods, with only herptiles showing statistically lower C values (Table 3b).
Concerning the eastern Aegean islands (Fig. 4), the power model parameters were similar for birds and snails, especially the z values. However, in herptiles, z values were significantly higher, whereas C values were significantly lower respectively from the other two animal groups (Table 4a and 4b).
Furthermore, the models that fit the data best -ΔAIC c <2 (following Burnham and Anderson 2002)-for each taxon in both islands, subgroups are presented (Tables 3a and 4a) and it is important to note that the best models for the specific subgroups for each taxon differ from those of the entire archipelagos (Table 2a).

Model
Abbreviation No. of parameters Equation

Model shape
Asymptotic

Discussion
Our results suggest that for the entire archipelago of the Aegean and the log transformed power model, the relationship between species and area was strong for snails, isopods, tenebrionids, and birds. In contrast, for herptiles and especially chilopods it is quite weak, given that area explains less than 50% of the variance in richness. Therefore, this could imply that, for these two taxa, other factors are more important than area in determining species richness, especially considering that chilopods in the Aegean are a well-studied group in both sampling and taxonomic level and thus the

Model comparison
Taking into account the evaluation of the different models, an individual model could not be described as the best to universally describe the SAR in the entire archipelago across all taxa (Table 2a and Table  S2). As a general pattern, we could assume that the power model seems to be performing better for most invertebrates, namely isopods, tenebrionids and also bees (from Kaloveloni et al. 2018), whereas for vertebrates there is more ambiguity and, especially for herptiles, sigmoid or models with an asymptote seem to perform better (in agreement with Gao and Perry 2016), possibly due to their relatively smaller total species number and/or lower dispersal ability.
Another interesting finding is that the best model for a given taxon in the entire archipelago is not necessarily among the best to describe the relationship for the different island subsets, which could potentially be attributed to intrinsic factors of the island group that is studied such as A max /A min or the number of islands. It is well documented that the linear model performs better as the number of islands included decreases, which is also observed for the invertebrates in Cyclades (Table 3a), and sigmoid models perform better for higher A max /A min ratios (Triantis et al. 2012, Matthews et al. 2016a). In addition, this could imply that different factors and underlying mechanisms possibly influence and define the SAR in the different island complexes even concerning the same taxon. In any case, this draws attention to the fact that model selection could be ambiguous even when it comes to a specific taxon in different subgroups of an archipelago. In our opinion, this further supports the evaluation and comparison of different models when studying SAR, especially when the aim is to produce guidelines and extrapolations for biodiversity management and protection in specific regions, since the translation of such insights into improved conservation guidance remains an important element for biogeography.

Biological significance and meaning of the power model parameters
A series of studies have tried to infer a biological meaning or explanation to the parameters of the power model and until today the discussion remains fruitful, since this constitutes an urgent matter for both pure and applied biogeography (Preston 1962, MacArthur and Wilson 1967, May 1975, Connor and McCoy 1979, Gould 1979, Coleman 1981, Martin 1981, Sugihara 1981, Rosenzweig 1995, 2004, Harte et al.1999, Lomolino 2000, Lomolino and Weiser 2001, Drakare et al. 2006, Fattorini 2007, Fattorini et al. 2017, Williams et al. 2009, Triantis et al. 2012, Matthews et al. 2016a, 2019b. A major critique is that interpretation of speciesarea curves derives from the post facto and ad hoc nature of the inferences and correlations drawn from them. This combined with the fact that the intercept parameter has been mainly overlooked as a quantity deserving biological or statistical explanation, or as a basis for biological inferences (Connor andMcCoy 1979, Gould 1979) and that C values show an increasing trend over scales from small to large, has led to scepticism that any biological significance could be attached to it and eventually to recommendations that C and therefore both parameters should be simply viewed as fitted constants devoid of specific biological meanings (Conor andMcCoy 1979, He andLegendre 1996). Another practice that probably creates more problems than it solves in the efforts to infer biological significance to the power model parameters is the comparison of C and z values from a mix of island SARs and species accumulation curves, which are two very different types of species-area curve (see Scheiner 2003, Matthews et al. 2016b).
The majority of SAR studies have mainly focused on the interpretation of z, the slope of the log transformed relationship. Z values are generally considered indicative of the processes establishing species richness and composition patterns (Triantis et al. 2012). Several factors have been theorised to affect the z value including, among others: isolation, with distant archipelagos having higher z values (Rosenzweig 1995), dispersal ability, with mobile taxa having lower z values (Wright 1981), scale of sampling, trophic rank, with species high up in the food web reporting higher z values (Holt et al.1999).
In the entire Aegean archipelago, tenebrionids, birds and herptiles have similar z values, significantly higher than for isopods and snails, which in turn have significantly higher C in comparison with the other taxa. At first glance, the higher slope values for herptiles and tenebrionids could be expected from their low dispersal ability, while for birds it could possibly be attributed to their highest place in the trophic rank. On the other hand, it is interesting and challenging to the idea that z values are mainly shaped by isolation and dispersal abilities as birds and herptiles, arguably the best and worst disperser respectively, have similar z values and, simultaneously snails and isopods, undoubtedly less mobile taxa than birds, demonstrate lower z values.
The intercept of the power model (C) has drawn significantly less attention than the slope (Gould 1979), with two main patterns having been reported for its values. First, logC values decrease progressively from inland to continental-shelf to oceanic systems (Triantis et al. 2012). With increasing distance from the possible species pool, dispersal is reduced and thus fewer organisms will be able to sustain viable populations through the rescue effect (Brown and Kodric-Brown 1977). Thus, in the most isolated islands, fewer species are expected than for the continentalshelf and inland archipelagos (MacArthur and Wilson 1967). Furthermore, logC values generally increase from vertebrates to invertebrates, which is consistent with the general expectation that vertebrates require more space to sustain viable populations than invertebrates. Therefore, C is considered to be indicative of the realised carrying capacity of the system per unit area (Triantis et al. 2012). Our results only partially agree with this pattern since herptiles in the Aegean have undoubtedly the lowest C in all the comparisons, whereas birds on the other hand do not demonstrate higher C values than tenebrionids (with which they have similar z) nor chilopods.

Comparison of z and C for different taxa in the two island subgroups
Within this frame of reference, we fitted and compared the estimated parameters of the log-log power model for two different island complexes in the central and eastern Aegean.
In the central Aegean islands, all taxa excluding herptiles (z=0.4), had similar z values, ranging from 0.16 to 0.24 (Table 3a) with none of the differences being statistically significant (Table 3b). This could be a useful reminder that what might appear to be notable differences in the power model parameters (i.e. Δz=0.24 and n=15), do not always represent statistically significant differences and should be approached cautiously. On the other hand, C values ranged from 1.5 for herptiles to 11.7 for isopods, with only herptiles showing statistically lower C values than all other groups (excluding chilopods). However, for most taxa, the relationship between species richness and area was not strong (R 2 <0.5), which complicates the comparisons (Table 3a). This could be partially attributed to the rather small number of islands (n=15) and to the fact that in order to compare between the same subset of islands, analyses included only islands where data for all 6 taxa were available, since for specific taxa in the entire Cyclades complex the correlation coefficient has been found to be notably higher (R 2 >0.8) for both isopods and birds (Sfenthourakis 1994, Simaiakis et al. 2012a. In the eastern Aegean islands we compared three taxa that differ in their dispersal abilities and in their biology in general snails, land birds and herptiles. Despite these differences, the z values were similar for birds and snails. On the other hand, for herptiles z and C values are significantly higher and lower respectively, with C being actually an order of magnitude lower than the other groups (Tables 4a and 4b). Following a traditional approach to the interpretation of the power model parameters we could point out that the lower values of C for the herptiles are expected, considering they are generally organisms that require numerous resources to be established on an island or ecosystem (Triantis et al. 2012), while their higher z value could potentially be attributed to their generally lower dispersal abilities (Williamson 1988, but see Aranda et al. 2012). However, as we intensely advocate for the model parameters to be examined and evaluated together and that separate interpretations of z and C should be rigorously discouraged, we concentrate upon the significantly inversely proportional deviation of the parameter values in herptiles from all other taxa, including birds, the other vertebrate taxon examined. These differences between herptiles and birds, also indicate that such a subdivision of organisms into broad groups (vertebrates), that several global scale analyses have implemented (Triantis et al. 2012, Matthews et al. 2016a in an effort to describe general patterns of variation in C and z values should be done with caution. Furthermore, on a more general note our results are suggesting that it is risky and potentially misleading to use a SAR model or parameter fitted to a specific taxon to infer patterns for a different (even closely related) taxon.
Nevertheless, regardless of statistical and conceptual shortcomings, analysis of the dependence of SAR slopes on taxon, environmental conditions and spatial scale have been proven useful, as differences in the z values convey information to underlying processes and their relative importance in a particular case (Drakare et al. 2006), and have undoubtedly resulted in our better understanding of island biogeography. However, we believe that the interpretation of the biological meaning of the power model parameters should generally be done with extreme caution, since it merely constitutes a correlation of a fitted parameter with complex biological processes that usually is not been inferred by a statistical/mathematical test. Besides, conclusions that have been drawn from separate evaluation of the two parameters should be taken with a grain of salt and ideally should be re-tested and possibly revisited. Recently, structural equation models (SEMs) were used to infer causality from different factors to the parameters of the power model further pinpointing, among others, the need to consider logC and z in tandem (Matthews et al. 2019b). Despite the recent advances in our knowledge about SAR, general consensus as to how specific factors and mechanisms contribute to SAR is still lacking among different taxa, environmental heterogeneity and conditions, spatial and temporal scales (Matthews et al. 2019b). Such a deeper understanding is crucial for SAR to eventually become a more powerful tool for describing, predicting and conserving biological diversity on islands and other isolated.

Conclusions
In the present study our aim was twofold. Firstly, to compare the power model's parameters for different taxa, especially between the exact same islands in two subgroups of the Aegean archipelago. In this respect, we got more information from the intercepts than from slopes (in agreement with Fattorini et al. 2017), and from the eastern rather than the central islands. As a general trend, snails and isopods display the highest intercept values whereas herptiles have by far the lowest. This could probably be a reflection to the area that these organisms need to establish viable populations and to their carrying capacity. On the other hand, the herptiles have the highest slopes which could be linked with their low Frontiers of Biogeography 2021, 13.4, e52929 © the authors, CC-BY 4.0 license 11 dispersal ability, but also to the fact that they have the lowest C values and number of species compared to every other taxon in the study. Secondly, to approach the species area-relationship of six animal taxa in the entire Aegean archipelago, a region where the SAR has been studied thoroughly, using a multimodel and information theoretic approach. In this regard we concluded that there is no universally best model to describe the SAR in the Aegean archipelago. Hence, in our opinion, given that the fit and comparison of different SAR models have been made computationally easy (Guilhaumon et al. 2010, Matthews et al. 2019a, the choice of the best model or combination of models for a specific situation and problem, there is no sound reason not to be preferred. Especially when the goal of the study is more applied, such as to set appropriate baselines for accurately assessing and protecting biodiversity on a specific isolate system and taxon.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Table S1. Islands and species number for each taxon included in the study Table S2a. Islands and species number for each taxon in the central Aegean Table S2b. Islands and species number for each taxon in eastern Aegean Table S3a. ΔAIC of models considering the entire Aegean archipelago Table S3b. ΔΑΙCc of models in the two Aegean island subgroups

Author Contributions
LM, MM, KV conceived and designed the study, LM, KV collected the data, LM analysed the data, LM, MM, KV discussed and interpreted the data, LM led the writing, all authors read and approved the final version of the paper.