Proceedings of the Royal Society B: Biological Sciences
You have accessResearch article

Biodiversity response to natural gradients of multiple stressors on continental margins

Erik A. Sperling

Erik A. Sperling

Integrative Oceanography Division, Scripps Institution of Oceanography, La Jolla, CA 92093-0218, USA

Google Scholar

Find this author on PubMed

Christina A. Frieder

Christina A. Frieder

Integrative Oceanography Division, Scripps Institution of Oceanography, La Jolla, CA 92093-0218, USA

Google Scholar

Find this author on PubMed

Lisa A. Levin

Lisa A. Levin

Integrative Oceanography Division, Scripps Institution of Oceanography, La Jolla, CA 92093-0218, USA

Center for Marine Biodiversity and Conservation, Scripps Institution of Oceanography, La Jolla, CA 92093-0218, USA

[email protected]

Google Scholar

Find this author on PubMed


    Sharp increases in atmospheric CO2 are resulting in ocean warming, acidification and deoxygenation that threaten marine organisms on continental margins and their ecological functions and resulting ecosystem services. The relative influence of these stressors on biodiversity remains unclear, as well as the threshold levels for change and when secondary stressors become important. One strategy to interpret adaptation potential and predict future faunal change is to examine ecological shifts along natural gradients in the modern ocean. Here, we assess the explanatory power of temperature, oxygen and the carbonate system for macrofaunal diversity and evenness along continental upwelling margins using variance partitioning techniques. Oxygen levels have the strongest explanatory capacity for variation in species diversity. Sharp drops in diversity are seen as O2 levels decline through the 0.5–0.15 ml l−1 (approx. 22–6 µM; approx. 21–5 matm) range, and as temperature increases through the 7–10°C range. pCO2 is the best explanatory variable in the Arabian Sea, but explains little of the variance in diversity in the eastern Pacific Ocean. By contrast, very little variation in evenness is explained by these three global change variables. The identification of sharp thresholds in ecological response are used here to predict areas of the seafloor where diversity is most at risk to future marine global change, noting that the existence of clear regional differences cautions against applying global thresholds.

    1. Introduction

    Continental margins cover approximately 11% of the ocean's seafloor, but are disproportionately important in providing goods and ecosystem services [1], including the biggest percentage-wise increase in fisheries catch in the last three decades [2]. Margins are also important to long-term carbon cycling as the largest sink for carbon produced on land and the shelf [3]. Their biodiversity is the source of these services, yet much of it remains unexplored [1]. Global analysis of deep-sea benthic communities have demonstrated that as biodiversity decreases, there can be concomitant declines in ecosystem functions including secondary production, transfer efficiency of detritus to higher food webs and organic matter recycling [4]—each with implications for fisheries and carbon-cycle processes. Benthic organisms on continental margins are currently facing ocean acidification, global warming and oxygen loss [57], ultimately resulting from increases in atmospheric CO2. Rates of change are unprecedented, raising questions about whether (and how) communities will adapt, and whether responses will reflect synergistic interactions among multiple stressors [810].

    Predicting how marine life may respond to future environmental change can be approached by (i) examining past changes through the fossil and/or sedimentary records [1113], (ii) conducting controlled experiments that allow for rigorous tests of specific mechanistic effects, and (iii) examining ecological shifts along naturally occurring environmental gradients [14,15]. This third approach can account for potential evolutionary adaptation and can examine responses to multiple stressors over a range of conditions [10].

    Many areas of the world's oceans, such as enclosed seas, upwelling zones and estuaries, have naturally or anthropogenically induced low oxygen and pH [16]. On upwelling margins sharp, persistent natural gradients in temperature, CO2 and oxygen allow study of adaptation to these stressors [1]. The upper and lower boundaries of oxygen minimum zones (OMZs) support strong zonation of communities, with rapid shifts in faunal diversity that are clearly responding to hydrographic changes, often in a threshold-like manner [17,18]. In many instances, the resident fauna have had long periods of time (millions of years [19]) to adapt to the extreme conditions found within OMZs, although they can be dynamic on glacial/interglacial (10 000 years) [20] or much shorter [12] time scales. In this study, we exploit naturally occurring hydrographic gradients across upwelling margins to examine the potential influence of environmental variables on the community structure of continental margin macrobenthos (polychaetes, crustaceans, molluscs, echinoderms and other invertebrates). Data were assembled for near-bottom hydrographic variables and sediment-dwelling macrofaunal organisms (those retained on a 0.3 mm sieve; a fundamental component of the food web and the size class for which data are most abundant [17]) from 94 margin stations in the eastern Pacific Ocean and Arabian Sea at depths ranging from 100 to 3400 m (median 800 m). Key to using natural hydrographic gradients to understand the relative predictive power of different variables is the general lack of correlation among environmental variables (including depth) across different oceanographic domains. On a local scale, there can be clear covariation between variables; for example, oxygen and carbon dioxide covary due to aerobic respiration in a given water mass (e.g. [2123]). These correlations break down, though, when comparisons are made across water masses—for instance, between North Pacific subarctic and subtropical gyres [24], or between the northern and southern California Current System [21]. Therefore, relatively little covariation in environmental parameters exists among sites that span multiple water masses and ocean basins, as in this study (electronic supplementary material, figures S2–S4).

    Here, we employ variance partitioning tools to consider the following three questions. (i) What is the relative explanatory power of temperature, oxygen and the carbonate system for diversity and evenness? (ii) Do threshold-type effects exist and at what range of values? (iii) At what level of one predictor does a secondary predictor become important? These latter questions are particularly relevant as many ecological responses appear to exhibit threshold-type effects at very low oxygen levels [2527]. These questions are addressed with regression trees and random forests, two non-parametric variance partitioning techniques well suited to determining the relative contribution of multiple predictor variables in situations where responses are nonlinear [28]. While the environmental parameters investigated have physiological significance, these analyses do not, of course, indicate that a given environmental parameter is specifically (mechanistically) driving diversity dynamics. Nonetheless, such analyses may usefully identify regions where communities are at risk of substantial losses of diversity or evenness due to relatively small-magnitude environmental changes.

    2. Material and methods

    (a) Dataset assembly

    Paired environmental and macrofaunal data were collected from published studies on upwelling margins from 78 sampling stations, and combined with newly obtained environmental and biological data from 16 stations on the Pacific margins of Costa Rica, Mexico and the USA (electronic supplementary material, table S1 and dataset 1). We focused on upwelling margins as these habitats are exposed to the sharpest environmental gradients and because macrofaunal studies from non-upwelling regions (e.g. North Atlantic) often do not publish associated environmental data. All studies used a sieve size of 0.3 or 0.5 mm, a size that captures the majority of macrofaunal diversity [29] (77 stations used a 0.3 mm sieve and 17 stations used a 0.5 mm sieve). We note that the area of seafloor sampled at each station differs between studies and even within studies (electronic supplementary material, dataset 1), which could differentially affect the accuracy of reported diversity metrics [30].

    Where they are available, we have considered only in situ environmental data rather than using gridded climatological data such as from the World Ocean Atlas (WOA). While such climatological data would allow for the standard deviation of a given predictor variable across time to be included in the analysis, in situ data are preferred because (i) the WOA is not especially precise with respect to oxygen as it uses 1° grid cells; (ii) WOA averages data over multiple years during which conditions are known to have changed, and macrofaunal organisms are known to respond to environmental change on timescales of less than 1 year [15,31]; and (iii) the included studies were published between 1991 and 2015, during which time there has been some global change with respect to oxygen, temperature and pCO2, and thus the in situ values are most relevant to the observed biological response. We only considered stations at water depths greater than 100 m, and consequently beneath the majority of wave mixing as well as the calculated winter mixed-layer depth for all localities [32]. Further, the majority of our sites are considerably deeper (median: 800 m) and thus expected to experience less variability [15,31]. For instance, in a seasonal study of the Pakistan margin [31], the 140 m station experienced relatively large changes to oxygen, temperature and salinity, whereas the 300 m station and below experienced variation that ranged from less than instrumental error to approximately 10%. CTD casts from that margin indicate that inter-seasonal differences in environmental parameters converge at more than 200 m depth [15], which is common on most margins (e.g. [33]). In our dataset, only approximately 7% of stations are above 200 m depth, and thus the impact of seasonal variability is expected to be limited.

    In general, most environmental data in published studies (including temperature, oxygen and per cent organic carbon) were collected at the time of benthic faunal sampling. The exceptions to this are nine temperature measurements that were calculated from climatological data (noted in electronic supplementary material, dataset 1) and most carbonate system parameters. In the past, it has not been routine to measure carbonate system parameters during benthic faunal field programmes, and thus this predictor variable must be estimated from global databases. Direct measurements were used where available (electronic supplementary material, dataset 1), but for the majority of the stations carbonate system parameters were obtained from the dataset of Goyet et al. [32]. This dataset contains gridded climatological fields for total dissolved inorganic carbon (DIC) and total alkalinity (TA) below the mixed layer as a 1° × 1° × 32 vertical-layer grid. For the Santa Barbara Basin, we used CalCOFI DIC and TA measurements taken during July 2009 at station 081.8046.9. pH, pCO2 and saturation state (with respect to calcium carbonate polymorphs) were calculated at in situ temperature, salinity and pressure using the seacarb package in R. Direct comparison between climatological [32] and measured [34] carbonate chemistry values at the same site suggests the level of difference (less than 1%) at bathyal depths is negligible compared with the wide range of values in the dataset (electronic supplementary material, figure S1).

    The degree of covariation between a given set of environmental parameters was investigated by using bivariate plots and calculating the coefficient of determination.

    (b) Analytical procedures

    Diversity metrics (H′(log2)) and J′) were calculated for each station using the software PRIMER. For J′ analyses, three stations from the Pakistan margin [15] with only one species present (J′ = 0/0, undefined) were removed. To maintain consistency, diversity metrics were re-calculated from published datasets or from datasets provided by the authors (all datasets used in this study are compiled in the electronic supplementary material, dataset 1).

    Regression trees recursively partition datasets based on an exhaustive search across all variables for the split that maximizes homogeneity of the two resulting subsets, with the model performance described by R2-values. Random forests build large numbers of bootstrapped regression trees using random selections of predictor variables at each node. The data not in the bootstrapped sample (‘out-of-bag’ data) are used to estimate the percentage of variation explained (pseudo-R2-values), with a negative or near-zero value indicating a very poor fit. Relative importance (explanatory capacity) of predictor variables is indicated by per cent change in mean square error when a predictor variable is randomly permuted. Partial dependence plots show the marginal effect of a given predictor variable on the response variable across the full range of values. We note that the predictor variables used here take different forms (e.g. linear and log units). While this would be an issue for many statistical analyses, it is not a problem for regression trees and random forests, as these analyses conduct exhaustive searches for the split that maximizes homogeneity in the two resulting subgroups and thus can accommodate all forms of categorical and continuous data [28,35].

    Regression tree analyses proceeded using the software package JMP. Regression trees can potentially be overfit, and thus we conservatively stopped splitting after four splits. To test whether this is still oversplit, we used the automatic stopping rule in JMP, which continues splitting until the fivefold cross validation R2-value is better than what the next 10 splits would obtain. For no analysis was four splits found to be oversplit. Random forest analyses proceeded using the package randomForest [36] in R using ntree = 1000 and mtry = 2. While the focus of this paper is on three major environmental parameters (oxygen, temperature and pCO2) whose changes are often hypothesized to provide the most pressing threats to marine life [6,37], a further random forest analysis for diversity and evenness was run using a broader range of possible predictor variables that included these three global change variables, all other parameters of the carbonate system, salinity, seafloor depth, mixed-layer depth, sieve size and area of seafloor sampled.

    3. Results and discussion

    (a) Correlation of environmental variables

    Tight correlation between pairs of environmental parameters would inhibit the ability of variance partitioning techniques to determine the relative predictive capacity of variables. Of the three global change variables, temperature is most strongly correlated with depth (R2 = 0.42), especially within ocean basins (eastern Pacific, R2 = 0.61; Arabian Sea, R2 = 0.77; electronic supplementary material, figure S2; all environmental correlations with depth had p < 0.05). Oxygen and pCO2 are moderately to poorly correlated with depth in the full dataset (R2 = 0.24 and 0.18, respectively) and eastern Pacific (R2 = 0.17 and 0.09, respectively; electronic supplementary material, figures S2 and S3). These two variables are more strongly correlated with depth in the Arabian Sea (R2 = 0.55 and 0.47, respectively; electronic supplementary material, figure S4), consistent with the relatively smaller geographical area and lower number of distinct water masses sampled. Overall, sufficient variation exists to distinguish the influence of depth from environmental variables, and depth was only ranked 8th of 13 in predictive capacity for diversity when all possible predictor variables were considered (electronic supplementary material, figure S6). Depth was ranked first in predictive capacity for evenness (electronic supplementary material, figure S7) but the complete model only explained a small percentage of variance (pseudo-R2 = 0.06), suggesting that the effect size of all investigated variables for evenness is small. In summary, while depth is recognized as an important correlate of diversity on non-upwelling margins and across large depth ranges [30], environmental factors appear to have more explanatory power across the strong hydrographic gradients present at upper continental slope depths on upwelling margins.

    With respect to pairwise correlations between the three global change variables, most variables were moderately to poorly correlated (e.g. R2 < 0.30), including all comparisons in the full dataset (electronic supplementary material, figure S2)—not surprising given the number of water masses considered. All pairwise correlations except for temperature and pCO2 in the full dataset and eastern Pacific, and temperature and oxygen in the Arabian Sea, had p < 0.05. The only highly correlated environmental variables were oxygen and pCO2 at the ocean basin level, specifically R2 = 0.57 for the eastern Pacific and R2 = 0.43 for the Arabian Sea. This relationship is to be expected given the biogeochemical linkage between these variables discussed above. While this suggests some caution should be applied to interpreting the differential effects of oxygen and pCO2, considerable variation does remain, thus allowing the individual effects to be discerned during variance partitioning.

    (b) Diversity: full dataset

    Machine learning methods for variance partitioning (regression trees and random forests) were first applied to better understand the potential for oxygen, temperature and pCO2 to explain variation in diversity, which was estimated with the Shannon diversity index (H′(log2)). In regression tree analysis, all three variables (T, O2 and pCO2) explain much of the variance in diversity (R2 on first three splits = 0.59; electronic supplementary material, figure S8). Oxygen concentration is the predictor variable that best explains variance, with diversity in communities living below 0.16 ml l−1 O2 (7 µM) around half that at higher oxygen concentrations (table 1; average H′ of 2.1 versus 4.0 in low- and high-O2 partitions). This oxygen level is similar—within available precision—to macrofaunal oxygen thresholds identified in the Pleistocene fossil record of the Santa Barbara Basin [12]. Further, this is the oxygen level identified on the Pakistan margin below which foraminifera-dominated processing of particulate organic matter occurs (relative to macrofauna), pointing to potential direct or indirect functional links between diversity and organic matter cycling [26]. It is emphasized though that the specific causal link between oxygen and diversity is unclear; it could be related to a direct physiological control [38,39], feeding efficiency [26], predator–prey dynamics [27], or a combination of these and other factors. Above 0.16 ml l−1 O2, the next recursive split in the regression tree is for temperature, with stations less than 7.1°C having higher diversity than those at higher temperatures (full regression tree in the electronic supplementary material, figure S8). These results are consistent with hypotheses suggesting that temperature may influence biodiversity differentially at low, medium and high temperatures [40].

    Table 1.Regression tree analyses of macrofaunal diversity (H′) and evenness (J′) with respect to temperature, oxygen and pCO2. Analyses were run for the full dataset, the eastern Pacific Ocean and Arabian Sea, and for the full dataset with sedimentary organic carbon (TOC) included as a fourth predictor variable. For each analysis, the predictor value for the first split is shown, as well as the average and standard deviation of the response variable for the subsets above and below that value. The R2-value represents the value for the first split, followed by the predictor variable for the second split. Graphical representations of the first four splits in regression trees and accompanying R2-values and cross-validated R2-values are shown in the electronic supplementary material, figure S8.

    response dataset stations first split predictor value response value (above; below) R2 first split second split
    diversity (H′) full 94 O2 0.16 ml l−1 4.0 ± 0.9; 2.1 ± 1.2 0.44 T
    diversity (H′) eastern Pacific 57 O2 0.098 ml l−1 4.0 ± 0.9; 1.9 ± 0.8 0.45 T
    diversity (H′) Arabian Sea 37 pCO2 1026 µatm 1.1 ± 0.9; 3.3 ± 1.1 0.56 O2
    diversity (H′) full with TOC 77 O2 0.17 ml l−1 4.0 ± 1.0; 2.1 ± 1.2 0.42 O2
    evenness (J′) full 91 T 9.8°C 0.68 ± 0.21; 0.82 ± 0.10 0.17 pCO2
    evenness (J′) eastern Pacific 57 T 10.2°C 0.61 ± 0.23; 0.81 ± 0.10 0.23 T
    evenness (J′) Arabian Sea 34 T 9°C 0.74 ± 0.16; 0.89 ± 0.07 0.29 pCO2
    evenness (J′) full with TOC 74 T 9.8°C 0.67 ± 0.21; 0.81 ± 0.10 0.18 pCO2

    These three variables also explain much of the variance in diversity in random forest analyses (pseudo-R2 = 0.49; figure 1) Oxygen concentration is again the predictor variable that best explains variance in the random forest (figure 1). The partial dependence plot from the random forest (figure 1) exhibits a sharp decrease through the same oxygen range indicated by regression trees, and suggests that above approximately 0.5 ml l−1 (22 µM) O2, oxygen plays little role in shaping diversity. pO2 (matm) is a measure that accounts for effects of pressure, temperature and salinity on oxygen availability to organisms, and when pO2 (matm) was used instead of [O2] (ml l−1), a nearly identical result was obtained with a regression tree first split at 6.6 matm and little relationship between of pO2 and diversity above approximately 20 matm (electronic supplementary material, figures S10–S11).

    Figure 1.

    Figure 1. Random forest analyses of macrofaunal diversity (H′) with respect to oxygen, temperature and pCO2. Predictive capacity (variable importance) analyses (a) measure the percentage change in mean square error when a predictor variable is randomly permuted. Partial dependence plots (bd) show the marginal effect of a given predictor variable on the response variable. Hash marks on the x-axis of partial dependence plots indicate the distribution of data (as deciles) across the variable range; areas with low data coverage may be less accurate.

    (c) Diversity: eastern Pacific and Arabian Sea

    The relative predictive power of temperature, pCO2 and O2 for explaining variance in diversity is different between ocean basins and regionally within the Pacific (figure 1; electronic supplementary material, figure S16). While the explanatory power of oxygen for diversity is high in both ocean basins, the sharpest drop in H′ in partial dependence plots is at approximately 0.17 ml l−1 (7.5 µM) O2 in the eastern Pacific and at approximately 0.5 ml l−1 (22 µM) in the Arabian Sea. Further, pCO2 emerges as the variable with the most predictive power in the Arabian Sea—with threshold behaviour at approximately 900 µatm in partial dependence plots and a first split at approximately 1000 µatm in the regression tree. By striking contrast, pCO2 appears less relevant to explaining variation in continental margin macrofaunal diversity in the eastern Pacific Ocean. Conversely, temperature has low predictive ability for diversity in the Arabian Sea, but has high predictive ability in the Pacific (figure 1).

    To investigate the effect of the carbonate system further, a separate random forest analysis was run using only pH, pCO2 and saturation states of calcium carbonate (Ω; aragonite and calcite). While carbonate chemistry parameters can covary on a local level, there is considerable variation based on differences in temperature, depth and salinity (electronic supplementary material, figure S5), allowing the effect of individual parameters to be partitioned. Results indicate that in the Arabian Sea pCO2 (values from 444 to 1140 µatm), and secondarily pH (values from 7.60 to 7.88), best explain variance in diversity (electronic supplementary material, figure S15). This suggests the possible driver of low extant diversity may be the physiological effects of increased pCO2 and decreased pH on organisms (through acid-base regulation from disturbances in extracellular pH [41]) rather than increased energetic costs of producing and maintaining a carbonate skeleton at undersaturation (note also both oceans are dominated by unmineralized polychaetes; electronic supplementary material, table S2). In the eastern Pacific, pCO2 levels extend far higher, and undersaturation is more extensive, yet the carbonate system parameters explain much less of the variance in diversity than in the Arabian Sea (figure 1; electronic supplementary material, figures S8 and S15).

    Why, then, does pCO2 appear to have high explanatory power for macrofaunal diversity in the Arabian Sea but not the eastern Pacific Ocean? The major taxonomic composition of the organisms is broadly similar (electronic supplementary material table S2), suggesting wholesale faunal replacement is not the cause. In the Pacific, long-term exposure of organisms to naturally high-CO2 waters [22,42] may lead to adaptations to high-CO2/low-pH conditions. By contrast, the Arabian Sea is characterized by much lower background pCO2 levels but higher temperatures. Here, the interaction of elevated temperature, which increases metabolic rate, could reduce tolerance to higher CO2 through energetic constraints (achieving appropriate aerobic scope) [41,43,44].

    (d) Evenness

    We further examined the extent to which evenness, as estimated by Pielou's J′ metric, is explained by these three global change variables. The regression tree indicates a first temperature split at approximately 10°C, with R2 = 0.17 (table 1); clearly, evenness is less well explained by these variables than diversity. Random forest analysis returned a negative pseudo-R2-value, although the analysis also finds temperature to be the predictor variable that best explains variance in evenness, and partial dependence plots also show a sharp drop in evenness through the 8–10°C range (figure 2). In permutation tests, oxygen performs no better than random. Thus, while temperature may have a small impact on evenness, other factors such as food supply may be equally important [25]. Notably, random forest analyses on the 74 stations with data available for sedimentary total organic carbon (a proxy for food) resulted in a positive, albeit small, pseudo-R2-value of 0.07 (with temperature remaining the variable with the best predictive capacity; electronic supplementary material, figure S14).

    Figure 2.

    Figure 2. Random forest analyses of macrofaunal evenness (J′) with respect to oxygen, temperature, and pCO2. Description of random forest plots as in figure 1. Eastern Pacific and Arabian Sea random forest shown in the electronic supplementary material, figure S9.

    (e) Predicting future ecological change

    Based on analysis of macrobenthos response to environmental gradients on modern continental upwelling margins, we predict disparate effects on ecosystem structure due to future global change. Diversity is likely to drop in response to projected decreases in oxygen levels and increases in temperature and pCO2, with variation between ocean basins depending on current adaptations and interactions of stressors. On the other hand, evenness may be relatively less affected by these investigated variables; continued research is needed to understand possible controls on evenness. An important message arising from these results is that while absolute thresholds can provide a guide, specific laboratory or empirical thresholds derived from one ocean basin may translate poorly to another. In particular, oxygen thresholds determined primarily from laboratory experiments on Atlantic species [45] are unlikely to be applicable to eastern Pacific species [46]; the same may also hold true for temperature and pCO2. Finally, the comparison of the eastern Pacific and Arabian Sea suggests that continental margin organisms may be well adapted to naturally occurring environmental stressors, but more sensitive to changes in stressors to which they are not normally exposed. If so, we predict that the biodiversity of relatively well-oxygenated and low-pCO2 North Atlantic margins (where the slight OMZ is characterized by pCO2 levels of approx. 500 µatm compared to approx. 1200 µatm in the North Pacific) will be strongly affected by deoxygenation and ocean acidification.

    Considering the basin-scale thresholds identified in this study, and making the assumption of niche stability [47,48], it is possible to predict areas of the seafloor most at risk of sharp biodiversity loss and attendant changes to carbon-cycle processes [4,26] by knowing stressor trajectories and identifying areas where present conditions are close to threshold values (figure 3). In general, upper to mid-slope depths (200–1500 m) are most at risk; depths more than 1500 m are not near threshold values of oxygen, temperature or pCO2. In the eastern Pacific, areas vulnerable to synergistic temperature and oxygen change are predicted for the southern tip of Baja to around San Francisco along the California coast, where critical oxygen and temperature levels co-occur in the approximately 400–800 m range depending on latitude, and on the Peru margin at approximately 600–800 m (figure 3a). In the Arabian Sea, large areas of the upper to mid-slope between Somalia and Iran are predicted to be vulnerable to increases in pCO2, with many areas also at risk from declining oxygen. The outer shelf on the Pakistan and Indian margins are also at considerable risk from these combined stressors, as are upper to mid-slope depths on the southwestern Indian margin (figure 3b). Diversity is of course only one factor influencing ecosystem function. Other responses such as identity and density of individual species, lifestyles, feeding modes, biomass or calcification may change at different levels than these identified diversity thresholds, with consequences for valued sediment services.

    Figure 3.

    Figure 3. Continental margin depths most likely to experience macrofaunal biodiversity loss due to deoxygenation, increasing temperature and increasing pCO2 for (a) eastern Pacific Ocean and (b) Arabian Sea. Depths along transects that are currently beyond the environmental threshold for the predictor variable with the most explanatory power are shown in black (natural low diversity). Depths along the transect that are close to diversity thresholds for the predictor variable with the most explanatory power—and most at risk of rapid diversity loss due to global change—are depicted in blue (eastern Pacific) and green (Arabian Sea). Depths close to diversity thresholds for the predictor variable with the second best predictive power are depicted with cross-hatching. Note that while the regression tree and random forest indicate decreases in diversity with respect to temperature at approximately 7°C in the eastern Pacific Ocean, the partial dependence plot indicates these effects are linear over a range of temperatures, and there may not be a specific threshold. Plotted diversity thresholds are based on figure resolution in Ocean Data View and the results of regression tree and random forest partial dependence plots. (Online version in colour.)

    4. Conclusion

    These analyses highlight the utility of modern ocean environmental gradients in providing a natural laboratory to study future oceans. Specifically, oxygen is identified as the variable that best explains variance in macrofaunal diversity, and suggests that continued research and forecasting efforts into the multiple (and complex) drivers of oxygen dynamics on upwelling margins [4951] and the resulting ecophysiological effects [39] should be a priority. This natural gradient approach accounts for adaptive responses arising from the evolutionary history of organisms, and illuminates responses over a complete range of variable space for multiple stressors in a manner not tractable in laboratory experiments. These results can also be used to design more relevant laboratory experiments, allowing a focus on threshold and vulnerable conditions. Together, the study of modern environmental gradients, the fossil record and targeted mechanistic experiments yield an integrated approach to predicting future responses to global change.

    Data accessibility

    The datasets supporting this article have been uploaded as part of the supplementary material and submitted to the Dryad Data Repository:

    Authors' contributions

    L.A.L. conceived the project, assembled the initial dataset, and contributed new macrofaunal data from the Costa Rica, Mexico, California and Oregon margins. E.A.S. and C.A.F. assembled the full dataset. C.A.F. calculated carbonate system parameters. E.A.S. conducted the analyses and wrote the paper in collaboration with L.A.L. and C.A.F. All authors gave final approval for publication.

    Competing interests

    We declare that we have no competing interests.


    E.A.S. and L.A.L. are supported by NSF-EAR 1324095.


    We thank B. Ingole, S. Sautya, D. Hughes, E. Vetter and B. Grupe for contributing raw data tables from published papers. We thank C. Rochman, N. Wegner, A. Beaudreau, G. Mendoza, M. Corrales, S. Choi and J. Gonzalez for assistance collecting and sorting samples appearing in our dataset as de novo data. We thank J. Bernhard and W. Berelson for assistance collecting new macrofaunal data, and A. Cau for assistance with initial data gathering. We thank Y. Takeshita, A. Knoll and V. Tunnicliffe for helpful discussion and critique, and M. Yasuhara, C. McClain and an anonymous reviewer for insightful comments during review.


    Published by the Royal Society. All rights reserved.