The role of island physiography and oceanographic factors in shaping species richness and turnover of nesting seabird assemblages on islands across the southeastern Pacific

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 (km2), island elevation (m), distance to the mainland (km), human density (individuals per km2), 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 log10(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 (log10 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).

研究区域
本研究的研究区域经度范围为西经130°44′至西经70°31′，纬度范围为北纬1°40′至南纬38°22′。本数据集涵盖6个群岛共计42个岛屿：皮特凯恩群岛（Pitcairn，4个岛）、拉帕努伊岛（Rapa Nui，4个）、德斯文图拉达斯群岛（Desventuradas，3个）、胡安·费尔南德斯群岛（Juan Fernández，3个）、加拉帕戈斯群岛（Galápagos，13个）以及智利沿海岛屿（15个）。本数据集包含的岛屿与其他岛屿及大陆的隔离程度各不相同，分布于不同的生物地理分区（Spalding等，2007；图1）。除皮特凯恩群岛位于太平洋板块（Pacific Plate）外，绝大多数大洋岛均位于纳斯卡板块（Nazca Tectonic Plate）。这些岛屿均属于火山岛及海山链。

环境变量
我们从文献及在线数据库中整理了岛屿特征信息。在去除高度共线性变量（即皮尔逊相关系数（Pearson correlation coefficient）>0.8）后，我们选取7个环境变量作为分析的预测因子：岛屿面积（km²）、岛屿海拔（m）、距大陆距离（km）、人口密度（individuals per km²）、海表盐度（PSS）、海表温度（℃）以及初级生产力（g·m⁻³·day⁻¹）（附录S1）。岛屿面积、海拔及距大陆距离来自文献，或通过R语言（R Core Team, 2016）及Google Earth（www.google.com/earth）中的工具估算。人口密度数据来自各国统计及普查机构提供的公共数据库（详见附录S1中的元数据）。所有海洋学变量均从Bio-ORACLE在线数据库的栅格地图中提取（Assis等，2018；Tyberghein等，2012），这些地图在赤道处的栅格分辨率约为9.2 km。我们以每个岛屿的中心为基准，创建50 km半径的缓冲区，使用R包‘raster’（Hijmans, 2017）计算每个海洋学变量的平均值。

海鸟群落
我们从文献、在线数据库及未发表报告中整理了东南太平洋沿海与大洋岛屿上53种海鸟的筑巢分布数据（附录S2）。所参考的文献大多发表于近30年，同时也纳入了早期出版物以佐证部分岛屿上已确认的海鸟记录（详见附录S2中的元数据）。对于多数智利岛屿，我们使用了1999-2018年（19年间）在智利沿海岛屿、2013-2018年（5年间）在智利大洋岛屿（即德斯文图拉达斯、拉帕努伊及胡安·费尔南德斯群岛）开展的多次科考的实地记录。所有海鸟数据均整理为以岛屿为行、物种为列的二元矩阵（1=存在，0=缺失）（附录S2）。

数据分析
海鸟α多样性——假设1与假设2
我们通过物种丰富度（即每个岛屿上筑巢海鸟物种的总数）来衡量各岛屿的海鸟α多样性。分析前，我们使用Shapiro-Wilk检验（Shapiro-Wilk test）检查所有环境变量（岛屿特征与海洋学变量）的正态性。对于不符合正态分布的环境变量，分别通过平方根变换或log10(x+1)变换进行校正。我们基于泊松分布与对数连接函数构建广义线性模型（Generalized Linear Model, GLM），以分析海鸟物种丰富度与环境变量之间的关系。模型选择与模型平均使用R包‘MuMIn’（Barton, 2018）完成。使用dredge函数基于最低二阶Akaike信息准则（Akaike Information Criterion, AICc）进行模型筛选，共生成128个不同的模型子集。使用model.avg函数进行全模型平均，纳入AIC权重（AICw）>0.05的最优模型子集。我们绘制了与每个环境变量相关的偏残差图，以展示其对物种丰富度的影响。最终模型纳入的环境变量包括：经log10变换的岛屿面积、经平方根变换的岛屿海拔、经平方根变换的距大陆距离、经平方根变换的人口密度、盐度、表层初级生产力及海表温度。最后，我们使用R包‘ape’（Paradis, Claude, & Strimmer, 2004）与‘ncf’（Bjornstad & Cai, 2018）中的残差距离相关图与莫兰I检验（Moran’s I test）程序，评估GLM结果的空间独立性。

海鸟β多样性——假设3与假设4
我们使用R包‘betapart’（Baselga & Orme, 2012）中的程序计算分区β多样性，该方法将β多样性拆分为周转组分与嵌套组分。β周转组分反映物种替代现象，而β嵌套组分反映物种丰富度的差异。这两个组分基于用于描述β多样性的相异性指数（常称为β总多样性）计算得出。我们使用各岛屿物种分布的Jaccard相异性指数（Jaccard dissimilarity index）计算β多样性。我们基于三次多项式拟合的线性模型，分析β多样性与地理距离、环境距离之间的关系。地理距离使用R包‘geosphere’（Hijmans, 2017）中的distm程序计算。环境距离以标准化（缩放与中心化）后的环境数据矩阵的欧氏距离（Euclidean distance）表示。我们使用R包‘ade4’（Dray & Dufour, 2007）中的Mantel检验（Mantel test）分析地理距离与环境距离之间的相关性。我们构建多元广义线性模型（multivariate Generalized Linear Model, multivariate GLM），以分析环境变量对海鸟物种组成（即空间周转）的影响。该多元GLM基于二项分布，纳入了单变量分析中全局模型所用的全部预测因子。分析使用R包‘mvabund’（Wang, Naumann, Wright, & Warton, 2012）中的manyglm函数完成。模型项的显著性通过偏差分析评估，显著性水平α=0.01。由于排序是多元数据可视化的有效方式，我们进行了主坐标典范分析（Canonical Analysis of Principal Coordinates, CAP，即基于距离的冗余分析），以展示不同岛屿上海鸟群落组成随环境变量的变化。CAP是一种基于模型的排序技术，可基于特定假设（即环境变量）展示物种组成的变化。该分析基于物种分布的Jaccard距离相似性矩阵与所有环境预测因子。尽管CAP未必对应多元GLM的拟合值，但我们发现R包‘vegan’（Oksanen等, 2009）中的capscale函数生成的CAP排序结果与我们的多元GLM结果一致（未展示数据）。排序图使用R包‘ggplot2’（Wickham, 2009）绘制。