Abundance and local-scale processes contribute to multi-phyla gradients in global marine diversity.

Among the most enduring ecological challenges is an integrated theory explaining the latitudinal biodiversity gradient, including discrepancies observed at different spatial scales. Analysis of Reef Life Survey data for 4127 marine species at 2406 coral and rocky sites worldwide confirms that the total ecoregion richness peaks in low latitudes, near +15°N and -15°S. However, although richness at survey sites is maximal near the equator for vertebrates, it peaks at high latitudes for large mobile invertebrates. Site richness for different groups is dependent on abundance, which is in turn correlated with temperature for fishes and nutrients for macroinvertebrates. We suggest that temperature-mediated fish predation and herbivory have constrained mobile macroinvertebrate diversity at the site scale across the tropics. Conversely, at the ecoregion scale, richness responds positively to coral reef area, highlighting potentially huge global biodiversity losses with coral decline. Improved conservation outcomes require management frameworks, informed by hierarchical monitoring, that cover differing site- and regional-scale processes across diverse taxa, including attention to invertebrate species, which appear disproportionately threatened by warming seas.

INTRODUCTION

Since the earliest voyages of discovery, the remarkable richness of tropical plants and animals has intrigued the greatest traveling scientists, including von Humboldt (), Darwin (), Wallace (, ), and Bates (). They reported an exceptionally high diversity of species at the regional scale, in terms of taxonomic collections, and also at the local scale, with their senses overwhelmed by the profusion of life observed in rainforests and on coral reefs. Since then, many studies have attempted to quantify historical observations and suggest causes of the “latitudinal gradient in biodiversity,” where species richness peaks in the tropics. Latitudinal gradients appear to be a global phenomenon for many groups of terrestrial and marine species but vary across taxa (, ). Patterns are scale-dependent (), with the latitudinal biodiversity gradients apparently stronger when reported as γ (regional) than α (local) diversity (). Although most studies assume that the underlying gradient should be unimodal and peak at the equator, a recent review consolidating many disparate surveys found that marine species richness was bimodal and asymmetric with latitude, dipping at the equator, and higher in the Northern Hemisphere (). However, subsequent analyses of 51,000 species that accounted for sampling effort showed that the latitudinal gradient was symmetric and bimodal () and that γ diversity was less sensitive to sampling bias than α diversity.

More than 30 different hypotheses have been proposed to account for global patterns of richness (, ), with approximately half favoring historical and evolutionary processes (), such as the accumulation of species in the tropics through rapid speciation rates or long-term environmental stability (). Other hypotheses invoke global variability in contemporary physical and ecological processes, including temperature and species interactions (Table 1) (, , ). Nevertheless, consensus on mechanisms that create and maintain global patterns of biodiversity is lacking. Among the many investigations on biogeographic gradients published to date, relatively few have extended beyond one continent or a single taxonomic group, and none among this group have accounted for the potentially confounding influence of sampling variation and species’ abundance at local scales, whereby richness is partly dependent on the number of observed individuals. A focus on a phylogenetically narrow set of taxa, particularly studies of single clades, risks confounding through clade-specific responses to environmental drivers that may not generalize across taxonomic groups (). Improved knowledge on how diversity is shaped by multiscale processes and variability in abundance is urgently needed for effective conservation because expanding human pressures are known to act at different scales of space and time.

Table 1

Hypotheses proposed to explain latitudinal biodiversity gradients.

Associated predictor metrics tested in this study and transformations applied in linear models are also listed. SST, sea surface temperature; Chl, chlorophyll; NO3, nitrate; PAR, photosynthetically active radiation; sqrt, square root.

Hypothesis	Metric	Measurement
TemperatureSpecies richness increases with temperature (76), perhaps dueto increased metabolic rates leading to faster generation timesand hence increased rates of mutation, evolution, and accumulationof species (30). Alternatively, niche conservatism may constrainthermal tolerances of most species to warm conditions (13, 31, 77),or temperature may indirectly link to diversity through its influenceon productivity.	SST (log)	Mean SST of sites, transformed to unitsof energy (−1/kT, where k is the Boltzmannconstant and T is in Kelvin) (30).
Primary productivityLight irradiance and nutrients drive primary productivity,which, at low to moderate levels at least, increases abundanceand species richness (32, 33).	Chl (log)NO3 (log)PAR (sqrt)	Satellite derived mean Chl, NO3, and PARmeasurements per 5–arc min (9.2 km) gridcells compiled in Bio-ORACLE (ocean rastersfor analysis of climate and environment) (71).
AreaSpecies richness increases with geographic area, independentof habitat diversity, in part through increased persistence ofanimals with large range requirements (34).	Shelf (sqrt)	Total continental shelf area (water depths,<200 m), from NOAA (National Oceanicand Atmospheric Administration) NGDC(National Geophysical Data Center) (78).
	Coral (log x + 1)	Total area of coral reef within ecoregion,assessed through “Reefs at Risk” (79).
Natural disturbance regimeSpecies richness decreases with, or peaks at intermediatelevels of, environmental disturbance.	Cyclone	Percentage of sites surveyed in an ecoregionwith a cyclone of at least category 3 passingwithin 50 km in the past 10 years, from IBTrACS(International Best Track Archive for ClimateStewardship) (80).
	SST range (log)	Difference between annual mean monthlymaximum and minimum SST, averagedacross years, from Bio-ORACLE (71).
Climatic stabilitySpecies richness is highest at locations with landscapeslittle affected by recent cataclysmic geological andclimatic events and thus have avoided major extinctions.	Stability (log)	Inverse of SD of mean water temperaturesince the last glacial maximum, estimatedfor three time periods: 0 ky (thousand years),8 ky, and 21 ky before present, fromPaleo-MARSPEC (Ocean Climate Layersfor Marine Spatial Ecology) (81).
Fragmentation and connectivitySpecies richness varies with rate of speciation, which isdependent on frequency of vicariant events, includingbreaks in population connectivity through isolation ofisland populations, coastal fragmentation, and changingdispersal pathways into adjoining regions.	Coast (log)	Total length of coastline (km) withinthe ecoregion.
	Connect	Number of abutting coastal ecoregions.
	Islands (log)	Number of separated coastline features onmap of the ecoregion.
Human disturbanceSpecies richness is negatively affected by human activities.	Pop (log x + 10)	Estimated population density, calculatedby fitting a smoothly tapered surface toeach settlement point on a year 2000world population density grid (82), usingthe quadratic kernel function describedin the study of Silverman (83). Populationswere screened for a density greater than1000 people per 0.04° cell, and thesearch radius was set at 3.959° (28).

The recent expansion of the Reef Life Survey (RLS) data set () to include the abundances of more than 4000 species across 11 classes of visually conspicuous mobile macroscopic organisms, in 44 countries (fig. S1), permits a rigorous evaluation of many of the long-standing hypotheses regarding the latitudinal biodiversity gradient based on a standardized methodology. A total of 13,858 transect blocks were censused during daylight at 2406 sites, generating records of 2516 actinopterygian fish (bony fishes), 66 chondrichthyan fish (sharks and rays), 8 reptile, 6 mammal, 105 holothurian (sea cucumbers), 104 echinoid (sea urchins), 180 asteroid (sea stars), 99 crinoid (feather stars), 805 gastropod, 28 cephalopod (octopi and cuttlefish), and 210 decapod and stomatopod crustacean (crabs, lobsters, and mantis shrimps) taxa. In total, the fauna considered represents <10% of all marine species () but >50% of large (>2.5 cm) species in regional checklists and field guides of shallow reef taxa [for example, the studies of Allen et al. (), Allen and Robertson (), Hickman (), Hickman and Finet (), Lamb and Hanby (), and Haussermann and Forsterra ()].

We calculated diversity at three scales for each of the 82 marine ecoregions of the world [of 233 worldwide ()]: by averaging observed species richness per transect within each site and then averaging the mean site richness values per ecoregion (mean observed species density per 50 m2 RLS transect, α1 diversity); by estimating species richness at the local scale (extrapolated asymptotic richness for RLS sites within a 12-km span, α2 diversity); and by estimating the total number of species inhabiting each ecoregion (γ diversity). The latter was calculated using an equivalent approach to capture-mark-recapture, by compiling a list of all shallow reef species recorded in the Ocean Biogeographic Information System (OBIS), FishBase, and SeaLifeBase for each ecoregion and then applying a correction for the proportion of unobserved species by cross-validation with the RLS data for the same families of organisms. Patterns in biodiversity at site, local, and regional scales are first described, and then, biogeographic and environmental covariates are introduced to characterize drivers that are likely behind these trends.

Results and Discussion

Distribution of inshore reef biodiversity

At ecoregion and local spatial scales, species richness of all taxa considered, vertebrates alone, and invertebrates alone was consistently maximal in the Indo-Malayan “Coral Triangle” region between Australia and Southeast Asia, declining rapidly from this area with latitude and more gradually with longitude (Fig. 1). Although also highest in this region, local-scale invertebrate richness differed slightly from ecoregion patterns in exhibiting moderately high values in some cooler regions (Fig. 1H). When mapped worldwide for both all taxa combined and vertebrates, site richness patterns were similar to patterns observed at local and ecoregional scales (Fig. 1). However, patterns of invertebrate richness at the smallest scale investigated (site) were reversed, with the highest values apparent in cooler regions (Fig. 1I).

Fig. 1

Patterns of global marine diversity.

Geographical distribution of ecoregion richness (species richness per ecoregion), local richness (rarefied estimates of total richness for sites <12 km from each other), and site richness (mean species richness per transect) for vertebrates, invertebrates, and all taxa. Legends indicate upper bounds for species richness bins. Ecoregion totals are predictively modeled using random forest procedures, with data trained using observed data from 82 ecoregions and 16 environmental covariates for which data were globally available.

Vertebrate richness showed an overall negative relationship with invertebrate richness at the site level (R2 = 0.24, P < 0.001; fig. S2) but a positive relationship across ecoregions (R2 = 0.45, P < 0.001; fig. S2). Thus, whereas all taxa (R2 = 0.43 for log model, P < 0.001) and vertebrates (R2 = 0.55 for log model, P < 0.001) showed strongly positive relationships between ecoregion and site richness (fig. S3), ecoregion richness was unrelated to site richness for invertebrates (R2 = 0.01 for log model, P > 0.05; fig. S3). Local-scale richness (α2 diversity) exhibited patterns intermediate between site and ecoregion richness. Notably, at this scale, invertebrate richness tended to increase slightly with vertebrate richness rather than show a strong positive or negative relationship (R2 = 0.10, 0.05 > P > 0.01; fig. S2).

Although global patterns for all taxa reinforced well-known centers of biodiversity (, ), substantial variation was evident among vertebrate and invertebrate classes at the site level, underpinned by variability in abundance. Bony fishes dominated both vertebrate and all taxa data sets (96% of vertebrate species and 64% of all species observed on transects); consequently, patterns for all taxa, vertebrates, and bony fishes were almost identical (Figs. 1 and 2). However, in contrast to the classical global pattern shown by bony fishes as well as patterns for ecoregion richness (), site-scale hotspots in different classes rarely overlapped (Fig. 2). Hence, patterns based on one taxon cannot be assumed for others when richness is considered at spatial scales smaller than ecoregion scales.

Fig. 2

Global patterns of site richness per 50-m transect for 10 investigated marine animal classes.

Data are smoothed using random forest predictions. Crustacea includes only species in orders Decapoda and Stomatopoda because other orders did not surpass the 2.5-cm minimum size requirement. Legends indicate the upper bounds for species richness bins.

At the site scale, most marine ecoregions worldwide had exceptionally high richness of at least one macroscopic animal class (Fig. 2). For example, crinoid species were concentrated in the Coral Triangle area but were absent or near absent from surveys in other tropical realms, including the eastern Pacific and the Caribbean. In contrast, local-scale richness hotspots for chondrichthyan fishes (sharks and rays) were positioned in the subtropical Pacific; reptiles were positioned off north-eastern and north-western Australia; large crustaceans and gastropods were positioned in subpolar realms; sea stars were positioned off southern South America and north-western North America; urchins were positioned in the warm-temperate south Pacific, eastern Pacific, Mediterranean, and Japan; holothurians were positioned in the eastern Pacific and Japan; and cephalopods were positioned off warm-temperate Australia, South Africa, and South America.

When species richness was related solely to latitude, ecoregion richness was maximal at low latitudes but with lower values on the equator than at latitudes of +15° and −15° (Fig. 3, A to C), as recently documented for a range of marine groups (). This decline reflected higher richness in the Caribbean than in equatorial Brazil, where richness may have been limited by the narrow continental shelf and turbid outflow of the Amazon, which in turn limited the habitable area for coral reefs (Fig. 1). In contrast, highest richness along the east African coast occurred near and just south of the equator. Similar relationships with latitude were also evident at site and local scales, albeit with little, if any, decline near the equator for local richness (Fig. 3, D to F), and with this trend inverted for visually conspicuous mobile invertebrates relative to vertebrates at the site scale (Fig. 3, G to I).

Fig. 3

Latitudinal trends in abundance and richness.

Bivariate plots relating ecoregion richness (A to C), local richness (D to F), site richness (G to I), and site abundance (J to L) to absolute latitude for all taxa, vertebrates, and invertebrates. Latitude was calculated as the mean latitude of sites investigated within ecoregion. Abundance y axes are shown as log scale. Best-fit polynomial curves to the third order were assessed using Bayesian information criterion (BIC) and R2 values; order of polynomial is shown in parentheses.

Latitudinal variation in the abundance of animals observed on transects mirrored patterns in site richness for both vertebrates and invertebrates (Fig. 3, K and L). These patterns canceled each other when total abundance of all taxa was considered, such that no consistent pattern was evident in overall abundance with latitude (R2 = 0.03; Fig. 3J). Total vertebrate abundance declined steeply toward polar latitudes (R2 = 0.66; Fig. 3K), whereas invertebrate abundance trended in the opposite direction (R2 = 0.30; Fig. 3L).

Low site richness of fishes in subpolar regions, and of invertebrates in the tropics, partly reflected the fact that only one or two individuals were counted in some transects, and thus, site richness could not exceed that number. These low abundance and richness counts apparently reflect real phenomena at the site scale, given that transect counts were standardized for constant area across all sites, but there is the possibility of varying behavioral responses of species to predation risk and thus diel activity (). Visual surveys may substantially underestimate the actual abundances of animals on transects, particularly invertebrates living concealed in complex habitats, such as coral reefs. However, we suggest that these biases are unlikely to account for the huge variation in abundance between tropical and temperate zones. For example, surveys undertaken along 50-m × 1-m transects day and night using the RLS method across 23 sites in Galapagos revealed that the numbers of large mobile macroinvertebrates nearly double at night (mean ± SE, 43 ± 2 m−2) relative to day (mean ± SE, 24 ± 1 m−2) as a consequence of animals moving out from concealment. Nevertheless, this difference was orders of magnitude lower than latitudinal variation in invertebrate counts, which were markedly high in cold waters.

Drivers of biodiversity patterns across scales

We evaluated the strength of evidence for hypotheses outlining drivers of latitudinal diversity gradients using two approaches. We first created generalized linear models (GLMs) to identify the most important correlates of species richness at different scales. We then constructed structural equation models (SEMs), which united the subset of informative variables identified in GLMs in a single hypothesized causal network, to disentangle the direct and indirect influences of the primary potential drivers across spatial scales and major taxa.

We considered the following environmental and geographic predictors that recur in hypotheses explaining the latitudinal biodiversity gradient and additionally corrected for local abundance at the site scale (Table 1): (i) temperature, through potential influences on metabolic rates, evolution rates (that is, mutation, evolution, and extinction), or phylogenetically conserved thermal preferences (, ); (ii) primary productivity (, ); (iii) area of continental shelf and coral reef (, ); (iv) natural disturbance regime (); (v) climatic stability since the last glacial maximum (); (vi) geographic connectivity and fragmentation (, ); and (vii) human disturbance (). GLMs indicated that annual mean sea surface temperature, coral reef area within ecoregion, and nitrate and chlorophyll levels were the most influential factors associated with patterns of site richness (Table 2). At the ecoregion level, shelf area (invertebrates) and connectivity to other ecoregions (sea stars) were also important in individual class models. Local invertebrate richness was influenced by coastline length and area of continental shelf and richness of all taxa by annual variability in sea temperature.

Table 2

Predictors of richness and abundance identified by GLMs.

Mean values per ecoregion provide replicate sample data in models. Marginal R2 values (% contribution) are shown for predictors that explain >5% of variation in best models identified using Akaike information criteria (AIC). Predictors trending in the opposite direction to associated latitudinal hypothesis are shown in italics; predictors with direction consistent with hypothesis are in bold. Additional details are provided in table S1. na, not assessed; −, poor model fit.

Taxon	Mean abundance	Site richness	Site richness (N-corrected)	Local richness	Ecoregion richness
All taxa	+NO3 (6)	+SST (56)+Coral (7)	+SST (60)+Coral (10)	+Coral (48) −NO3 (7)−SST range (6)	+Coral (42)−NO3 (15)
Vertebrates	+SST (68)	+SST (85)	−NO3 (21)+Coral (10)	+SST (73)	−NO3 (45)+Coral (15)
Invertebrates	−Coral (33)+Chl (6)	−SST (33) +Chl (6)−Shelf (6)	−Stability (9)	+Coast (8)+Shelf (7)	+Coral (20) +Shelf (8)−NO3 (6)
Actinopterygii	+SST (68)	+SST (85)	+Coral (25)−NO3 (7)	na	+Coral (44)−NO3 (16)
Asteroidea	−SST (37)+Chl (14)	−SST (41)+Chl (9)	−	na	+Connect (7)
Echinoidea	+NO3 (23)−Coral (8)	+NO3 (9)−Island (9)	−Shelf (14)	na	+Coral (9)
Gastropoda	−Coral (34)	−SST (28)	−	na	+Coral (19) +Chl (7)
Crustacea	−Stability (31)+NO3 (8)	−Stability (42)	−Stability (9)+Coral (9)	na	+Coral (16)+SST range (6)

Thus, latitudinal gradient hypotheses other than temperature, productivity, and habitat area received inconsistent or no empirical support from the GLMs. About half of the modeled metrics with influences on richness at the site scale responded in the opposite direction to predictions from associated hypotheses (Table 2). Coral reef area and temperature were highly correlated (R2 = 0.67; table S2), and at least one of these factors positively contributed to nearly half of the models. Both predictors contributed significantly to all taxa site richness, indicating that area of coral reef added species to site richness above the contribution of temperature alone. Sea surface temperature was a significant predictor primarily in site-level models, and coral reef area was a significant predictor in ecoregion models.

Richness of the investigated invertebrate classes generally showed poorly defined relationships with predictors, whereas models explained a very high proportion of variance in bony fish richness at site (R2 = 0.87) and ecoregion (R2 = 0.67) scales (table S1). Despite the most important predictors differing between classes, the model fit was good for all taxa combined (R2 = 0.66 for sites and 0.62 for ecoregions; table S1).

The GLM models that excluded the influence of abundance (by considering the residuals of the relationship between abundance and richness) generated an unexpected result. We found a much better model fit for site richness of all taxa (R2 = 0.73; table S1) than the model fits for any of the classes individually (best R2 value of 0.36 for any single class—bony fishes; Table 2 and table S1). Thus, within a standardized sample of individuals at a site, total richness was much more predictable than fish richness, which was more predictable than the numbers of species in other classes investigated. Given that the animal classes have different site-scale distributions globally (Fig. 2), we expected the opposite, with a reasonable fit for some classes to spatial and environmental covariates but a poor overall fit when aggregating across classes into total richness because of the likelihood of patterns canceling each other; therefore, we expected species numbers to generate a noisy overall pattern when summed. Ecological interactions that limited species richness at the site level were presumably occurring; if one mobile fish or invertebrate class has elevated richness within the community, then low richness was needed in others to maintain a predictable sum.

Nevertheless, not all classes contributed with equal weight to the predictable richness totals; fish richness probably affected invertebrate richness much more than invertebrate richness affected fishes. The strength of the relationship between bony fish richness and environmental variables at the local scale was so remarkable (R2 = 0.87), and thus, the unexplained residual variation was so low that the influence of other factors, including contributions of other animal classes, could not be great. This is consistent with a strong predatory influence of fishes on invertebrate species [see below, the studies of Bertness et al. (), Leleu et al. (), Ogden and Lobel (), Hay (), and Hughes ()] and with reduced activity and consumption among invertebrates when predatory fishes are present. Among herbivores, the mobility of fishes may also confer a competitive advantage over invertebrates in some situations by allowing rapid access to patchy algal resources. However, contrary competitive outcomes are known, with Caribbean urchins reported to deplete algae below levels needed by herbivorous fishes (–). Experimental research on cross-phyla interactions between grazers and predators is urgently needed to ascertain the context dependency of competitive mechanisms limiting local invertebrate richness.

Although we found no significant effect of human population density on the biodiversity metrics studied, this result relates to a single index of human impact (population density), and to richness only, and therefore should not be taken as indicating that fish communities were unaffected by human activity. Our results contrast with published studies focused on other reef community metrics, notably total fish biomass, where a substantial human footprint has been detected (, ). Our lack of correlation with richness metrics may also have concealed major changes in the species’ relative abundances and composition associated with human activity that can occur regardless of constancy of species richness totals (). Alternatively, effects of exploitation and habitat degradation may extend across nearly all sites, obscuring an effect of human population density on richness.

The GLM results suggested that water temperature, nitrate, and coral reef area were consistent predictors of biodiversity across scales, but they did not indicate whether these variables were acting directly on richness or indirectly, such as through abundance. A structural equation modeling approach revealed drivers operating in different directions at site and ecoregion scales and only minor connections between these scales (Fig. 4). The direct link between site and ecoregion richness for all taxa was weak despite a highly significant correlation between these two metrics (R2 = 0.35, P < 0.001; fig. S3), which are both positively associated with the same extrinsic factor, coral reef area.

Fig. 4

Most likely causal network showing hypothesized links between environmental drivers and diversity across scales.

Width of arrows indicates the magnitude of standardized linear coefficients in SEMs. Black arrows are positive coefficients, and red arrows are negative coefficients. Dashed arrow between local and ecoregion richness for all taxa indicates a link at the margins of significance (P = 0.066), albeit with a moderately high linear coefficient value (2.1). Model fit is good in all cases, indicated by nonsignificant χ2 values.

For all taxa, temperature exerted its strongest influence on richness through coral reef area, which directly enhances site and ecoregion richness and indirectly influenced local richness through aggregation of site richness (Fig. 4A). Changes in sea temperature, occurring with climate change, are thus likely to act directly at site scales on species abundance and richness and indirectly at broader scales because of effects on coral reef availability, which in turn contributes habitat structure and complexity to the seascape (). A world with few coral reefs would have a much lower mean ecoregion richness and probably global richness because about 30% of all marine fish species are associated with coral (), although the number of those dependent on coral is unknown.

Temperature also showed a direct positive effect on site richness (Fig. 4A), thus adding to the growing empirical support for a metabolic-related contribution to spatial variability in biological diversity (, , , ). Nevertheless, no direct influence of temperature was evident at the ecoregion scale when effects of reef area were also considered, suggesting that faster rates of evolution and extinction recently described for fish lineages from tropical relative to extratropical regions may not be reflected in contemporary latitudinal patterns ().

Nitrate showed opposing relationships to site abundance and ecoregion richness (Fig. 4A). Among sites, high free nitrate may reflect high inputs or low primary production due to other factors limiting plant growth (for example, iron, light, and temperature). If high, inputs would be expected to fuel primary productivity, which in turn increases total invertebrate abundance (Fig. 4C) and hence total site richness. Alternatively, whereas localized upwellings may occur within nutrient-depauperate systems, potentially supporting both high nutrient– and low nutrient–associated species interspersed within the broader seascape, low-nutrient sites (and their associated species) are unlikely within high-nutrient ecoregions. Further work is needed to clarify the functional basis of broad-scale relationships involving nitrate, including interactions with plant productivity, other macro- and micronutrients, and coastal pollution (for example, runoff of nutrients from watersheds).

Patterns of vertebrate richness were overwhelmingly related to temperature, which influenced site richness both directly and indirectly, the latter through an extremely strong influence on abundance and also through coral reef area (Fig. 4B). Multiple lines of evidence suggest that temperature positively contributes to richness through physiological processes over ecological time scales. First, sharp changes in species richness were observed over short distances with steep spatial temperature gradients, such as the 45% decrease in site richness over the 100-km distance between warm tropical and cool upwelling waters in Galapagos (). Second, model fits linking temperature with site richness (R2 = 0.90; Fig. 4B) and local richness (R2 = 0.84; Fig. 4B) are exceptionally high, leaving little variation for unexamined historical factors to account for. Regardless of regional rates of speciation and extinction (), vertebrate species must rapidly redistribute in response to temperature at the local scale, as also evident in ecological studies at locations with rapidly changing climate (). This presumably occurs through metabolic rate–dependent mechanisms that influence physiology, growth, dispersal, demography, habitat selection, species interactions, and, ultimately, abundance and occupancy patterns. In a rapidly warming world, a close relationship between temperature and vertebrate richness suggests that climate change may quickly overshadow historical constraints on diversity, a phenomenon already evident in numerous poleward range extensions ().

Nitrate exhibited little relationship with patterns of vertebrate richness and showed opposing relationships with invertebrate richness at site and local scales. Coral reef area also had opposing relationships with invertebrate richness at small and large scales, including an extremely strong influence at the ecoregion level. A negative relationship between coral reef area and site abundance may be due to undercounting by divers of invertebrates that are responding to fish predation pressure through concealment in a structurally complex coral habitat (). For example, sea urchins and other invertebrate grazers have been observed to adopt a more cryptic behavior and decrease in abundance in the presence of predators in Australia, Mediterranean, North East Pacific, and New Zealand [reviewed in the study of Leleu et al. ()]. In addition to higher rates of fish predation on invertebrates at tropical temperatures, fishes have also been reported to dominate herbivory on coral reefs worldwide because of their mobility and multiple foraging strategies, except where numbers are reduced by fishing (–, ). Where herbivorous fishes occur in warm-temperate seas, their impact on algae is generally patchy or seasonal ().

An overarching conceptual model of the latitudinal diversity gradient for large mobile reef-associated taxa that accounts for the above findings needs to consider the following: (i) negligible variation in total animal abundance with latitude but large variation in response to site-scale influences (such as predation, competition, topographic complexity, and wave exposure); (ii) fish abundance increasing with temperature; (iii) abundance of visually conspicuous mobile invertebrates increasing with nutrients but decreasing with site fish abundance (presumably due to negative ecological interactions); (iv) total site richness increasing with site abundance, temperature, and coral reef area; (v) local richness increasing with site richness; and (vi) ecoregion richness increasing with coral reef area, continental shelf area, historical evolutionary processes, and, to a minor extent, local richness (Fig. 5). The strong latitudinal gradient for all taxa at the site level, despite little consistent latitudinal pattern in total abundance, was due to a steeper slope of the abundance/richness relationship for fishes compared to other classes. For a hypothetical site with 1000 individuals, the predicted total species richness was 52, 21, 11, and 3 if only bony fishes, sea stars, gastropods, and echinoids, respectively, were present, based on regressions shown in fig. S4. Thus, the pronounced decrease in all taxa richness with latitude was directly related to decreasing proportional abundance of the most speciose group, bony fishes.

Fig. 5

Proposed model of global marine diversity.

At the site scale, temperature and nutrients influence abundance, which affects site richness, which in turn strongly influences local richness. Fishes control abundances of large mobile invertebrates through predation, generating a negative relationship between vertebrate and invertebrate richness at the site scale. At the ecoregion scale, species richness is influenced by local richness, the extent of coral reef, and biogeographic factors. [Top and bottom photos by G.J.E. and middle photo by R.D.S.-S. (University of Tasmania)].

Our SEM models depict upscaling processes, with ecoregion richness dependent on local richness, which is dependent on site richness. Richness is thereby viewed here as driven from small spatial scales reflecting the aggregation of individual sampled sites to provide regional totals. A more common paradigm in macroecology is that local communities are filtered subsets of broader species pools, with regional richness driving site richness (). We statistically compared these two alternative frameworks. Upscaling drivers generated a more likely configuration than downscaling drivers in SEMs (AIC = −357 versus −351, respectively), explained a higher proportion of the variance in most cases, and produced overall stronger relationships among observed variables based on standardized regression coefficients (fig. S5). The downscaling “all taxa” and “vertebrate” models also significantly departed from good model fit (df = 4; χ2 = 9.53, 21.94, and 7.18; P = 0.049, P < 0.001, and P = 0.126, respectively, for all taxa, vertebrate, and invertebrate SEMs). Nevertheless, differences in model diagnostics were relatively slight, with both downscaling and upscaling processes probably operating concurrently.

Conservation implications

The paths, and presumably the mechanisms, through which temperature influences biodiversity of mobile reef fauna greatly differed between ecoregion and local scales. Temperature appears to primarily affect species richness at ecoregion scales through the addition of coral reef area in warm latitudes and at local scales through metabolic processes. These differences highlight a need for investigations linking changing biodiversity to climate to consider the complicating effects of scale. Caution is clearly needed when inference on the rate of change in regional species pools is based on changes measured through time across a set of small sites, particularly when rates and scales of global change are estimated [for example, the study of Dornelas et al. ()]. Observed changes in richness at sites are likely to occur more rapidly than changes at the ecoregion scale, and rates of species gain relative to loss may also differ with geographic scaling ().

Our study indicates that ecological interactions affect the relative abundance and richness of at least visually conspicuous mobile invertebrates at local scales. Thus, monitoring of biodiversity needs to encompass the range of guilds and trophic levels in a community and not be limited to “indicator” groups that are strongly taxonomically constrained, such as fishes. The fact that interactions between classes were not visible at the ecoregion level may be a consequence of species richness at this spatial scale reflecting evolutionary drivers, such as temperature, more than local factors. Studies encompassing wider habitat and species diversity should further inform the patterns reported here, as should studies clarifying the influence of potential biases, including whether visual census methods result in systematically low estimates of mobile invertebrate richness in the tropics due to predator avoidance and nocturnal foraging.

Differences in the processes operating at large and small spatial scales indicate a need for comprehensive conservation planning that takes these scale differences into account. In addition to ecoregional representation, conservation practice would benefit from recognition of local-scale patchiness, whereby biodiversity elements can be concentrated in scattered locations rather than coinciding with the boundaries of ecoregion maps. Examples of integrated practice are BirdLife International’s “Important Bird Area” and “Endemic Bird Area” approach (), the Nature Conservancy’s “coarse” and “fine” filters (), and the International Union for Conservation of Nature’s “Key Biodiversity Areas” (), which all cover ecoregional representation while also including species with small distributions and thus patchiness at the local scale, which may otherwise be overlooked.

Although species loss is likely to be slower at ecoregion than site scales, rapid regional loss in richness can potentially occur with loss of coral reef area, which is indicated by SEM output to play a fundamental role in modulating biodiversity patterns at ecoregion scales. Ecoregions with large areas of coral reef have substantially higher total richness than ecoregions with a similar climate but little coral. From this, we confirm that the increasing frequency of catastrophic coral bleaching events (, ), which now extend globally, represents a major immediate threat to global patterns of shallow reef biodiversity.

Although a primary focus of conservation efforts on tropical coral reefs appears warranted given the importance of coral area as a global structuring agent, we conclude that initiatives are additionally needed to safeguard the large proportion of global biodiversity that lies in temperate and polar realms, particularly for invertebrates. Hotspots of local-scale biodiversity occur outside the tropics for many animal classes (Fig. 2) yet remain inadequately protected from anthropogenic pressures (). In particular, mobile macroinvertebrates were negatively associated with temperature at the site scale and are perhaps more susceptible to warming trends than previously thought. A potential effect of warming may be the suppression of visually conspicuous mobile invertebrate diversity due to increasing richness, activity, and abundance of predatory and herbivorous fishes. With “business as usual,” we predict mobile invertebrates to suffer substantial population declines and site-scale species losses, at least, with global warming. Given the geographic extent and magnitude of ongoing change, integrated transnational conservation strategies encompassing the global extent of biodiversity need urgent consideration.

MATERIALS AND METHODS

Reef vertebrate and invertebrate data

Standardized data were obtained by trained scientific and recreational divers participating in the citizen science RLS program, through underwater visual census of reef fishes, higher vertebrates, and mobile invertebrates along 50-m transect lines. Full details of fish census methods are provided in previous studies (, , ), and an online methods manual (www.reeflifesurvey.com) describes all data collection methods, including that for invertebrates. Data quality and training of divers are detailed in the study of Edgar and Stuart-Smith (). All observed vertebrate species were counted in 5-m-wide transect blocks (250 m2), and cryptic fishes and mobile invertebrates (total length, >2.5 cm) were counted in 1-m-wide transect blocks (50 m2). Data from the two methods were combined for a 50-m2 block area by randomly subsampling 20% of all individuals observed and removing any cryptic species present (defined by family in the online methods manual) in the 250-m2 blocks to avoid double counting, then adding these data to the corresponding 50-m2 data. Site richness was the total number of species per 50 m2.

A total of 13,858 transect blocks were censused at 2406 sites in 82 of 232 marine ecoregions worldwide (); 2516 bony fish, 66 sharks and rays, 8 reptile, 6 mammal, 105 holothurian, 104 echinoid, 180 asteroid, 99 crinoid, 805 gastropod, 28 cephalopod, and 210 malacostracan crustacean taxa were recorded.

Estimation of local richness

Local richness was assessed using rarefaction methods for groups of sites located within a 12-km maximum span, a classification that resulted in local groups with a mean separation of 4 ± 2 (SD) km between sites. We initially applied two rarefaction methods, calculating the Chao2 estimator for replicated incidence data () and extrapolation to asymptotic richness (1000 samples) using coverage-based rarefaction (). Calculation were made using the vegan () and iNEXT () packages in R. Incidence- and coverage-based estimators showed near 1:1 agreement for all taxa (R2 = 0.99). Consequently, we applied the incidence-based method because the coverage method produces less consistent results when richness is low.

However, estimates of local richness were affected by the number of transects surveyed within each location (R2 = 0.24, 0.11, and 0.37, for all taxa, vertebrates, and invertebrates, respectively). Assumptions underlying the rarefaction methods were contravened, as is the case generally with field studies, where regional richness is estimated through extrapolation. The pool of species within each location is not bounded, species are not randomly dispersed within each location, and different species have differing probabilities of observation when present. Consequently, as survey incidence increased, so did total estimated local richness due to continued addition of rare species, particularly vagrant species at edges of latitudinal and longitudinal ranges and depth and habitat boundaries. To minimize this source of bias, we first calculated the linear regression of local richness versus the number of surveys in locations. We then subtracted predicted richness based on the regression from observed richness for each site grouping and added the predicted richness for groupings with 200 sites surveyed. Thus, our local richness estimates are standardized to represent total richness for 200 sites within a maximum 12-km span.

Estimation of ecoregion richness

Ecoregion richness values were estimated using a capture-mark-recapture process by cross-validation of tallies of recorded species in online databases with RLS species counts. Species in online databases were analogous to initial captures, whereas RLS records for the same ecoregion comprised independent recaptures (). Species lists for each ecoregion were compiled using online presence records filtered by excluding species categorized as pelagic, deepwater (not recorded if shallower than 30 m), freshwater, and ecoregions without RLS surveys. To align online data with shallow reef habitat, we also excluded species in families not recorded during RLS surveys. Databases used comprised the OBIS (www.iobis.org/), FishBase (www.fishbase.org/), and SeaLifeBase (www.sealifebase.org/). The OBIS data used comprised 52,932 records of 14,336 species when duplicate records of species in an ecoregion are excluded. The FishBase (fishes) and SeaLifeBase (invertebrates) data comprised a total of 52,575 records of 11,807 species. The three sources were combined, and duplicates were removed to provide a data set of 18,452 species. RLS data comprised 18,602 ecoregion records of 4127 species. Raw ecoregion tallies derived from online databases are greatly affected by sampling effort. A correction was applied by calculating the species lists for each ecoregion using the OBIS, ReefBase, and SeaLifeBase data sets combined, then multiplying the total species count by an independent estimate of the proportion of species not observed. This proportion was calculated using the RLS data, as the inverse of the percent of species recorded by RLS divers in each ecoregion that were not included in the combined OBIS/FishBase/SeaLifeBase data set. Thus, estimates of species richness within each ecoregion were independent of RLS totals because the contribution of RLS data was in the form of proportions and thus unrelated to total counts.

Predictive mapping of global biodiversity

Random forest models (), which are available in the “extendedForest” packages for R (https://r-forge.r-project.org/projects/gradientforest), were developed to predict the distribution of species richness in inshore habitats globally. Procedures were similar to those detailed in the study of Edgar et al. (). Each random forest consisted of 2000 regression trees, where each tree was fit to a bootstrap sample of the biological data.

Relationships were identified between the mean richness of major taxonomic groupings observed per site in 82 marine ecoregions and associated values for the absolute value of latitude, longitude, hemisphere (N or S), human population density, eight remotely sensed oceanographic covariates available through Bio-ORACLE () [nitrate, phosphate, silicate, chlorophyll, photosynthetically active radiation, sea surface temperature annual range, mean sea surface temperature, and SD of sea surface temperature; see the study of Edgar et al. ()], and five map-derived variables (continental shelf area at a depth of 50 m, coral reef area, island frequency, coastline length, and isolation index; Table 1). Ecoregions with a value of zero for a diversity metric (for example, local reptile richness) were removed from analysis and treated as missing values when generating predictive models. To estimate predictor importance (percent change in accuracy) and model fit (R2), we used cross-validation where observations that were not selected in bootstrap samples for trees were compared to predictions. Only the six most important predictor variables were used in the predictive models, assessed by the percent change in accuracy when a predictor variable was removed from a model. The fit of the model for each diversity metric was as follows (R2 for site, local, and class levels, respectively): Actinopterygii (88%, na, 69%), Chondrichthyes (21%, na, 34%), Reptilia (15%, na, 34%), Malacostraca Crustacea (25%, na, 14%), Asteroidea (47%, na, 0%), Crinoidea (51%, na, 32%), Echinoidea (51%, na, 1%), Holothuroidea (36%, na, 11%), Cephalopoda (10%, na, 31%), Gastropoda (39%, na, 26%), vertebrates (80%, 72%, 63%), invertebrates (4%, 27%, 23%), and all taxa (75%, 66%, 56%). Local richness was not assessed (na) at class level because of low total counts in many ecoregions. Mammals were only sighted at sites in five high-latitude ecoregions; hence, their distribution was not modeled.

Assessing influence of predictors using GLMs

The influence of explanatory variables representing hypotheses driving latitudinal biodiversity gradients of abundance, site richness, local richness, and ecoregion richness was tested within the framework of GLMs (Gaussian family) using data for 2406 sites and 80 ecoregions, after poorly sampled outlying ecoregions in southern New Zealand and North and East Barents Sea were excluded (four sites) from the 82 ecoregions with RLS data. Response and independent variables were tested for normality using the Shapiro test and a priori–transformed using log (+1 where necessary) or square root transformations. Best models were identified with forward stepwise model building procedures based on BIC (). The presence of extreme outliers and nonlinear relationships was examined using biplots before each independent variable was included in the model. The proportion of variance explained by the best model, and the increase with the addition of each independent variable (marginal R2), was calculated as R2 = 1− (model deviance/null deviance). The influence of abundance on patterns of local richness was determined through analyses of the residuals of a linear regression model of abundance on richness.

Structural equation modeling

We used SEMs to statistically test the relationship between environmental variables and diversity indices across scales (site, local, and ecoregion) and between broad taxonomic groups (all taxa, vertebrates, and invertebrates). SEMs unite multiple variables in a single hypothesized causal network and thus provide a powerful tool to test multivariate hypotheses under a unified framework (). We also constructed competing models to test whether the influence of variables propagates upward from local to regional, or vice versa. Data used in analyses were averaged at the ecoregion level, and ecoregions with fewer than n = 3 observations were removed, yielding a total of N = 74 replicates used in analyses.

To statistically evaluate the SEMs, we used global variance-covariance estimation procedure implemented in the lavaan package () in R (). Variables were assessed for collinearity before model fitting but were considered to be sufficiently uncorrelated as to not bias estimates of SEs. Individual model fit was assessed using a χ2 goodness-of-fit test, where the model is considered to reproduce the data well if the resulting significance value is P > 0.05. All variables were log10-transformed before fitting to improve model fit. Standardized regression coefficients are reported. AIC were used to compare models of different structures (site → local → regional, regional → local → site), and models were considered to be significantly more likely if the difference between their AIC values was >2.