Mutualisms weaken the latitudinal diversity gradient among oceanic islands – Nature

Datasets

Regional plant distribution data along with explanatory variables were extracted from the GIFT v3.0 database³ using the GIFT R package⁴². The GIFT database integrates plant distributions from floras and regional checklists with geographic, environmental and socioeconomic data; it includes data for nearly 3,400 regions, over 350,000 plant species, and about 4 million species-by-region occurrences. This database is unique in that it includes many island regions in addition to mainland regions. Here we used all GIFT regions for which checklists of native angiosperms were available. We removed islands for which island geology (that is, volcanic, floor, shelf, fragment and others) was undetermined, and condensed island types into two major types: oceanic (including atoll, floor or volcanic) or non-oceanic (shelf or fragment). Oceanic islands represent newly formed land masses, assumed to be colonized de novo, whereas non-oceanic island have some historical connection to mainlands, and therefore direct contact with these source populations. Explanatory variables included absolute latitude and longitude of the region’s centroid, area (in km²; log-transformed for all analyses), and elevational range⁴³ (difference between lowest and highest elevation in m above sea level). When elevation range was unknown or reported as zero from aerial elevation maps, we assigned an elevation of 1 m as a minimum necessary elevation for data collection (that is, the island must be above sea level to be detected). For islands, we also included island isolation, or distance to the nearest mainland⁴⁴ (km). Prior to analyses, we removed islands smaller than 6 km², to remove any small island effect^45,46,47 (that is, small islands exhibiting a weaker species–area relationship than larger islands).

We considered three mutualisms—pollination syndrome, mycorrhizal status and N-fixing status—as they are known to influence plant establishment and fitness^48,49,50 and have the most complete data. Determination of mutualist status for pollination, mycorrhizal, and N-fixing mutualisms followed similar approaches. Pollination syndrome status for each plant species was extracted from the GIFT database as either biotic or abiotic pollination. Pollination syndrome was assigned to a given species first by matching species-level data^51,52 (using https://ecotraits.landcareresearch.co.nz/, http://tropical.theferns.info and https://www.floraweb.de); if not possible, we relied on family level assignment^53,54. Mycorrhizal status was assigned using the FungalRoot database³⁸. We assigned mycorrhizal status based on species-level status when possible. However, if species-level assignment was not possible, we applied genus-level assignment to the given species. We included four major mycorrhizal types: AM, EM, orchid mycorrhizal (ORC) and NM plant species. We assigned ambiguous (AMNM) or dual mycorrhizal types (AMEM)—those plant species that have been found as exhibiting either (AMNM or AMEM) or both mycorrhizal statuses (AMEM)—to AM as they are capable of forming these mycorrhizal types and may therefore experience the AM filter¹⁹. Our assignment of pollination status and mycorrhizal status to unknown species based on family and genus level, respectively, is valid as these mutualism types have been shown to be phylogenetically conserved at these taxonomic scales^23,55. Finally, N-fixing status was assigned using the most recent N-fixing plant database from Werner et al.⁵⁶. We first assigned N-fixing status (either N-fixing or non-N-fixing) at the species level; if this was not possible, we assigned a family level proportion to the remaining species in a given assemblage²². In total, we assigned 256,499 out of 880,176, 388,505 out of 1,816,146, and 7,645 out of 318,926 species to species-level assignment for pollination, mycorrhizal and N-fixing status, respectively. Note that the total number of species vary depending on whether information for genus or family was available for each mutualism status; we removed unmatched species for which family- or genus-level mutualist status was unknown.

In order to address how differences in resolution of assignments might affect our main conclusions, we conducted a sensitivity analysis. Specifically, for each mutualism type, we relabelled 20% of the mutualist-associated plant species as non-mutualists, with 20% of each flora relabelled from biotically pollinated to abiotically pollinated, N-fixing to non-N-fixing, or arbuscular mycorrhizal to non-mycorrhizal. We then reran models testing for island species deficit and contribution to species deficit (see ‘Statistical analyses’ for more details). Our major conclusions are robust to this sensitivity analysis, with the contribution to species deficit remaining highest for mutualism-associated plant species at low latitudes, and the relative importance of the mutualism filter emerging as a driver that is at least as important as abiotic drivers (Extended Data Fig. 5).

Statistical analyses

We tested for the plant LDG and evaluated whether this pattern was weaker on oceanic islands than mainlands (hypothesis (1)). We modelled species richness within each mutualism type (pollination syndrome, mycorrhizal and N-fixing) for each mutualist status (biotically or abiotically pollinated; AM, EM, ORC or NM; N-fixing or non-N-fixing) on mainlands based on absolute latitude. We then used these models to predict expected species richness for each mutualist status for each island based on mainland communities at the same latitude. We used these mainland-based predictions and observed island floras to calculate total species deficit, proportional species deficit, and contribution to species deficit for each mutualist status on each island. We then tested for drivers of the difference between island and mainland species richness at equal latitude (species deficit, extracted from mainland model), including both abiotic variables and an aggregate metric of proportion of mutualistic plant species on the mainland (mutualism filter strength, hypothesis (2)). Finally, we tested for the role of mutualist status within each mutualism type, latitude, and other biogeographical variables in predicting both proportional species deficit and contribution to species deficit (hypothesis (3)). All models included a residual autocovariate to correct for spatial autocorrelation⁵⁷, with neighbourhood distance set to 2,000 km and weighted by inverse distance. All analyses were done in R 3.4.1⁵⁸, with linear mixed models constructed with the lme4 package⁵⁹; significance for these models was tested using lmerTest⁶⁰, using Satterthwaite’s approximations for t-test and corresponding P values, with a P value of <0.05 used as the threshold for significance.

To model the LDG across our dataset, we used a generalized linear model (GLM) with negative binomial errors and a log-link function. The response variable was species richness (number of species); the fixed effect predictors were land type (mainland, oceanic island or non-oceanic island), absolute latitude, and an interaction between absolute latitude and land type. This analysis tested for a significant weakening of the LDG on each island type relative to the mainland (interaction between absolute latitude and land type in initial model). Results did not change qualitatively if we also included area in these models. However, we directly assessed the effect of island area in our subsequent deficit models; therefore, we did not include it here. We then reran models within each land type to examine the strength of the LDG in each land type separately.

We implemented two null approaches to confirm that the dampened oceanic island LDG cannot arise from a filter applied equally across all latitudes. We first generated an oceanic island null model to directly test using the model described above. To do this, we took individual mainland values of richness and multiplied them by the ratio of island to mainland species richness at the lowest latitude island (ratio = 0.12 at 0.02° latitude). This approach forces both intercepts of the oceanic island and null oceanic island to be equal; therefore, deviation from the null model can only be detected at higher latitudes. We found that oceanic islands exhibited a significantly weaker LDG than this null oceanic island model (Supplementary Table 1). Second, we examined the ratio of island to mainland species richness predicted from the LDG model across latitude. We find that the ratio increases with increasing latitude, more than doubling over 60° of latitude (Fig. 2a and Extended Data Fig. 2), confirming that the cumulative effect of all filters on island plant species richness increases with proximity to the Equator.

Our three metrics of island (1) species deficit, (2) proportional species deficit and (3) contribution to species deficit were determined by comparing observed island richness to the predicted species richness of an equal latitude mainland community. We used generalized additive models⁶¹ to model the relationship between latitude and species richness for total mainland species richness and for each mutualistic status within each mutualism type (biotically or abiotically pollinated; AM, EM, ORC, or NM; N-fixing or non-N-fixing). We predicted mainland species richness using a smoothed term of absolute latitude (with no limits on k) with a negative binomial error distribution and a log-link function. We used these models to estimate the expected species richness of a mainland community at the same absolute latitude of each island for each group. We then used observed and expected values of total species richness for each island community to calculate the island species deficit (expected species richness based on mainland generalized additive model minus observed total species richness), proportional species deficit (the proportion of expected species richness that was not observed on each island, within each of the eight mutualist statuses; for example, AM proportional deficit = (AM expected – AM observed)/AM expected), and the contribution to total species deficit (species deficit of each mutualist status divided by total species deficit of a given island; for example, AM contribution to species deficit = (AM expected – AM observed)/(total species richness expected – total species richness observed)). Species deficit (1) represents the raw difference between island and equal latitude mainland species richness, proportional species deficit (2) the proportion of mainland species lost from an equal latitude island within a particular mutualist status, and the contribution to species deficit (3) represents the relative contribution of a particular mutualist status to the total species deficit.

To model species deficit, we used a weighted GLM with Gaussian errors and an identity link (i.e. no transformation); sampling units were weighted by their sample size (expected species richness based on latitude). We included the fixed effects of absolute latitude, area, distance, elevation range, precipitation, and mutualism filter strength. The mutualism filter strength variable represented the proportion of mainland plant species that are mutualistic with at least one of the three mutualist types we assessed (biotic pollinators, AM fungi, or N-fixing bacteria). To assess relative importance of different variables, we subsequently used (1) AIC model averaging to determine the significance of each variable included in our model and (2) variance partitioning for the final model from model selection using the relaimpo package^62,63. Finally, we ran additional models to test interaction terms that we hypothesized may be important a priori. Specifically, we tested interactions between mutualism filter strength and each abiotic driver (area, distance, precipitation, and elevational range) as well as absolute latitude, followed by interactions between absolute latitude and each abiotic driver. We accomplished this by sequentially adding interactions, and we only retained terms that improved model fit, as determined by AIC values.

Our approach to estimating proportional and contribution to species deficit produced implausible estimates for a few islands with unusually high species richness (occurring in pollinator, mycorrhizal, and N-fixing analyses on 17, 41, and 2 out of 212 islands, respectively). Rather than removing these observations, we constrained these extremes to the smallest and largest values that should occur for proportional species deficit and contribution to species deficit (0 and 1; ‘constrained response’ models). We report results from this approach because it provided the best fit to the data. However, we also confirmed that this approach did not influence our results by rerunning the analyses including these extreme values (‘extreme response’ models) and using a linear model with Gaussian errors and the same predictors as we report in this manuscript. All results are reported in Supplementary Tables 3–5.

To model proportional species deficit, we used a weighted GLM with Gaussian errors and an identity link. A separate model was run for each mutualism type (pollinator syndrome, mycorrhizal and N-fixing). Fixed effects included absolute latitude, area, distance, elevation range, precipitation, mutualist status and the interaction between absolute latitude and mutualist status. Region was specified as a random effect to control for repeated measures of each island (that is, the separate measures of each level of mutualist status). These models allowed us to test for a consistent biotic mutualist filter within each mutualism type and evaluate whether this varied with latitude. To do this, we used post hoc contrasts using pairwise comparisons of model coefficients with least squares means.

To model contribution to species deficit, we used a weighted GLM with Gaussian errors and an identity link. We included the fixed effects of absolute latitude (as a polynomial function), mutualist status (either pollinator syndrome, mycorrhizal type or N-fixing type), area, distance, elevation range, and the interactions between mutualist status and absolute latitude. Region was again specified as a random effect. To confirm that these patterns were not influenced by interactions between latitude and other biogeographical variables, we separately modelled each mutualist status within each mutualist type with a more complex set of predictors. We created full models that added interactions between absolute latitude and area, distance, precipitation, and elevational range, and we used stepwise model reduction to simplify these models based on significance (P < 0.05) to obtain the most parsimonious model. For both proportional species deficit and contribution to species deficit, we confirmed that alternative modelling approaches (such as beta regression or binomial errors with a logit link) did not produce a better fit to the data.