For seabirds, food supplies and nest sites are largely driven by oceanographic gradients and island habitats, respectively. Research into seabirds’ ecological roles in insular ecosystems is crucial to understanding processes that structure seabird nesting assemblages. We examined the influence of island physiography and oceanographic factors on the spatial variation in α and β-diversity of nesting seabird assemblages.

Location
Southeastern Pacific Ocean.

Taxon
Birds

Methods
We compiled data from 53 seabirds breeding on 41 coastal and oceanic islands using different sources: our field records, online databases, environmental reports, and literature. We used generalized linear models (GLM) to describe the effect of island physiography (area, elevation, and isolation) and oceanographic factors (surface temperature, salinity, and primary productivity) on seabird species richness (α-diversity). We applied multivariate GLM to test the effects of physiographic and oceanographic predictors on species composition (β-diversity). We used Jaccard dissimilarities on species occurrences per island to calculate β-diversity partitioned into turnover and nestedness. Polynomial models allowed us to model these metrics against geographical and environmental gradients and so analyze patterns in seabird β-diversity across spatial scales.

Results 
Species richness was highest in Galápagos, Pitcairn, and Rapa Nui. Changes in seabird α-diversity across islands were determined by island area and distance to South America but not by oceanographic variables. Physiographic and oceanographic factors were significant in determining β-diversity. Changes in β-diversity were mostly due to species replacement (β-turnover) across three major island Systems (Galápagos Archipelago, Chilean coastal islands, and oceanic islands of the southeastern Pacific). The contribution of β-nestedness was restricted to small scales (within archipelagos).

Main conclusions
Physiographic and oceanographic factors explain species diversity of seabird assemblages on islands of the southeastern Pacific. Oceanographic variables did not affect species richness but significantly influenced species composition. Change in species composition reflects gradients across three marine biogeographical realms: Temperate South, Eastern Indo-Pacific, and Tropical Eastern Pacific. The low degree of species nestedness may reflect multiple evolutionary origins.

Study area

Our study area extended longitudinally from 130°44'W to 70°31'W and latitudinally from 1°40'N to 38°22'S. The dataset included 42 islands of six archipelagos: Pitcairn (4 islands), Rapa Nui (4), Desventuradas (3), Juan Fernández (3), Galápagos (13), and Chilean coastal islands (15). Thus, the islands of our dataset have varying degrees of isolation from other islands and the mainland and are distributed within different biogeographic provinces (Spalding et al., 2007; Figure 1). Most of the oceanic islands are located on the Nazca Tectonic Plate, except for the Pitcairn Archipelago which is located on the Pacific Plate. These islands are all part of chains of volcanic islands and seamounts.

Environmental variables

We compiled information on island characteristics from the literature and online databases. After removing highly collinear variables (i.e., Pearson correlation coefficient > 0.8), we considered seven environmental variables as predictors in the analyses: island area (km²), island elevation (m), distance to the mainland (km), human density (individuals per km²), sea surface salinity (PSS), sea surface temperature (ºC), and primary productivity (g.m^-3.day^-1) (Appendix S1). Island area, elevation, and distance from the mainland were obtained from the literature or estimated using tools available in ‘R’ (R Core Team, 2016) and ‘Google Earth’ (www.google.com/earth). We derived data on human density from public databases provided by national institutes of statistics and censuses (see metadata in Appendix S1). All oceanographic variables were extracted from the raster maps available in the ‘Bio-ORACLE’ online database (Assis et al., 2018; Tyberghein et al., 2012). The raster resolutions of these maps were approximately 9.2 km at the equator. We used a 50 km radius buffer around the center of each island to calculate the average values of each oceanographic variable using the R package ‘raster’ (Hijmans, 2017).

Seabird assemblages

We compiled data on nesting occurrences of 53 seabird species (Appendix S2) in coastal and oceanic islands of the southeastern Pacific from the literature, online databases, and unpublished reports. Most of the consulted literature was published in the past 30 years, although we also included older publications to reference confirmed reports on some islands (see metadata in Appendix S2). For the majority of the Chilean islands, we used information on seabirds from our field records from multiple expeditions conducted over 19 years (1999-2018) on Chilean coastal islands, and five years (2013-2018) on Chilean oceanic islands (i.e., Desventuradas, Rapa Nui, and Juan Fernández). All of the seabird information was organized in a binary matrix of species occurrence (1 = presence, 0 = absence) for each island (Appendix S2).

Data analysis

Seabird α-diversity - Hypotheses 1 and 2

We measured the α-diversity of seabirds of each island by quantifying species richness, measured as the sum total of nesting seabird species on each island. Before the analyses, we used the Shapiro-Wilk test to check the normality of all environmental variables (i.e., island characteristics and oceanographic variables). Environmental variables that differed from the normal distribution were adjusted by applying either square root or log₁₀(x + 1) transformations. We analyzed the relationships between seabird species richness and environmental variables by fitting generalized linear models (GLM) based on Poisson distribution and log-link function. Model selection and averaging were performed using the R package ‘MuMIn’ (Barton, 2018). The function dredge was used to perform a model selection routine based on the lowest second-order Akaike Information Criterion (AICc) and applied to 128 different model subsets. We applied the function model.avg to perform full model averaging, considering the best model subsets (AICw > 0.05). We plotted the partial residuals associated with each environmental variable to depict their effect on species richness. The final model considered the environmental variables island area (log₁₀ transformed), island elevation (square root transformed), distance from the mainland (square root transformed), human density (square root transformed), salinity, superficial primary productivity, and sea surface temperature. Finally, we used residual-distance correlograms and the Moran’s I test routines in the R packages ‘ape’ (Paradis, Claude, & Strimmer, 2004) and ‘ncf’ (Bjornstad & Cai, 2018) to assess the spatial independence of the GLM results.

Seabird β-diversity - Hypotheses 3 and 4

We applied the routines in the R package ‘betapart’ (Baselga & Orme, 2012) to calculate the partitioned β-diversity, which separates the turnover and nestedness-resultant components. The β-turnover component reflects species replacement, while β-nestedness reflects differences in species richness. These components are derived from the dissimilarity index chosen to describe β-diversity, which is frequently referred to as β-total. We used Jaccard dissimilarities on species occurrences per island to calculate β-diversity. The relationship between β-diversity and geographical and environmental distances was analyzed using linear models based on a third-degree polynomial fit. Geographical distances were calculated using the distm routine in the R package ‘geosphere’ (Hijmans, 2017). Environmental distances were represented as Euclidean distances of a scaled and centered environmental data matrix. The correlation between geographical and environmental distances was analyzed using the Mantel test in the R package ‘ade4’ (Dray & Dufour, 2007).

We fitted multivariate GLM models to analyze the effect of environmental variables on seabird species composition (i.e., spatial turnover). The multivariate GLM was based on a binomial distribution and included the same predictor variables used in the global model for univariate analyses. The analysis was performed using the function manyglm of the R package ‘mvabund’ (Wang, Naumann, Wright, & Warton, 2012). The significance of the model terms was assessed using analysis of deviance considering α = 0.01. Since ordinations are a good way to represent multivariate data graphically, we performed a canonical analysis of principal coordinates (CAP, or distance-based redundancy analysis; Legendre & Andersson, 1999; Anderson & Willis, 2003) to depict changes in seabird assemblage composition across islands in relation to environmental variables. CAP is a model-based ordination technique that can represent changes in species composition given specific hypotheses (i.e., environmental variables). The analysis was based on a resemblance matrix of Jaccard distances of species occurrences and all of the environmental predictors. Although CAP does not necessarily represent the fitted values of multivariate GLMs, we observed that the function capscale in the R package ‘vegan’ (Oksanen et al., 2009) produced CAP ordinations that were consistent with our multivariate GLM results (data not shown). The ordination was plotted using the R package ‘ggplot2’ (Wickham, 2009).