Predicting species presence and richness on islands is important for understanding the origins of communities and how likely it is that species will disperse and resist extinction. The equilibrium theory of island biogeography (ETIB) and, as a simple model of sampling abundances, the unified neutral theory of biodiversity (UNTB), predict that in situations where mainland to island migration is high, species-abundance relationships explain the presence of taxa on islands. Thus, more abundant mainland species should have a higher probability of occurring on adjacent islands. In contrast to UNTB, if certain groups have traits that permit them to disperse to islands better than other taxa, then phylogeny may be more predictive of which taxa will occur on islands. Taking surveys of 54 island snake communities in the Eastern Nearctic along with mainland communities that have abundance data for each species, we use phylogenetic assembly methods and UNTB estimates to predict island communities. Species richness is predicted by island area, whereas turnover from the mainland to island communities is random with respect to phylogeny. Community structure appears to be ecologically neutral and abundance on the mainland is the best predictor of presence on islands. With regard to young and proximate islands, where allopatric or cladogenetic speciation is not a factor, we find that simple neutral models following UNTB and ETIB predict the structure of island communities.
Several interacting extrinsic and intrinsic factors may control the presence of species on islands. The equilibrium theory of island biogeography (ETIB) predicts that both island distance from continents and island size influence species richness, where extinction resistance is higher on larger islands and colonization rates increase on nearer islands [1–3]. When communities tend towards equilibrium, it is expected that colonization decreases and extinction increases with time as the environment becomes saturated with species. Additionally, the age, type and habitat complexity of an island may be predictive of species richness [3,4].
Intrinsic factors specific to the taxa themselves, such as the ability to disperse from the mainland to islands and maintain population sizes that are resistant to extinction, may also influence species presence [3,5–10]. Species able to disperse widely throughout the mainland and across the ocean may be better candidates for establishing populations on islands compared with those taxa that are poor dispersers or never maintain large population sizes. Therefore, all taxa from the mainland may not share the same immigration and extinction resistance abilities . Dispersal itself into new communities may be an evolutionarily conserved trait [12,13], which suggests that island communities should be phylogenetically clustered relative to the mainland source pool. Two studies examining phylogenetic structure of mammal communities on islands generally found mixed results, with overdispersion and underdispersion occurring in a smaller fraction of cases, indicating that idiosyncratic ecological processes mostly generated these communities [14,15]. Generally, given enough time, habitat stability and distance, both allopatric and within-island in situ diversification may also be predictors of island richness and degree of endemicity [3,16–18].
In contrast to phylogenetic assembly and diversification, communities on the mainland and species presence on islands may simply follows predictions from the unified neutral theory of biodiversity (UNTB) , where niche or phylogenetic differences have little bearing on community composition. Rather, as an extension of ETIB, neutral theory assumes that trophically similar groups of species are functionally equivalent with regard to niche, and demographically identical with respect to birth, death, immigration and speciation [19–21]. As a model of zero-sum ecological drift, UNTB predicts that a death in a community of a fixed number of individuals is replaced randomly by the offspring of any species, which effectively eliminates the roles of ecology or potential negative competitive effects experienced from all other species in the community [20,22].
Although criticized for being unrealistic and sometimes not fitting data [22–24], UNTB provides a null starting point for addressing presence on islands. The theory predicts then that for young and proximate insular regions, where allopatric or in situ speciation is irrelevant, abundance on islands simply reflects similar species-abundance structures on the mainland, which is consistent with the curvilinear version of extinction and immigration expectations from ETIB [25–27]. Therefore, if species differences with respect to colonization or phylogenetic history are negligible, it is expected that neutral mainland community abundances will predict species presence on islands; abundant taxa on the mainland simply have a higher probability of dispersal and successful colonization. In contrast, if there are measurable phylogenetic differences in dispersal, population size maintenance or competition, then island communities may be structured differently and not simply as a random sample from the mainland communities. In the situations where islands recently became disconnected from the mainland as sea level rose, extinction resistance might be more important than colonization ability for predicting presence on islands.
We ask if the presence of snake species on islands in the Eastern Nearctic (ENA) is better predicted by phylogenetic assembly rules, which could show evidence of habitat filtering or competition, or by ecological neutral theory given mainland community abundances. Snakes in the ENA provide a particularly interesting test of these hypotheses: (i) they are a species-rich group (n = 84 taxa), (ii) community richness responds to changes in environment, and (iii) communities are structured by both traits and local habitats . Importantly, trait filtering may be relevant for some clades like the watersnakes (Thamnophiini), which are associated with aquatic habitats [29,30] and therefore may be better adapted for dispersing over water for small distances to colonize islands. We used a complete phylogeny for ENA species to test for phylogenetically structured turnover, and we assembled species-abundance distributions (SADs) for 54 mainland communities with associated presence/absence on adjacent island communities to test for neutral community assembly. For these young insular systems where turnover of snake species from the mainland is likely to be high, we predict that community assembly will generally follow predictions from ecological neutral theory where mainland abundance predicts species presence on islands, and that phylogenetic structure will show little evidence of over- or underdispersion.
2. Material and methods
(a) Community data
We obtained published surveys for the presence of snakes on 54 islands from the ENA (figure 1; electronic supplementary material, tables S1 and S2). These islands probably originated after the Pleistocene, when they became disconnected from the mainland as sea levels rose (island distances from the mainland range from 0.32 to 19.6 km) or formed as barrier islands during the Holocene [32,33]. While it is important to consider that other nearby islands along with the mainland may act as colonization sources , we note that only 20% of the islands are closer to another island than the mainland (electronic supplementary material, table S2) and of these only five are closer to a larger island. All species on these islands are a subsample of those on the adjacent mainland (electronic supplementary material, table S1); there are no known endemic taxa on these islands and we have no evidence that species were introduced onto the islands or in the mainland communities.
To obtain presence and abundance data on the mainland, we surveyed counties on the mainland adjacent to each island using the museum record interface HerpNet (www.herpnet.org, 3–5 June 2014). From these, we obtained 18 295 occurrence records and found a total of 73 taxa in these coastal communities (out of 84 species found in ENA), and confirmed that all of these taxa were present according to previous sources [35–44]. We grouped mainland source communities into counties or combinations of counties adjacent to each island, resulting in 30 unique county combinations paired to each of the 54 islands (electronic supplementary material, table S1). We identified a total of 51 species distributed across these 54 islands. While it was not possible to standardize collecting methods to estimate local abundances across museum records, these samples are likely to represent the relative frequency of occurrence of each species placed into collections, have the advantage of detecting uncommon taxa, and should be suitable for testing UNTB.
(b) Richness tests and island biogeography
To understand the basic predictors of relative species richness on islands (species richness on each island/species richness on the adjacent mainland), we used the following variables: island area, straight-line distance from the mainland, number of habitats as generated using GAP Land Cover data (which includes land use and vegetation types at 30 m resolution ), latitude, longitude, the 1st, 2nd and 3rd principal components of the Bioclim variables covering temperature (Bioclim 1–2, 4–6, 8–11; removing correlated variables) and precipitation (Bioclim 12–19)  after centring. Using Bayesian model averaging (BMA), we assessed the relative importance of predictors over the entire model space, which includes the effect of each predictor variable assessed independently or cumulatively with other variables . The best-fitting model using a single parameter or combination of parameters was chosen as the one with the highest posterior probability distribution using the R  package BMS . We identified the relative importance of variables over the entire model space  and predicted which combinations of these variables yield the highest cumulative probability. To determine if overall richness is predicted by area, habitat numbers and latitude, we also fitted simple linear and typical species-area power function models (species richness = intercept × areaslope, S = cAz; see ) with latitude and habitat numbers as covariates, and ranked these tests using AIC.
(c) Phylogenetic data
To determine if island communities deviated phylogenetically from the adjacent mainland pool, we assembled the dated phylogeny of 144 species from the United States comprising all 73 taxa found on the mainland and island communities generated by Burbrink & Myers  for the entire United States. This phylogeny was well supported and dates were similar to those from previously estimated studies [50–52].
(d) Phylogenetic tests
We determined if presence on islands reflects overdispersion or underdispersion relative to the adjacent mainland source, the latter suggesting that traits from certain groups permit species to either colonize islands or resist extinction [28,53,54]. We note that underdispersion could be interpreted as competition, particularly in situations where competition results in unidirectional trends in morphology [55,56]. Here, we used a metric of richness-independent phylogenetic turnover, PCD , modified from the R package Picante . To contrast this with other metrics of phylogenetic turnover, we also used mean nearest taxon distance (MNTD or MNND in Picante), which is sensitive to concentrated changes at the tips of the phylogeny, mean phylogenetic distance (MPD), which is sensitive to concentrated changes deeper in the phylogeny, and the Dab statistic, which corrects problems in the former metrics where a difference of 0 in identical communities is occasionally not recovered [59,60]. We first determined if dissimilarities between islands and their paired mainland communities were different from an overall null, where PCD = 1.0 indicates that differences between the paired island–mainland communities were random.
We tested if PCD over all of these pairs was significantly different from 1.0 using a t-test. We also assessed the probability of under- or overdispersion by subsampling randomly 1000 times from the mainland communities to generate test island communities of the same size as our real island communities, and then calculated the overall phylogenetic turnover of this null community. We compared the z-score distance between our observed metrics and this null distribution of metrics for each island/mainland community pair, and estimated the probability that our observed estimates were significantly different from the null, which generates the typical standardized metric commonly used in other studies . We determined if the z-scores among the four metrics were correlated and whether they generally provide similar results.
Using simple shuffling techniques and calculating metrics like phylogenetic species variability (PSV) may increase type I errors, where overdispersion is more often predicted because of failure to account for allopatric speciation (closely related species do not coexist due to competition but rather failure to colonize the same area [61,62]). To account for this, we used the DAMOCLES model in R , which assesses colonization, extirpation and speciation to test neutral assumptions about community assemblage while controlling for allopatric speciation and diversification outside of the focal area. Here, using estimates of colonization and extinction for each island and adjacent mainland separately from the pool of Nearctic snakes , we simulated 100 datasets under null models of assemblage, estimated PSV, MNTD and MPD, and compared our observed values for these metrics for each community to these nulls. We determined if the z-scores among the three metrics were correlated and provide similar results (e.g. the standardized effects of these tests). We also compared the difference between PSV from each paired mainland and island community with the expectation that the overall differences were equal to 0 (using a t-test), indicating neither island nor mainland communities were more over- or underdispersed than the other.
We also applied BMA in the R  package BMS  to determine which variables alone or in combination predict PSV within islands or the mainland or phylogenetic turnover (PCD) between islands and adjacent mainland . We identify the most important predictors of phylogenetic assembly from the same set of variables to predict richness.
We also developed a statistical test in R to determine if PSV was biased by the presence of one or more taxa being present on the island, such that the absence of a single taxon could reduce PSV and PCD. This test detects outlier taxa by estimating the difference between PSV with and without the presence of each taxon over the absolute sum of all PSV differences across all taxa . All code developed for the island–mainland phylogenetic community assemblage and neutral predictions is available in the Dryad data repository (http://dx.doi.org/10.5061/dryad.dp157).
(e) Tests of local and metacommunity neutrality and island presence
We tested neutral models predicting ecological equivalence under UNTB for local and metacommunity assembly within our 54 mainland sites (reduced to 30 independent sites; 24 were surveyed from the same pool of mainland counties) and across these sites given that they may be connected to a broader metacommunity with variable immigration rates. Determining that the local communities cannot reject neutrality permits us to then make inferences about mainland species abundance and associated island presence . We used the approach by Harris et al.  implemented in the program NMGS (https://github.com/microbiome/NMGS), which applies Bayesian inference to fit models that distinguish and test for neutral community assemblage under the local and full multi-site model. This method, which is different from previous genealogical approaches [65,66], models UNTB as a hierarchical Dirichlet process (HDP) for large populations, where neutrality and many finite yet panmictic populations are linked by rare migration events. Using this UNTB–HDP approach, we produced the full posterior distribution over all parameters, Θ (biodiversity parameter) and Μ1:54 (all immigration rates [21,64]) for our real communities. From this, we generated neutral estimates from the local and the full neutral model using our fitted parameters from the observed data. We then compared the generated likelihoods under the local and full neutral models with our observed metacommunity estimates, which were used to produce a Monte Carlo significance test where the proportion of samples from the simulations that exceeded our observed values represent a pseudo-p-value. We generated model parameters and neutral simulations for each community and used a burnin of 25 000 samples, thinning every tenth estimate. We produced 50 000 samples for testing neutrality.
Given that mainland SADs were predicted by UNTB, we determined if local mainland species abundance is positively related to the presence of species on islands [25,27] using a mixed effects logistic regression model in the form of a generalized linear mixed model (GLMM) fitted by maximum likelihood. This determines if island presence or absence for each mainland species was predicted by the abundance of the same species on the adjacent mainland given the random effects of island differences (island identity) and species differences (species identity) using the standard logit (natural log of odds) transformation with binomial responses for binary (presence or absence) outcomes.
(a) Relative species richness
Using BMA, we found that size of the island was the only predictor of relative species richness (figure 2). This measure alone had the highest posterior inclusion probability (99%), suggesting that most of the posterior model mass contained variables that would include the area of the island; area was the only variable contained in the model with the highest posterior probability (14%). We also find that the power function model (S = cAz) given latitude best predicts these data (ΔAIC = 9.25) with an r2 = 0.34 and a slope (z) = 0.20.
(b) Phylogenetic assemblage models
All metrics that assess phylogenetic beta diversity were significantly correlated (correlation coefficient 0.27–0.73; p = 0.04–2.05−10), therefore we discuss results using PCD (for all statistics, see electronic supplementary material, table S2). Values for PCD among all island–mainland pairs cannot be rejected from the null hypothesis that the phylogenetic differences between these communities are random (mean = 1.04, s.d. = 0.16; t-test: t = 1.87, p = 0.07). Using phylogenetic assemblage models to compare our observed PCD estimates of phylogenetic turnover against null shuffling we found that none of the islands showed underdispersion and 5.5% showed overdispersion (figure 3). Using BMA, the highest posterior predictive model (3.5%) contained no predictor variables. This suggests that none of the measured variables were predictive of the mean differences between our observed and simulated PCDs (figure 2).
All metrics that assess PSV were significantly correlated (correlation coefficient 0.78–1.00; p = 8.20−11–2.20−16) and we therefore report the statistics for PSV (for all statistics, see electronic supplementary material, table S2). Differences among all island–mainland pairs for PSV were not significantly different from 0 (mean = −0.003, s.d. = 0.05; t-test: t = −0.40, p = 0.69), suggesting that islands or mainland communities are not more or less over- or underdispersed than one another. Only five communities showed that the removal of a single species of viperid (Crotalus adamanteus in three communities and Agkistrodon piscivorus in two communities) would significantly alter PSV. We also found using the DAMOCLES model, where a type I bias for overdispersion is reduced, that no communities were overdispersed on the islands, and only 7% were underdispersed on the mainland and on the island communities, though this may simply be that islands with reduced richness provide little statistical power to detect deviations from the null. None of the communities that were underdispersed in the mainland were also underdispersed on the adjacent islands. Island BMA shows that the highest posterior predictive model (15.6%) included only the predictive variable species richness (posterior mean 0.99), which demonstrates the effect of the relative size of the island community for predicting PSV (figure 2). In contrast, BMA for the mainland with the highest posterior predictive model (32.1%) indicated that species richness, location (latitude and longitude) and precipitation were most predictive of PSV.
(c) Neutral models and tests
Comparisons of likelihoods under the full and local neutral models suggest that the observed likelihoods (llh) were not significantly different from simulated neutral values (full p = 0.91, median simulated llh = −3206.74, median observed llh = −2656.01; local p = 0.871, median simulated llh = −3553.83, median observed llh = −3656.01; figure 4), indicating that both local and metacommunities are not significantly different from predictions under UNTB. We found that the GLMM for predicting island presence, given frequency of mainland taxa with both island and species random effects, to be preferred over those without random effects or with only a single random effect (ΔAIC = 144.21), yielding a significant positive coefficient for mainland species abundance associated with presence on islands (expressed as an odds ratio estimate = 1.78; lower and upper confidence interval = 1.44–2.11; p < 5.01 × 10−11; figure 4).
As expected under UNTB, we show that abundance on the adjacent mainland is the best predictor for the presence of snake species on 54 islands in the ENA (figure 4). Abundance structures on islands mirror those on the mainland for young islands and those particularly close to the mainland, which is consistent with predictions from UNTB and similar theories . Every log-scale increase in the abundance of a particular species on the mainland doubles the probability for the presence of that species on the adjacent island. Though the mainland abundance structures are consistent with UNTB, this does not mean that ecology fails to play any role with regard to the actual assemblage of species, particularly at a more regional scale. For example, both traits and phylogeny are important for structuring communities of snakes across North America . At least for islands in the ENA, where the time scale is more recent and area is within the fixed boundaries of islands, it is likely that particular species' presence on islands is related to mainland abundance given ecological neutrality.
While community structure may influence phylogenetic relationships among taxa locally, such that overdispersion is predicted under competition and underdispersion is predicted for trait filtering , we find little signal for this against null PCD simulations measuring mainland to island turnover (figures 2 and 3), though it is possible that competition does not exclude species with similar traits or close phylogenetic relationships . Additionally, it is possible that there is a lack of statistical power in areas with low species richness, though to our knowledge there are no studies that have examined the power of phylogenetic community assembly methods to detect types of dispersion given source and local community sizes. Our results contrast with those studies where traits that are important for colonization reflect strong phylogenetic signal [12,13].
Even though watersnakes (Thamnophiini) might be able to preferentially colonize islands given their predilection for aquatic habitats [29,30], we see no signal of filtering for any phylogenetic group. For watersnakes though, having a consistent source of freshwater on the islands may be more important for establishing populations. Dispersal has been shown to be limited within certain watersnakes , such as in the mangrove saltmarsh snake (Nerodia clarkii), a species associated with coastal habitats , though this species is the only natricine to have colonized the West Indies . It is also possible that certain unmeasured traits with low phylogenetic signal are important for colonization; phylogenetic measures of turnover would provide no signal. In light of the island rule, though, where morphology is known to vary in the unique ecological contexts of the islands and mainland , future studies would need to account for differences in traits from mainland/island community pairs prior to understanding if there are key traits for colonization.
When accounting for biogeographic history and allopatric speciation, we also fail to detect overdispersion, where it has been shown that using random species shuffling to produce null communities artificially inflates the probability of detecting overdispersion . Throughout the entire Nearctic, there is strong phylogenetic community signal associated with deep biogeographic structure in snakes . Trait filtering associated with changes in climate along steep ecological gradients was detected, though not at the local levels shown here, and not with regard to mainland and island dispersal. Similarly, Cardillo et al.  did not find strong evidence of filtering in mammals across islands. They note that results could be affected by the negative impacts of human activity on the coastal regions of the world, particularly in the southern United States . The possibility remains that coastal communities of snakes in the ENA may have been disrupted by human influence, therefore prehistoric reconstructions of communities may be important for establishing pristine assemblages on the mainland.
With respect to neutral models, ETIB is related to UNTB, though the former does not consider abundance structures [19,20,26]. We find that richness on islands, as a fraction of the mainland pool, is simply related to area and not, surprisingly, to habitat diversity or mainland distance (figure 2). This has been demonstrated in other studies [6,72] where area (and not island habitat diversity or distance) is the greatest predictor of richness. Better quantification of species associations for a particular habitat may help determine if habitat is yet a driver of richness when controlling for area , though it is likely that colonizing species exist as generalists on islands. Unfortunately, habitat specificity of snakes on all these islands is largely unknown. It is also expected that, with greater distances (and lower migration), mainland and island SAD will differ markedly, particularly when cladogenetic speciation drives richness . While this cannot be tested in the ENA, because no unique species occur on these islands, this area of study remains open for studying snakes in the tropics, where examples of cladogenetic diversification occur, such as in Madagascar, with approximately 99% endemicity of colubroid snakes .
While the presence of species on islands can be predicted given neutral theory, it is uncertain whether these island or mainland communities are at equilibrium with respect to colonization (or speciation) and extinction, and, in particular, the zero-sum assumption of UNTB. Equilibrium is an important part of UNTB, and the probable lack of equilibrium has been used to criticize UNTB or show that predictions are inconsistent given time and other factors [22,73]. However, we point out that many community studies are consistent with UNTB and that although equilibrium may not have been demonstrated in all cases, it is likely that some relaxation of this constraint and the zero-sum assumption is permitted [74,75] or is unnecessary to predict SAD [76,77]. Therefore, immigration is simply a sampling process, and if the mainland community has a neutral species-abundance distribution, then so should the adjacent island .
We note that with regard to UNTB, local abundance can be linked to dispersal ability . Therefore, if dispersal ability drives local mainland abundance, then it is possible that it also accounts for island presence. Related to this, we show a significant random effect given species identity in our GLMM, suggesting that while there is not a discernable phylogenetic assembly effect, certain taxa are more likely to colonize islands given abundance or migration ability. For instance, we find that racers (Coluber constrictor), one of the most abundant mainland species, were found on 74% of the islands surveyed, and the seven species found on more than 30% of the islands are from phylogenetically distant groups of colubroids, including colubrines, lampropeltinines, natricines and crotalines. Furthermore, snakes on several islands may have been present when the area became disconnected from the mainland, which excludes the numerous barrier islands that formed independently from the mainland and discounts the effect of persistent natural disasters over the Holocene. In those cases, extinction resistance, which is again related to abundance, would be the strongest indicator of species presence on the islands and produce essentially the same outcome as predicted by our GLMM models.
Finally, several authors have examined the connection between UNTB and expectations on phylogenetic trees [79,80]. It is possible to link these concepts at the level of individuals in communities through models which account for selection, given that neutrality is not likely over large time scales . We examined UNTB at the community level, and separately assessed if phylogenetic community structure exists with respect to over- or underdispersion. The way forward would be to model predictions about assembly using abundance, UNTB and phylogeny simultaneously. This would therefore require estimation of biogeographic area of origin, modes of speciation, ecological community interactions and ecological drift over time. Phylogeny along with community structure could then be tested against both null UNTB predictions and against ecological drift, which would ultimately help shed light on processes operating across time scales from local community assemblage to regional diversification and macroevolutionary deep-time processes.
All research followed the College of Staten Island Institutional Animal Care and Use Committee. Only databases and previously published literature (Burbrink & Myers ) were used for this study.
All data and code for this paper are included in the Dryad data repository (http://dx.doi.org/10.5061/dryad.dp157) and the electronic supplementary material.
All authors helped conceive the study and plan strategies for analyses. F.T.B. wrote code, generated the database, ran computational analyses and wrote the manuscript. E.A.M. helped write the manuscript. A.D.M. generated GAP Land Cover data and helped write the manuscript. R.A.P. contributed to database development and helped write the manuscript.
We have no competing interests.
This project was supported in part by L. Clampitt, D. Rosenberg and S. Harris, and a US National Science Foundation grant (DEB 1257926) to F.T.B. and R.A.P. (DBI-0905765 and DEB-1441719).
We thank K. Krysko, T. Kahn, W. Gibbons, M. Wead, J. Boundy, S. Christman, K. Dodd, C. Lechowicz and C. Sheehy for help assembling island communities. For useful discussions regarding UNTB metrics and DAMOCLES, we thank the scientists that were down from day one: L. Harmon, J. Rosindell and R. Etienne.
- Received July 16, 2015.
- Accepted October 26, 2015.
- © 2015 The Author(s)