Species–area relationships and additive partitioning of diversity of native and nonnative herpetofauna of the West Indies

Abstract To evaluate the regional biogeographical patterns of West Indian native and nonnative herpetofauna, we derived and updated data on the presence/absence of all herpetofauna in this region from the recently published reviews. We divided the records into 24 taxonomic groups and classified each species as native or nonnative at each locality. For each taxonomic group and in aggregate, we then assessed the following: (1) multiple species–area relationship (SAR) models; (2) C‐ and Z‐values, typically interpreted to represent insularity or dispersal ability; and (3) the average diversity of islands, among‐island heterogeneity, γ‐diversity, and the contribution of area effect toward explaining among‐island heterogeneity using additive diversity partitioning approach. We found the following: (1) SARs were best modeled using the Cumulative Weibull and Lomolino relationships; (2) the Cumulative Weibull and Lomolino regressions displayed both convex and sigmoid curves; and (3) the Cumulative Weibull regressions were more conservative than Lomolino at displaying sigmoid curves within the range of island size studied. The Z‐value of all herpetofauna was overestimated by Darlington (Zoogeography: The geographic distribution of animals, John Wiley, New York, 1957), and Z‐values were ranked: (1) native > nonnative; (2) reptiles > amphibians; (3) snake > lizard > frog > turtle > crocodilian; and (4) increased from lower‐ to higher‐level taxonomic groups. Additive diversity partitioning showed that area had a weaker effect on explaining the among‐island heterogeneity for nonnative species than for native species. Our findings imply that the flexibility of Cumulative Weibull and Lomolino has been underappreciated in the literature. Z‐value is an average of different slopes from different scales and could be artificially overestimated due to oversampling islands of intermediate to large size. Lower extinction rate, higher colonization, and more in situ speciation could contribute to high richness of native species on large islands, enlarging area effect on explaining the between‐island heterogeneity for native species, whereas economic isolation on large islands could decrease the predicted richness, lowering the area effect for nonnative species. For most of the small islands less affected by human activities, extinction and dispersal limitation are the primary processes producing low species richness pattern, which decreases the overall average diversity with a large among‐island heterogeneity corresponding to the high value of this region as a biodiversity hotspot.

Despite its broad applicability, there are several limitations of the SAR. First, for some organisms, samples from fixed areas are difficult to obtain, and species richness can only be expressed in units of sampling effort (Crist & Veech, 2006). Second, differences in the relative abundances of species are often ignored (Gotelli & Colwell, 2001). Third, as we demonstrate in this study, species are generally not homogeneously distributed among islands, and the heterogeneity in biodiversity within (α) and among (β) islands cannot be revealed by the SAR. What's more, β-diversity among islands is only partly explained by island area; thereby, SAR alone cannot quantify how much of the total β-diversity is due to area (β area ) and how much is due to other factors (β replace ; Golodets, Kige, & Sternberg, 2011;Zajac, Vozarik, & Gibbons, 2013). However, the comparison of the diversity within (α) and among (β) islands and the contributions made by area (β area ) and other factors (β replace ) are important for strategic conservation planning because a low α-diversity with a high β-diversity suggests that species assemblages are heterogeneous among islands and species are often endemic to individual islands, while a high α-diversity with a low β-diversity indicates that species assemblages are homogenous and species within each island are a subsample of the same species pool (Francisco-Ramos & Arias-González, 2013); a low β area with a high β replace suggests that factors such as in situ speciation and human introduction may have a greater role in influencing patterns of β-diversity, while a high β area with a low β replace indicates that species richness varies in a more predictable manner determined by factors such as habitat diversity, carrying capacity, and extinction/immigration dynamics as envisioned by MacArthur and Wilson (1967), Lomolino & Weiser (2001), Triantis & Bhagwat (2011).
Area can only explain a portion of species richness. Other factors, such as competition, dispersal, colonization, and speciation, also play important roles in influencing species richness (Crist & Veech, 2006).
But as the increase in study interests on small islands (Lomolino, 2000;Triantis et al., 2006) and the improvement of calculation methods were represented by the applications of piecewise regression and information theory (Dengler, 2010;Gentile & Argano, 2005;Matthews, Steinbauer, Tzirkalli, Triantis, & Whittaker, 2014), scholars have argued that the species-area pattern comprises of more than one SAR, whereby processes operating at different spatial scales lead to different Z-values (Gao & Perry, 2016;Lomolino & Weiser, 2001;Losos & Schluter, 2000;Morrison, 2014;Rosenzweig, 2004). Lomolino and Weiser (2001) and Rosenzweig (2004) even proposed three biological scales of species-area curve with three corresponding dominant processes of species addition and Z-values ranges: (1) Stochastic extinction forces structure insular communities on small islands with Z-values ranging from 0.10 to 0.20; (2) extinction/immigration dynamics on islands of intermediate size with Z-values ranging from 0.25 to 0.45; and (3) speciation on relatively large islands with Z-values higher than 0.60. And this proposal has been proved by Gao and Perry (2016), who recently detected the small island effect (SIE) of herpetofauna of the West Indies using a three-segment piecewise regression approach.
Thereby, we hypothesize that the overall Z-value is an average of the three Z-values of different scales, and the Z-value calculated from intermediate to large islands could lead to an overestimation of the overall Z-value.
However, there are several imperfections in previous studies concerning SAR. First, small islands are ignored, and most of these studies covered only a portion of the region, fewer than 200 islands in each case. Second, the Power model is used by default, without the evaluation of other possible candidate models. Third, the comparison of the diversity within and among islands and the contributions made by area and other factors remain blank. Fourth, although humans now become a major geological and environmental force (Corlett, 2015), studies for nonnative species are still insufficient. Here, we evaluate biogeographical patterns of native and nonnative herpetofauna in the West Indies, on a scale not previously attempted, comparing (1) SAR models, (2) Z-values, and (3) the heterogeneity in biodiversity within and among islands. The aims of the study were to investigate the following: (1) Whether sigmoid instead of convex models provide a better performance; (2) whether the Z-value is overestimated in previous studies and how Z-values vary among groups; and (3) whether area plays the same role toward explaining the among-island heterogeneity for native and nonnative species.

| Study area and data
The West Indies comprises over 3,000 islands, cays, and emergent rocks belonging to three main island groups: Bahamas, Greater Antilles, and Lesser Antilles. We derived complete herpetological species lists for each island from Powell and Henderson (2012), who recorded more than 1,000 species on 749 islands. We digitized islands using base maps in ArcMap 10 and ArcGlobe 10 (ESRI, Redlands, CA, USA), including not only the 749 islands included in Powell and Henderson (2012) but also hundreds of small explored islands that have no herpetofaunal species, for a total of 1,668 islands varying in area by over 10 orders of magnitude, from 3.9 × 10 −5 km 2 to 1.1 × 10 5 km 2 ( Figure 2). The resulting map was projected by a UTM_18N coordinate system with WGS_1984 datum.
Power and Exponential are by far the best known models, whereof Power model is the most frequently applied to species-area data in the literature. Recent work indicates that the shape of SAR curves in arithmetic space is often sigmoid rather than convex and has an upper asymptote, and thus attempts to model the SARs using various functions have emerged. Dengler (2009) listed as many as 23 possible models, among which eight are commonly used and summarized by Guilhaumon et al. (2010). We investigated eight candidate models (Power, Exponential, Negative exponential, Monod, Rational, Logistic, Lomolino, and Cumulative Weibull), including both convex and sigmoid models, to analyze SARs in each group. The Akaike's information criterion (AIC) was applied as a criterion for model selection (Burnham & Anderson, 2002). We also calculated Akaike weights (ω) for all models to evaluate each model's probability of providing the best explanation of the data. Moreover, we calculated the first and second derivatives of the Lomolino and Cumulative Weibull to test the shape of the maximum-likelihood curves fit for each taxonomic group, and we also applied a bootstrapping approach which resampled the data points 10,000 times and gave bootstrapped estimates of parameters and the probability that inflexion point occurs within the range of observed island size.
We also analyzed the linear function of the Power model (logS/ logA) to compare Z-values among different species groups and compare the Z-value of all herpetofauna obtained here with that calculated by Darlington (1957). Because the logarithm of zero is undefined, and any data transformation such as log(S + 1) biases C-and Z-values when comparing among studies (Williamson, 1981). To avoid biasing the parameters, only islands with recorded species were analyzed. Although removing observations is not favored, we believe sacrificing some data in exchange for an unbiased result was worthwhile (Russell, Clout, & McArdle, 2004).

| Additive diversity partitioning
Alpha diversity is the diversity of a single island or location; β-diversity is the number of species in a region that are not observed on an island due to multiple localities, and γ-diversity is the total diversity of a region. Because typically α-diversity is measured with respect to some fixed spatial scales, rather than a scale that varies by several orders of magnitude, we replaced α-diversity with the average diversity of islands in our dataset (hereafter referred to as "a" for brevity) and βdiversity with among-island heterogeneity (hereafter referred to as "b" for brevity). We used additive diversity partitioning to quantify the heterogeneity in biodiversity by comparing the diversity within (a) and among (b) islands and by comparing the contributions made by area (b area ) and other factors (b replace ). In the additive approach, diversity can be explored across spatial scales (Gering & Crist, 2002), and γdiversity (regional scale) is partitioned into the sum of the average diversity of islands (a) and among-island heterogeneity (b). Because a, b, and γ-diversity are measured using the same units, their relative importance can be quantified (Crist & Veech, 2006).
Although understanding the spatial scale at which diversity is generated is important, we also must know the reasons that cause one or more species to be missing from an island. When a species is missing from an island, one reason might be that the island is not large enough.
So, we used additive diversity partitioning combined with SAR to partition b into b area , which represents the average difference between a and the diversity predicted for the largest sample (S max ) and b replace , the average number of missing species that are not explained by area.
In order to unify the calculation method, we estimated S max using the Lomolino model (we use Lomolino because Cumulative Weibull did not provide parameter estimation for nonnative crocodilians) for each group.
We performed all analyses using R 2.15.2 (R Development Core Team 2012) and used the mmSAR package (Guilhaumon et al., 2010) for modeling the SARs. Species richness showed a phase of rapid rise across small-island areas followed by a subsequent flattening toward an asymptote. The model with lower AIC value indicates a stronger evidence for being better over the others. However, according to Burnham and Anderson (2002), 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 10 or more have no support. So, in our study, the model with the lowest AIC and models with ∆AIC no higher than two were all considered equally the best; for example, for amphibian species group (Table 1), the second best-fit model had a ∆AIC of 30, which was >2, so only one relationship (Cumulative Weibull) was selected for this group. However, for all species and all nonnative species groups, the second best-fit model had a ∆AIC ≤ 2, so two relationships were considered equally for these two groups. Exponential and Monod were never selected as best models; Negative exponential was selected only once; Rational and Logistic were selected three  Table 1).

Species richness increased as island area increased
The Lomolino regressions displayed sigmoid curves for native amphibians, native frogs, and nonnative turtles, but convex curves for the other taxonomic groups, whereas the Cumulative Weibull regressions displayed sigmoid curves only for native amphibians and native frogs, but convex curves for the rest of the groups (Table 2).
Moreover, the Lomolino regressions had a higher probability that an inflexion point occurs within the range of observed island size than the Cumulative Weibull regressions for each taxonomic group (Table 2).   (Triantis et al., 2012).
According to the additive diversity partitioning for each group ( Figure 5), a explained only 0.12%-0.67% of the variation in island species richness (mean = 0.25%). On the contrary, b explained more than 99% of the variation in island species richness. Thus, in this region and for these taxa, b ≈ γ-diversity.
We also compared b area and b replace for each group ( In the native species groups, the proportion of b area to the total b ranged from a low of 20.9% in native snakes to 49.7% in native crocodilians, with an average of 32.5%. In the nonnative species groups, the proportion of b area to the total b ranged from a low of 3.9% in nonnative snakes to 13.7% in nonnative amphibians, with an average of 8.7%. Thus, the area effect played a much less important role in nonnative species than in native taxa.

| DISCUSSION
We found the Cumulative Weibull and Lomolino models were better than the traditional Power and Exponential models to our data.
Compared with the C-value for the main taxonomic groups summarized by Triantis et al. (2012), the overall C-value in this study was high among vertebrate studies. The high C-value demonstrates that amphibians and reptiles require less space to sustain viable populations than other vertebrates, given the small body size and distinctive physiology of many species. For instance, frogs in the genus Eleutherodactylus (Eleutherodactylidae), which comprise the dominant frog fauna of the West Indies, contains terrestrial-breeding frogs that lay eggs on land or tree leaves, and these eggs later hatch into miniatures of the adults, bypassing the tadpole stage (Hedges, 1993), and that reproductive mode can occur in a cave, on a mountain top, or high in a tree without direct dependency on water; this greatly enhances viability in a water-deficient and space-limited environment.
If small islands are excluded from the calculation, Z-value will only represent the SAR from intermediate to large scale. Thereby, the exclusion of small islands could lead to an overestimation of the overall Z-value. Compared with the Z-value ranges pointed by MacArthur and Wilson (1967) and Rosenzweig (1995), the negative value of −0.01 for nonnative crocodilians may due to the small number of islands (n = 9), an inadequate sample size for robust linear modeling (Chase & Bown, 1997). However, the remaining Z-values (0.02-0.20) were also very low.
Over 1,300 species have been recorded in the region, with as many as 407 on the single largest island; however, a explained only 0.12%-0.67% of γ-diversity. The small a (0.81) indicates low species evenness, and this is because most of the islands across a large spatial range are small, and over 900 (about half) are not known to support any herpetofaunal species. The low Z-values and small a may partly be due to the fact that the "zero" islands (islands with no recorded species) are removed when calculating Z-values but included when calculating a. Moreover, the low Z-values and small a paradoxically reveal that (1) the reduction in area leads to a significant loss of species especially on small islands, and (2) (Powell & Henderson, 2012). The trade in live animals (e.g., pets and food) also has played an important role Powell et al., 2011). Second, many small islands have no recorded species, reducing the average diversity, whereas other small islands harbor relatively many species, lowering Z-values especially for nonnative species. At least some of these discrepancies likely result from human activities, which continuously supply new colonists to some islands but also transform previous land-use types into anthropogenic biotopes such as cultivation and settlements (Sfenthourakis, 1996), hosting some herpetofaunal species that coexist with humans (Henderson & Powell, 2001;Raxworthy & Nussbaum, 2000).
The increase of Z-values from lower to higher taxonomic groups is an inevitable result of the positive SARs. Because when a lower taxonomic group is combined with another one, species increments become larger and larger as the increase in island size, leading to a higher resultant Z-value.
Herpetofaunal species may be removed with a reduction in area (Sfenthourakis & Triantis, 2009;Triantis et al., 2006). This situation is particularly detrimental to habitat specialists (Sfenthourakis & Triantis, 2009). Moreover, many very small islands are subject to episodic instability and catastrophic events, such as devastation by storms or inundation by tidal surges, which can lead to total species extinction (MacArthur & Wilson, 1967). Extinction event rather than colonization is thought to be the primary process producing species-area patterns (Losos, 1996;Perry, Rodda, Fritts, & Sharp, 1998;Rand, 1969).
However, in this study, we found colonization can still be an important process on small islands, and this conclusion is consistent with that of Russell et al. (2004), who suggested that human activities can be a decisive factor affecting species richness.
Area explained only around 22.8% of the total b. The proportion of b replace to the total b was as high as 77.2%, suggesting rapid withinisland speciation is a main source of new species in this region, for example, the adaptive radiation of anoles (Losos, 2009). Our findings are consistent with those of Losos and Schluter (2000), who indicated that very few species present on small islands cannot also be found on nearby large islands. Apart from within-island speciation, human introductions, although minor, have become a new mode of entering the region for some species, many of which are not native to the West Indies, enlarging the among-island heterogeneity. For instance, Anolis carolinensis (native to USA) has arrived on Anguilla with the development of tourism (Eaton, Howard, & Powell, 2001).
Area is even weaker at explaining the among-island heterogeneity for nonnative herpetofaunal species, because nonnative species are introduced via human activities that might be directly caused by the economic development instead of island area (Powell et al., 2011). Normally, the larger the island, the larger human population and richer natural resources there would be, so economic development may have a positive correlation with island size. However, this is not always the case. For example, the US-Cuban trade embargo strongly increases Cuban economic isolation (Helmus et al., 2014), lowering the estimated S max and therefore decreasing area effect toward explaining the among-island heterogeneity for nonnative species. On the other hand, larger islands have lower extinction rate, higher colonization, and more in situ speciation than smaller islands, maximizing the estimated S max for native species.
Suitable herpetofaunal habitats are deteriorating on many (maybe most) islands as a consequence of human activities Powell & Henderson, 2012). Consequently, the probability of alien species becoming established is increased, and Powell et al. (2011) strongly recommended that an early detection and monitoring system for alien species be established in this region. If current trends continue, as in other regions of the world (e.g., Perry et al., 1998;Whitaker, 1973), we suggest that a number of endemic West Indian herpetofaunal species will be extirpated from large islands and might survive solely (if at all) on tiny offshore islets. A Z-value of around 0.3 obtained by Darlington (1957) suggests if the area is reduced by 90%, then the number of species it supports will be halved. However, our study demonstrated that Z-value is scale dependent, so the loss rate of species against area reduction is also scale dependent rather than constant across scales, and native species on large islands are most sensitive to area reduction.
Thereby, we also recommend that conservation priority should be given to large and habitat-rich islands with the most species-rich community to maximize the number of endemic species preserved.

| CONCLUSIONS
The Cumulative Weibull and Lomolino models are reported to be sigmoid (Dengler, 2009;Guilhaumon et al., 2010;Lomolino, 2000;Simaiakis et al., 2012;Tjørve, 2009). However, in this study, we F I G U R E 4 A comparison of the Z-value outcomes among 748 islands conducted in this study, seven islands used by Darlington (1957), F I G U R E 5 (continues)