Controls of plant diversity and composition on a desert archipelago

Study area: the Midriff Islands

The Midriff Islands in the Gulf of California, Mexico (Fig. 1) are an isolated, arid, and long-inhabited archipelago. This set of 14 islands that span the Gulf of California, or Sea of Cortés, from the Baja California peninsula to the Mexican mainland occur within the Central Gulf Coast subdivision of the Sonoran Desert (Shreve, 1951). The islands range in size from 0.2 km² (Isla Cholludo) upward to 1,223 km² (Isla Tiburón, the largest island in Mexico), reaching highest elevations of 1,316 m on Isla Ángel de la Guarda and 885 m on Tiburón (Felger, Wilder & Romero-Morales, 2012). They have been inhabited or visited by the sea faring and hunter-gatherer Comcaac (Seri people), or their ancestors or predecessors, for millennia (Felger & Moser, 1985; Bowen, 2000). Geologically, the Gulf of California dates to ca. six Ma (Ledesma-Vázquez & Carreño, 2010) and is one of the world’s most productive oceans (Douglas et al., 2007). Today the islands are a UNESCO World Heritage site, a Mexican natural protected area (DOF, Diario Oficial de la Federación, 2000), and since the 1950s uninhabited. For a comprehensive review of the Gulf of California islands, see Case, Cody & Ezcurra (2002). Collections were made under Mexican federal collecting permit NOM-126-SEMARNAT-2000 issued to Dr. Ezcurra.

Species-area model

Like many other scaling phenomena in biology, the relationship is usually described by a power function of the form S = kA^z, where S is the number of species on an island, A is island area, z is the power exponent, and k is a scale coefficient. In practice, the species-area curve is an accumulation function similar to those used to evaluate the completeness of herbarium collections (Soberón & Llorente, 1993) and is based on Preston’s canonical model (Preston, 1962). The theory on biological scaling has dramatically expanded since 1967 (Brown & West, 2000), but the basic species-area model, based on MacArthur and Wilson’s “equilibrium theory” of island biogeography is still one of the most robust models of ecological theory (see Gould, 1979 for a review).

To estimate the power coefficient z of the area–diversity relationship for the Midriff Islands, we used generalized linear models (GLM), a novel approach in species-area relationships. Because our dependent variable, the number of species on each island, is a count variable that is expected to have a Poisson error distribution, we used Poisson regression (also known as log-linear regression), in which the fit of the model to the data is measured as a log-likelihood function with a χ² distribution of the error (McCullagh & Nelder, 1989).

The derivation of our model is as follows (see Appendix S2 for an alternative derivation): The species-area model for islands S = kA^z can be rewritten in log-log form: log S = log k + z log A. This is a linear model where log k, the intercept, can be also expressed as a constant parameter b = log k. Raising the exponent of both sides we get: S = e^{(b+z log A)}, a log-linear function in which the count variable S is expressed as the exponent of the linear function b + z log A. The GLM algorithm identified the values of b and z that minimize the error of the fit. The error in turn is calculated as a log-likelihood deviance function: $\text{ε}=2{{\displaystyle \sum }}^{\text{}}{S}_i\text{log}\left(S_i/{\widehat{S}}_i\right)$ that is numerically convergent to the more familiar chi-squared statistic: $X^2={{\displaystyle \sum }}^{\text{}}\left[{\left(S_i-{\widehat{S}}_i\right)}^2/{\widehat{S}}_i\right]$ (in both equations, S_i is the number of species observed in island i, and ${\widehat{S}}_i$ is the number of species predicted by the model for island i). In short, we fitted a Poisson GLM to the species richness of the islands, using the logarithm of the area as the predictor. The slope and standard error of the fitted function gave us the value of z and its error.

Because in a Poisson model the variance is equal to the mean, each individual Chi-squared error term $\left[\left(S_i-{\widehat{S}}_i\right)/\sqrt{{\widehat{S}}_i}\right]$ is distributed as a normal-distribution deviate $\left[z=\left(x-\text{μ}\right)/\text{σ}\right]$ that can be tested against the critical values of the Normal curve (i.e., errors are Pearson residuals, see Duffy, 1990). Using this property, we were able to identify islands that had significantly more or significantly fewer species than predicted by the equilibrium model of island biogeography.

Predictors of species richness

After confirming that some islands depart significantly from the species-area model, we tested additional predictors of species richness for the archipelago (Table 1), using a step-wise process on the same GLM procedure. The predictors were: (a) presence of archeological artifacts (yes or none found) as a surrogate for prehistoric human presence, (b) habitat diversity (or number of distinct habitats; a number between 1 and 14), (c) topography (subtle, moderate, or rugged relief), (d) distance to the nearest mainland, and distance to the nearest island (in km), (e) island type (oceanic or landbridge), (f) seabird dynamics (presence/absence of nesting seabirds), and (g) precipitation. Some of our predictors were continuous variables while others were categorical.

Data for the predictors are as follows (full sources can be seen in Table 1): (a) The presence of humans on the islands was incorporated into our biogeographic analysis via archaeological investigations undertaken on all islands in question (Bowen, 2009). (b) While area and habitat diversity are correlated (Sfenthourakis & Triantis, 2009), and hence not truly independent, they may jointly have a significant additive effect (Hortal et al., 2009). For example, small islands can have a high degree of topographic complexity despite their size. In Sonoran Desert ecosystems, topographical diversity is strongly related to habitat and plant diversity (Búrquez et al., 1999). Following Hortal et al. (2009) we defined habitat diversity as the number of physiognomic vegetation types on each island as described by Felger, Wilder & Romero-Morales (2012). (c) Topography was taken into consideration by assigning each island to one of the three categories of relief based on highest elevation and our experience on the islands. (d) Isolation by distance was taken into account via an island’s proximity to the nearest mainland (Baja California peninsula or Sonora) and in several cases neighboring large island(s) that may serve as sources for immigration. (e) Geologic history of the island was incorporated via island type: landbridge or oceanic. (f) Given the preponderance and importance of nesting seabirds in the Gulf of California, especially on small islands (Velarde et al., 2005), we classified each island as a seabird island if there is annual presence of over ca. 5,000 breeding individuals and characteristic guano white-washed rocks and soil following Anderson & Polis (1999). (g) Because reliable and long-term climatological or meteorological data do not exist for the Midriff Islands, we used statistical extrapolations of long-term precipitation data from six coastal meteorological stations maintained by the Mexican government. Three stations were selected on the Gulf coast of Baja California (N–S: San Felipe, Bahía de los Ángeles, and Santa Rosalía, Ruiz Corral et al., 2006a, 2006b) and three from the Sonoran coast (N–S: Puerto Libertad, CONAGUA, 2010; Bahía de Kino, Ruiz Corral et al., 2005; and Guaymas, García et al., 1973). Precipitation for individual islands was calculated using linear interpolation between the nearest three stations to each island.

Island plant diversity was modeled using GLM in R (R Core Team, 2016). At each step of the regression analysis we tried all the variables that could be added, then selected the best fit using Akaike’s Information Criterion. We continued to add variables until no variable could significantly improve the fit.

Species composition

To examine patterns in plant species composition across the 14 islands, a checklist of the flora of the Midriff Islands (Wilder, 2014; Appendix S1) was converted into a 476 (species) × 14 (islands) presence-absence matrix. Both species richness and species composition analyses include non-native species (only 12 non-native taxa occur in the flora). A non-standardized principal component analysis (PCA) was conducted on this matrix. Jackson’s (1993) broken-stick test was used to determine significant axes. To interpret the floristic variation captured by each axis in terms of external driving variables, we regressed the site scores of each axis against the same seven predictors (a–g) used for species richness. Five small islands were omitted from the distance-to-the-source analysis because the majority of their flora is formed by widespread coastal halophytes that obscure other patterns.

Using Noy-Meir’s (1973) divisive polythetic classification method, we partitioned the species along the first two significant PCA axes into floristic groups with minimum intragroup and maximum intergroup variation, numerically searching for partition thresholds that minimized intragroup variance. To visualize the distributional similarities of the floristic assemblages obtained, the centroids of these groups in the first two PCA axes were then clustered into a hierarchical dendrogram using the procedure hclust in R with average-linkage clustering and Euclidean distance (R Core Team, 2016).

Following Baselga (2010) we performed an evaluation of floristic nestedness on the island data to test whether small islands contain a flora that is distinct from that of the large islands, or if, alternatively, the small island flora is a nested subset of that in the large islands. We also tested whether the flora of small islands, if nested, is a random subset of the larger islands or a non-random collection of large island species. For this analysis, we used Monte–Carlo simulations randomly extracting two vectors with a probability of containing a mean of 19 species per simulated island, and comparing the Sorensen similarity among them. We did this 1,000 times, and the mean similarity and standard deviation of the simulated island pairs was compared to the true mean similarity among small islands in our data set.

Finally, to test whether certain families were overrepresented in some islands, we used a X² analysis of the islands × families matrix with the number of species present in each family and each of the Midriff Islands. Probability values were subject to a Bonferroni correction to adjust for bias in multiple comparisons. Family classifications follow the Angiosperm Phylogeny Group IV system (Angiosperm Phylogeny Group, 2016) and recent work by Stevens (2019), reflecting current knowledge of evolutionary relationships.

Predictors of species richness

The species-area model (Fig. 2) yielded a Poisson regression line with r² = 0.85 and a slope z = 0.303 (s.e. ± 0.01). As is frequently the case in island ecosystems, the value of the estimated species-area exponent (0.303) was significantly higher than 0.263, Preston’s canonical value expected for random species accumulation in terrestrial ecosystems (t = 4.0, d.f. 12, P = 0.0009). The analysis of the residuals identified four islands with more plant species than expected based solely on their area: Tiburón, Dátil, Alcatraz, and Cholludo, all part of the Tiburón landbridge archipelago. Six islands were shown to have significantly less species than expected: Ángel de la Guarda, San Pedro Mártir, Partida Norte, Patos, Salsipuedes, and Rasa. Apart from Ángel de la Guarda, each of these species-poor islands can be characterized as bird guano islands. Furthermore, all of them, except Patos, are true oceanic, not landbridge, in their geologic origin. Finally, four other islands, San Esteban, Las Ánimas, San Lorenzo, and Nolasco fell within the predicted range of the species area curve (Table 2).

After fitting the effect of island area on species richness, four additional predictors accounted for a significant amount of the residual variation (Table 3): (1) habitat diversity, 34%, (2) seabird rookery, 23% (3) island type (oceanic vs. landbridge), 21%, and (4) topography, 14%. Together, all these factors accounted for 92% of the observed variation in species richness, and the residual variation (8%) was not significant, matching the assumptions of a random Poisson distribution. Two other factors, distance to the source and human presence, had no significant effect in accounting for the observed variation.

Species composition

The PCA analysis had two significant axes as determined by the broken-stick distribution model (Jackson, 1993). The first axis explained 41% of the overall floristic variation in the dataset and the second 19% (broken-stick expected values were 23% and 16%, respectively). The PCA simultaneously analyzes the relationship between islands, each with their component species (Fig. 3), or the distribution patterns of species among islands (Fig. 4). When islands are assessed, the first axis is strongly correlated with species richness (r² = 0.99; Fig. 3B), while the second axis displays a floristic-composition gradient in the large islands from west (Baja California, Ángel de la Guarda) to east (Sonoran mainland, Tiburón; r² = 0.74; Fig. 3C); the floristic composition of the small islands was independent of distance to the source. This second PCA axis was also correlated with precipitation (r = −0.61, P = 0.04) highlighting the fact that the Baja California flora west of the Gulf occurs in drier environments than the Sonoran flora in the east.

When the distribution patterns of the 476 plant taxa are analyzed, five distinct groups of species can be recognized from their loadings in the first two PCA axes (Fig. 4): small or single island species, generalist species, and large-island species can be distinguished along axis 1, while Sonoran-side species, and Baja-side species clearly separate along axis 2.

There was a high level of nestedness of the small island flora in the large islands, especially with Tiburón where Baselga’s (2010) β_NES was 0.68–0.80. Indeed, with the exception of two taxonomically suspect species (an undescribed Mammillaria on Dátil and Cholludo islands and Johnstonella grayi var. cryptochaeta in Salsipuedes, reported in Case & Cody (1983a) but not supported by herbarium specimens), all the rest of the small island flora is nested within the plant list of the large islands. Under the hypothesis of random colonization from the large islands, the Monte–Carlo test predicted a mean Sorensen similarity of 0.04 (s.e. ± 0.04) between small islands, while the true mean similarity between small islands is almost ten times higher (0.34 se ± 0.02). The significant difference (t = 6.7, d.f. 3, P = 0.003) between the observed similarity among small islands and the much lower value predicted under the random colonization hypothesis endorses the idea that the group of species colonizing small islands is not a random subset of the flora of the larger islands but rather part of a smaller group of species that are able to tolerate the harsh environmental conditions of small islands, like saline sea-spray and eutrophic guano substrates.

The Cactaceae showed a significantly higher richness on the smaller islands compared to the rest of the flora (X² = 71.9, d.f. 13, P < 0.0001). Most of the small bird-guano islands (Alcatraz, Cholludo, Las Ánimas, Patos, Rasa, Salsipuedes, and San Pedro Mártir) showed a very high proportion of cacti and halophytic succulents (as defined by Wilder, Felger & Romero-Morales (2008)) pooled together, which differed significantly from the overall pattern of species richness (X² = 44.8, d.f. 13, P < 0.0001). The family Polygonaceae had no species on small and medium islands and only one species (Eriogonum inflatum, the common desert trumpet) on large islands (Tiburón, Esteban, and San Lorenzo), but eight species on the otherwise species-poor Ángel de la Guarda, a significantly higher number compared to the other large islands (X² = 13.2, d.f. 3, P = 0.004).

Influences on species richness

Plant species richness patterns among the islands confirms the well-established species-area power function. The ecological heterogeneity of the islands, expressed both as the number of different habitats and topographic variation, accounted for the largest proportion of variation in species richness not accounted by area alone (Table 3). The combination of these two factors in supporting diversity has been long recognized (Darwin, 1860; Preston, 1962; Hortal et al., 2009), including prior analyses for these islands (Felger & Lowe, 1976; Cody, Moran & Thompson, 1983).

The geologic history of the Gulf of California is reflected in the island type: oceanic islands were either formed from the eruption of deep underwater volcanoes or were derived from segments of the peninsular crust that became separated from the peninsula by deep faults millions of years ago. In contrast, landbridge islands are really topographic mounds separated from the coast by shallow channels less than 100 m deep, which were connected to the mainland during most of the Pleistocene glaciations and only became isolated toward the end of the last glacial period, some 10,000 years ago (Graham, Dayton & Erlandson, 2003; Fig. 1). In agreement with the region’s geologic history, landbridge islands have significantly more species than oceanic ones.

An interesting case study is Isla Ángel de la Guarda, an extremely arid oceanic island rifted off the Baja California peninsula ca. two Ma (Nagy & Stock, 2000). Ángel harbors only 217 species, 54 less than predicted by our species-area model and significantly less than the 349 species in the similarly-sized landbridge Tiburón Island that receives almost twice as much rainfall (Felger, Wilder & Romero-Morales, 2012). The effects of isolation on Ángel de la Guarda are also evident in its relatively high level of endemic reptiles, rodents (Lawlor et al., 2002), and plants (Moran, 1983).

Although the geologic nature of the island (landbridge vs. oceanic) and xeromorphic gradient are strongly related to species richness, the distance to the mainland per se did not show a significant relationship with plant species richness. This was true even among oceanic islands (Table 3), likely because of the narrowness of the Gulf. The importance of distance in the Gulf of California is taxon-dependent and varies with colonization and persistence rates (Case & Cody, 1987). Taxa that have high colonization rates, such as birds, show few effects of insularity (Cody & Velarde, 2002), while mammals and reptiles are less dispersive and more often evolve unique forms in isolation (Soulé & Sloan, 1966; Lawlor et al., 2002). Plants, which also have high colonization rates and often possess long distance dispersal ability, show a relatively low level of endemism in these islands, except among Cactaceae (Rebman, 2002).

Finally, the presence of seabird rookeries explains a significant amount of the variation in species numbers: small guano islands harbor less species than similarly-sized islands that lack rookeries. Recent studies have shown that excess nitrogen can significantly reduce the diversity of terrestrial ecosystems (Simkin et al., 2016).

The relatively high value of the species-area exponent (z) is the statistical result of the existence of many small guano islands with low species richness. The species composition of these small islands is quite different from that of larger islands or from that of the mainland source (see next section). Porous land-sea boundaries facilitate the movement of nutrients from the ocean into the islands. The importance of marine inputs in regions with high primary productivity on small islands is accentuated by their inherent high perimeter-to-area ratio (Anderson & Wait, 2001; Talley, Huxel & Holyoak, 2006). Accordingly, marine subsidies (in this case guano), have a disproportionate importance as a driver of biological diversity on small islands (McCauley et al., 2012).

Influences on composition

The small bird islands in our analyses are linked to the productive waters of the Gulf of California through a long chain of ecological interactions. These isolated islands are free of predatory mammals (Lawlor et al., 2002), which combined with shallow relief provide ideal habitat for some of the largest aggregations of breeding seabirds in North America (Velarde et al., 2005). Seabird diets consist mostly of pelagic fish that occur in abundance in the Midriff region and feed on microscopic plankton in extremely productive upwelling zones (Velarde, Ezcurra & Anderson, 2013).

Seabirds deposit copious amounts of marine derived nutrient rich guano, high in nitrogen (N) and phosphorus (P; Hutchinson, 1950), in much greater concentration on bird islands than on non-bird islands in the Gulf of California (Anderson & Polis, 1999). The presence of elevated levels of N and P on seabird islands act as a strong filter that selects for specific plant species, as shown in the high mean Sorensen similarity value of small islands and the identification of a clustering of small/single island species in the PCA and species cladogram analyses (Fig. 4). It is primarily these, and only these, species that occur on the small islands, resulting in the depressed plant diversity observed on these islands. However, the plant species that can tolerate these nutrient loads, such as cacti and halophytic succulents, occur in greater proportion and remarkable abundances (see the insets of Cholludo and San Pedro Mártir islands in Fig. 1). The percentage of succulent and halophilous flora of bird islands in the Midriff region is relatively greater than that on large islands, suggesting a difference in establishment ability within the flora and providing an outstanding example of marine-based nutrients controlling terrestrial diversity patterns in small islands.

Analysis of the insular floras as a matrix reveals the importance of the proportional distance of the larger islands to either the Baja California peninsula and the Sonoran mainland, as well as the E–W precipitation gradient of the Midriff as a driver of species composition (Felger, Wilder & Romero-Morales, 2012). The majority of the species encountered on Ángel have peninsular affinities and are adapted to extreme aridity. Persistent cold-water upwelling and location of Ángel on the western edge of the North American monsoon, as well as the southeastern edge of Pacific-derived winter storms result in the western Midriff Islands being the most arid portions of the Gulf of California. In addition, ca. 30 plant taxa characteristic of the northern mediterranean region of the Peninsula reach their southern limit on or near Ángel (Moran, 1983). This is exemplified in our analysis by the preponderance of the buckwheat family Polygonaceae, a family that is strongly associated to the California Floristic Province, on Ángel and its near absence on Tiburón.

The human dimension

Humans have been a broad-scale determinant of species diversity on islands worldwide, often through habitat disruptions and species extinctions (Burney & Flannery, 2005) and cultural use (Bye & Linares, 2000; Heinsohn, 2003). Portions of the world thought to be pristine, without human modification, have been shown to be otherwise (Denevan, 1992; Gómez-Pompa & Kaus, 1992). While our analysis is limited to the presence/absence of archaeological elements on the islands, careful ethnographic work with the Comcaac (Seri people) on the coast of Sonora has helped illuminate their profound interaction with the island and coastal environments (Felger, 1976; Nabhan, 2003; Wilder et al., 2016). The consumption, utilitarian use, harvesting, and transport of plants across the region and between islands by the Comcaac and or their ancestors has been shown repeatedly (Felger & Moser, 1985; Nabhan, 2002).

The unique culture of the Comcaac grew and evolved in the Gulf of California, and is endemic to the region. Their language, Cmiique Iitom is a language isolate (Moser & Marlett, 2010) and their knowledge of the natural world is extensive (Felger & Moser, 1985; Marlett, 2014). They lived in and visited the larger Gulf of California region for thousands of years, with archaeological evidence of their presence on many islands in the Midriff (Table 1; Bowen, 2000, 2009). Yet, any environmental changes wrought by the Comcaac have been integrated with and are indistinguishable from the background ecological dynamics of the Midriff Islands. This is most likely because human population pressure was low and varied from ca. 180–3,500 individuals through time (Nabhan, 2002), limited by the scarcity of fresh water (Felger & Moser, 1985).

One manifestation of the non-altered environment of the Midriff Islands is the low frequency of non-native plant species (Felger & Moser, 1985; Felger, Wilder & Romero-Morales, 2012). Just 2.5% (12 species) of the Midriff Island flora is non-native (Wilder, 2014), relative to non-native floristic compositions of 14% among the Baja California Pacific Islands and 28% among the California Channel Islands (Ratay, Vanderplank & Wilder, 2014). The diminished presence of foreign cultures in the homeland of the Comcaac greatly reduced the opportunities for the establishment of non-native species.