Scientific Reports

Sci Rep. 9: 9303

Latitudinal and bathymetrical species richness patterns in the NW Pacific and adjacent Arctic Ocean

Hanieh Saeedi

hanieh.saeedi@senckenberg.de

¹²³, http://orcid.org/0000-0003-2362-0328Mark J. Costello⁴, http://orcid.org/0000-0002-8747-2451Dan Warren⁵, Angelika Brandt¹²

1. Department of Marine Zoology, Senckenberg Research Institute and Natural History Museum, Senckenberganlage 25, 60325 Frankfurt am Main, Germany

2. 0000 0004 1936 9721grid.7839.5FB 15 Biological Sciences, Institute for Ecology, Evolution and Diversity, Goethe University Frankfurt, Max-von-Laue-Str. 13, 60438 Frankfurt am Main, Germany

3. OBIS data manager, Deep-Sea Node, Frankfurt am Main, Germany

4. 0000 0004 0372 3343grid.9654.eInstitute of Marine Science, University of Auckland, Auckland, 1142 New Zealand

5. Senckenberg Biodiversity and Climate Research Centre, Senckenberganlage 25, 60325 Frankfurt am Main, Germany

© The Author(s) 2019

DOI: 10.1038/s41598-019-45813-9

Published online: 26 June 2019

Abstract

Global scale analyses have recently revealed that the latitudinal gradient in marine species richness is bimodal, peaking at low-mid latitudes but with a dip at the equator; and that marine species richness decreases with depth in many taxa. However, these overall and independently studied patterns may conceal regional differences that help support or qualify the causes in these gradients. Here, we analysed both latitudinal and depth gradients of species richness in the NW Pacific and its adjacent Arctic Ocean. We analysed 324,916 distribution records of 17,414 species from 0 to 10,900 m depth, latitude 0 to 90°N, and longitude 100 to 180°N. Species richness per c. 50 000 km² hexagonal cells was calculated as alpha (local average), gamma (regional total) and ES50 (estimated species for 50 records) per latitudinal band and depth interval. We found that average ES50 and gamma species richness decreased per 5° latitudinal bands and 100 m depth intervals. However, average ES50 per hexagon showed that the highest species richness peaked around depth 2,000 m where the highest total number of species recorded. Most (83%) species occurred in shallow depths (0 to 500 m). The area around Bohol Island in the Philippines had the highest alpha species richness (more than 8,000 species per 50,000 km²). Both alpha and gamma diversity trends increased from the equator to latitude 10°N, then further decreased, but reached another peak at higher latitudes. The latitudes 60–70°N had the lowest gamma and alpha diversity where there is almost no ocean in our study area. Model selection on Generalized Additive Models (GAMs) showed that the combined effects of all environmental predictors produced the best model driving species richness in both shallow and deep sea. The results thus support recent hypotheses that biodiversity, while highest in the tropics and coastal depths, is decreasing at the equator and decreases with depth below ~2000 m. While we do find the declines of species richness with latitude and depth that reflect temperature gradients, local scale richness proved poorly correlated with many environmental variables. This demonstrates that while regional scale patterns in species richness may be related to temperature, that local scale richness depends on a greater variety of variables.

Introduction

The latitudinal and bathymetrical gradients of marine species richness have been widely studied at both regional¹–⁵ and global scales⁶–¹⁰. Recent studies showed that the global latitudinal richness gradient in most marine species follows a bimodal pattern correlated with sea surface temperature⁶,⁸,¹¹–¹⁴. That is, richness was highest in the tropics but it dips at the equator. In general, present marine species richness gradients decline from mid to high latitudes and from shallow to deep sea in many taxa¹,⁷,¹¹,¹⁵–¹⁷. However, diversity in some deep-sea taxa such as gastropods and nematodes increases from the continental shelf to the bathyal and abyssal zones due to increased environmental stability¹⁸–²⁰. The deep sea is almost two-thirds of the Earth, and over 84% of the ocean area is deeper than 2,000 m⁷,²¹,²². In contrast to shallow waters, in the deep sea chemical energy and carbon flux mostly control the species diversity, and temperature does not predict variation in rarified diversity in many taxa (e.g., Bivalvia and Gastropoda)²²,²³. Climate change can alter deep-sea latitudinal diversity gradients, even at tropical latitudes²⁴. The latitudinal gradients in species richness in the deep sea were generally present for the last 36 million years, but were weakened or absent during glacial periods²⁴. Thus, considering both latitude and deep-sea gradients of marine species richness together in the same geographic region seems overdue.

Long-term global and more recent regional processes are likely to drive the marine species richness. For example, major historic events (e.g. glaciation and plate tectonics) and evolutionary processes such as origination, dispersal (range expansion), and extinction, are considered important driving factors shaping the current latitudinal patterns of richness of marine species¹⁶,²⁵–²⁹. Global climatic constraints resulting from plate tectonics modulate ocean circulation, resulting in changes in surface water characteristics as well as altering connectivity between populations³⁰. In addition to continental drift and sea level change, recent latitudinal marine species richness analyses have considered light (as photosynthetically active radiation (PAR)), sea surface temperature, and habitat (e.g., continental shelf) in shaping the latitudinal gradients for shallow water marine species (e.g., bivalves and gastropods less than 200 m; recent and fossil marine zooplankton)⁸,¹¹,¹⁶,¹⁷,³¹–³⁴. Because light and temperature directly influence biomass and/or abundance, diversity may then increase as a result of secondary population dynamics and/or evolutionary processes²⁹,³⁵,³⁶. Temperature, as a proxy for thermal energy, also enhances the utilization rate of chemical energy by organisms. Temperature may also influence diversity by allowing a greater range of energetic lifestyles at warmer temperatures (the metabolic niche hypothesis)³⁵,³⁶. Tropical warmer climates have thus increased metabolic scope and biodiversity by fostering greater population size and extinction resistance²⁹. This allows more species to inhabit specialized niches as a result of greater available energy, and generates faster speciation and/or lower extinction rates²⁹,³⁵.

The tropical and subtropical areas of the west Pacific host the highest number of marine species worldwide³⁷–⁴¹. They also have high topographic complexity, including large semi-enclosed seas, many islands, and deep-sea trenches. The regions high species richness may thus be due to high rates of speciation due to warm temperatures and repeated separations and reconnections of populations due to changing sea levels and continental drift. However, despite the high species richness and uniqueness of the NW Pacific, its latitudinal and depth gradients of marine species richness and their potential causes have not been studied. Here we show how species richness changes with latitude and depth, and consider potential explanatory factors including temperature, oxygen (dissolved and saturated), primary productivity, chlorophyll, current velocity, salinity, nitrate, ocean area, and sampling effort. We also considered the adjacent Arctic Ocean of the NW Pacific to discover how these patterns change towards the highest latitudes. If these latitudinal and depth gradients are largely temperature correlated it would suggest that other variables, including topographic complexity, had negligible influence on the evolution of the fauna.

Methods

Our study area included the NW Pacific and its adjacent Arctic Ocean from latitude 0 to 90°N, and longitude 100 to 180°N including 14 sea basins (Fig. 1). All geographic distribution records were extracted from Ocean Biogeography Information System (OBIS) (www.iobis.org) and Global Biodiversity Information Facility (GBIF) (https://www.gbif.org) (for citations of the datasets used, please see SI, Table S1). The extracted data were merged and duplicates excluded. All species names were matched against the World Register of Marine Species⁴² and synonyms reconciled. Distribution records were manually checked for suitability, and dubious records were either corrected (e.g., reversing latitude and longitude fields, duplicate records) or removed (e.g., fossil records). The final dataset consisted of 324,916 distribution records of 17,414 marine species (1,792 families) from 0 to 10,900 m depth (SI, Fig. S1). Moreover, all the species were categorized to shallow-water and deep-sea benthic and pelagic groups (SI, Fig. S2). The statistical software R 3.4.4 and ArcMap 10.5.1 were used to analyse the data and plot the graphs.

Figure 1View larger versionFigure 1. Study area located in the NW Pacific and adjacent Arctic Ocean (latitude: 0 to 90°N in latitude and longitude: 100 to 180°N) showing 14 Sea Basins with different colours. ArcMap 10.5.1 was used to create this figure (https://support.esri.com/en/products/desktop/arcgis-desktop/arcmap/10-5).

We used three measures of species richness for each c. 50 000 km² hexagonal cell, latitudinal band, and depth interval: (1) the total number of species per hexagon (alpha diversity); (2) the total number of species per latitudinal band (gamma diversity); and (3) the estimated richness per sampled hexagon using rarefaction (ES50). The rarefaction method was used to reduce the effect of sampling effort on species richness patterns by counting the number of species in a constant number of random samples. We considered each sample as a unique combination of date and location where one or more species were recorded. We used ES50 in the ‘vegan’ R package to estimate the species richness in 50 samples per hexagon⁴³. Analyses considered samples at spatial scales of hexagons and 5° latitudinal bands. The former capture the heterogeneity in the underlying data best, but contain data collected by a variety of methods and sample sizes which are not necessarily comparable (e.g., plankton, whales, turtles, sharks, tracked animals, sediment cores, trawls, SI, Table S1). Hexagonal species richness analysis introduces considerable variability and bias between hexagons. The aggregation of samples into latitudinal bands aims to smooth out these effects because the number of species in a 5° band would tend to reach an asymptote.

We extracted environmental factors including average temperature (°C), dissolved oxygen (mol.m⁻³), primary productivity (g.m⁻³.d⁻¹), chlorophyll (mg.m⁻³), current velocity (m⁻¹), saturated oxygen (%), salinity (PSS), and nitrate (mol.m⁻³) for shallow water records (0–500 m); and average temperature, dissolved oxygen, chlorophyll, current velocity, salinity, and nitrate for deep-sea records (>500 m) from Bio-ORACLE (http://www.bio-oracle.org/)⁴⁴,⁴⁵ and Global Marine Environment Datasets (GMED) (only saturated oxygen) (http://gmed.auckland.ac.nz/)⁴⁶. All the extracted environmental layers were at a 5 arcmin (c. 9.2 km) spatial resolution

We used generalized additive models (GAMs) to examine the impact of environmental predictors on number of species and ES50 on a per-hexagon basis. The mid-point of each c. 50 000 km² hexagonal cell was calculated and collated with the spatial resolution of the environmental variables. Due to the high incidence of zeros in our species count data, models were built using the negative binomial error distribution. We fitted models via restricted maximum likelihood, using the automatic predictor selection implemented in the mgcv package⁴⁷ to control the complexity of smooth terms. For each analysis we fitted an intercept-only model, which represented the null hypothesis that response variables were not explained by environment, spatial sampling bias, or spatial autocorrelation. For models built using number of species as the response variable, we used the total number of records for each locality as an estimate of sampling effort. Models using ES50 as a response variable excluded number of records as an explanatory variable, as ES50 calculations are themselves intended to control for sampling effort. To model the effects of spatial autocorrelation on predictor and response variables, we used a two-dimensional spherical spline on latitude and longitude of sampling sites⁴⁷.

For models using number of species as a response variable, we fitted one model using only spatial sampling bias, one using only spatial autocorrelation, and one using both sampling bias and spatial autocorrelation. We also fitted a model for each environmental predictor separately, and one that represented the combined impacts of all environmental predictors. Models built using environmental predictors also included sampling effort and the effects of spatial autocorrelation. Candidate sets of models using ES50 as a response variable were the same as for number of species, except for the exclusion of sampling effort. The GAMs relating species richness and ES50 to environmental predictors on a per-hexagon basis, as well as tables presenting detailed model selection results, are in SI, S1 (Species Counts and Environment Based on Hexagonal Cells).

Generalizing per-hexagon data to 5° latitudinal bands substantially reduced sample size, which limited our ability to fit complex functional responses. Therefore, for this data, we fitted generalized linear models (GLMs) using a Poisson error distribution, and used number of records per band to control for differences in sampling effort. Per-hexagon and 5° models were evaluated using the small sample size corrected Akaike Information Criterion (AIC⁴⁸), a statistical method used to choose models with optimal fit to the data while controlling for over-parameterization⁴⁹,⁵⁰. The models with lower AIC scores are those that demonstrated a better compromise between model fit and model complexity. A difference in AIC value (deltaAIC) of less than two is considered to be inconclusive when comparing models. GLMs relating species richness and ES50 values to environmental predictors are given in SI, S2 (Species Counts and Environment Based on 5° Latitudinal Bands) and Fig. S4.

Results

Species richness

Most records in the study area were in shallow depths, and only 8% of the hexagons had no reported samples (Figs 2 and S2). Of these, almost all (7% of hexagons) were in the tropical and subtropical NW Pacific (between 0 to 30°N). In contrast, for the deep sea, 62% of the hexagons had no data. The highest sampling effort in both shallow (11,000 records) and deep sea (1,000 records below 500 m) in the study area was in the Philippines around Bohol Island (10.05°N, 124°E; ocean area: 3,000 km²) (Fig. 2).

Figure 2View larger versionFigure 2. Total number of records, species, and ES50 calculated per c. 50 000 km² hexagonal cells in depths above and below 500 m. ArcMap 10.5.1 was used to create this figure (https://support.esri.com/en/products/desktop/arcgis-desktop/arcmap/10-5).

The highest alpha species richness (shallow water: 8,800 species; deep sea: 800 species) was in the Philippines around Bohol Island (Fig. 2). The ES50 species richness index ranges from 0 to 50 and accounts for the effect of the number of samples. The highest ES50 (49) was also observed in Bohol Island in shallow water. However, the highest ES50 (48) in the deep sea was recorded around Luzon Island (16.14°N, 122.40°E; ocean area: ~1,250 km²) in the Philippines. In general, the Bering, Japan, Philippine, Sulu and Celebes Seas, including some areas of subtropical NW Pacific Ocean, had the highest alpha species richness in the shallow and deep NW Pacific. The Arctic Ocean also had high alpha species richness, with some values of ES50 > 40.

Latitudinal gradients in number of records, alpha species richness, gamma richness, and ES50, increased from the equator to latitudes 5 or 10°N, then decreased, and further reached another peak at the highest latitudes (75 to 90°N) (Figs 3 and 4). Latitudes 5 to 10°N hosted the highest species richness and ES50 which was mostly belonged to the Philippine Sea (average sea surface temperature = 28 °C and bottom temperature = 15 to 18 °C) (Figs 3 and 4). That further peaks in alpha and gamma richness were not seen in average ES50 which indicates that those peaks were due to the number of samples (Fig. 4). The latitude 60 to 70°N had the lowest species richness and ES50, as this is the smallest ocean area in this region (about 500,000 km²). Apart from this decrease related to ocean area, ES50 was very similar across latitudes (Fig. 4).

Figure 3View larger versionFigure 3. The number of sampling records, alpha species richness (number of species per hexagon), and ES50 ± SE calculated per hexagon against latitude for both shallow water (0–500 m) and deep sea (>500 m) species. Solid lines shows order 6 polynomial trend lines regression trends over the latitude.

Figure 4View larger versionFigure 4. Gamma species richness (number of species per 5° latitudinal band), average ES50 ± SE, average sea surface (SST) and bottom (BT) temperatures (10⁴ ± SE), and ocean area against 5° latitudinal bands.

Both shallow-water and deep-sea Anthozoa and Chordata had their highest gamma species richness in the subtropical areas from 5 to 10°N (Fig. 5). For Mollusca, highest shallow-water species richness was around 5°N and from 30 to 40°N in deep sea. Surprisingly, both shallow-water and deep-sea Polychaeta had their highest gamma species richness in the Arctic from 70 to 80°N, but not the NW Pacific Ocean. Most of the distribution records for higher taxa were observed between latitude 20 to 40°N (SI, Fig. S3). Mollusca, Porifera, and Bryozoa were mostly distributed around 10° latitude. Cnidarians, Pisces, and crustaceans had their higher distributions from latitude 20 to 30°N.

Figure 5View larger versionFigure 5. Total number of species (Gamma species richness) against latitude for five selected taxa.

Environmental variables

For shallow waters, the model using salinity as the sole environmental predictor for number of species had the lowest AIC value, followed closely by the model with all environmental predictors (deltaAIC of 0.32) (SI, S1). Model selection tables for the response variables are included in SI, S1. Several other models were within deltaAIC < 2 from the best model, including models using temperature, productivity, chlorophyll, current, and oxygen saturation. However, with a deltaAIC of only 1.8 between the top model and the one containing only spatial autocorrelation, none of these results is particularly strong. Considering ES50 for shallow communities yields stronger results; the model with all environmental predictors had the lowest AIC value, with a deltaAIC > 25 between it and the model containing only spatial autocorrelation. A model for ES50 in shallow waters that contained only the effects of productivity had a deltaAIC of 0.06 when compared to the one with all environmental predictors, indicating that much of the predictive power of the environment for determining ES50 in shallow waters likely comes from the effects of productivity. No other model for ES50 in shallow waters had a deltaAIC < 2 when compared to the top model.

In deep waters, both species richness and ES50 per hexagon were correlated with the environmental predictors. For both response variables, the model with the lowest AIC was the one containing all environmental predictors, comparing those models to models with only spatial autocorrelation results in a deltaAIC of 10.0 for number of species and 6.7 for ES50. Much of the explanatory power in these models seems to come from a combination of salinity, temperature, and dissolved oxygen; additional models with deltaAIC < 2 included salinity and temperature as predictors for number of species and oxygen as a predictor of ES50. The difference in deviance explained by all models containing environmental predictors compared to models containing only spatial autocorrelation is not large. However, this does not indicate that these predictors are not important; many of the predictors show substantial spatial autocorrelation, and as such the effects of the two-dimensional smoother fit to latitude and longitude may tend to assume some of the predictive power of the environmental predictors.

Much of the deviance in ES50 and species richness across 5° latitudinal bands was explained by sampling effort (number of records), but model selection nevertheless shows significant effects of the environmental predictors in some cases (SI, S2 and Fig. S4). Model selection tables for the response variables are included in SI, S2. Salinity was the best predictor of number of species in shallow water, but the intercept-only model was within 0.42 deltaAIC of the top model for ES50 in shallow water, indicating little explanatory power of the environment. In deep-water communities the story is similarly ambiguous when looking over 5° latitudinal bands; current was the top predictor of species counts in deep waters, while nitrate was the top predictor of ES50. We do caution that the severely reduced sample size and lack of explicit spatial autocorrelation terms in these models likely makes these results less reliable than those for individual hexagons.

Depth

From a total number of 324,916 species distribution records, 83% were from shallower than 500 m: 0–50 m = 48,809; 50–100 m = 10,869; 100–200 m = 10,670; 200–500 m = 23,538; and for > 500 m there were 22,104 records (7%). All the species were divided into four groups including shallow water and deep-sea pelagic and shallow water and deep-sea benthic species. A total number of 352,969 distribution records (the number of records here are higher than stated before because some species were grouped in both pelagic and benthic, or/and shallow water and deep sea categories) of four groups were mapped (SI, Fig. S1). About half of all the records (174,182 records) belonged to shallow-water pelagic species and only ~3% (12,176 records) were classified as deep-sea pelagic species. About 21% (73,282 records) and 26% of (93,329 records) the records belonged to shallow-water and deep-sea benthic species, respectively. The total number of records per hexagon per 100 m depth intervals was highest from 0 to 500 m among all depths, and then decreased sharply with depth (Fig. 6).

Figure 6View larger versionFigure 6. The relationship of sampling effort (total number of records), alpha species richness (average number of species), gamma species richness (total number of species), and average ES50 per hexagon per 100 depth intervals. The solid lines show the logarithmic trends.

In the deep sea, sampling effort generally decreased with depth in 100 m depth intervals with the exception of two peaks around 2,000 and 4,000 m. Average ES50, alpha, and gamma species richness per hexagon decreased with depth (100 m depth intervals) and had the highest values from 0 to 2,000 m (Fig. 6).

Discussion

There were eight times less data available for the deep sea than coastal depths in our study area. Nevertheless, the data still showed a decline in species richness with depth when adjusted for sampling effort (ES50). Similarly, near coast areas had greater richness as ES50 than offshore areas (Fig. 2).

The most species rich area in the world ocean is the Coral Triangle, which includes the Philippine Sea (latitudes from ~5 to 35°N)⁷,³⁰,⁵¹,⁵². This is the most species-rich region for bivalves⁸,³⁷,³⁸,⁴⁰,⁵³, bony fish⁵⁴,⁵⁵, sharks⁵⁶, crustaceans⁵⁷, ascidians⁵⁸, anemones and corals⁵²,⁵⁹, benthic marine algae⁶⁰, endemic fish⁶¹ and other marine organisms⁶². Our data support these findings. We found that the Philippine Sea around Bohol Island (10.05°N, 124°E) had the highest alpha species richness of both the shallow and deep NW Pacific. This area is where the Philippine Sea has its highest number of islands and is a meeting point for the species rich Sulu, South China, and Celebes Seas. Even when corrected for sampling effort, this area and Luzon Island (16.14°N, 122.40°E) had the highest expected species richness in shallow water and deep sea, respectively. However, specific taxa might show different diversity patterns. For example, Ophiuroidea species richness in continental shelf to upper-slope areas peaks at tropical Indo-west Pacific and Caribbean (0–30°) latitudes, following the water temperature²². In contrast, deep-sea ophiuroid species show maximum richness at higher latitudes (30–50°), in regions close to continental margins where carbon export flux is high²².

Many coral reef groups reach their greatest diversity in the Coral Triangle. In addition to higher coral reef diversity in this area³⁰,⁶³, there are more reasons why the Coral Triangle is so species rich. Being tropical it has not suffered from glaciation driven extinctions. The tropics have higher rates of potential speciation due to warm temperatures increasing mutation rates and decreasing generation times (see Costello and Chaudhary 2017 for a recent review). Warmer sea surface temperatures, high productivity, and habitat availability and heterogeneity are likely important factors responsible for high tropical shallow species richness⁶,⁸,³⁹,⁶⁴. Our richest spots for species were in the Philippine Sea with a mean SST from 24 to 28 °C and ~20 °C in winter⁶,⁸,⁶⁴. The NW Pacific is relatively rich in nutrients from land and rainfall, enabling high productivity, large population sizes, and intra-specific competition that can drive speciation. The many islands and deep-sea areas create a diversity of bathymetric and oceanographic conditions to which species can evolve to adapt, and past fluctuations in sea level will have isolated populations in coastal areas that have subsequently become reconnected. In addition to present and past environmental conditions, this region serves as a meeting point for the biota of Asia and Australia, and the Indian and Pacific Oceans, including species that may have evolved in Australia as part of Gondwanaland³⁰,⁶⁵. Although the area has never been glaciated, it was definitely affected by Quaternary sea level cycles which could create a diversity pump mechanism, as a result of the relationship between Pleistocene sea level changes and the complex geography of the area, or as a result of accumulation, reflecting the exceptional environmental features of the region³⁰.

While looking at our data set in isolation, one may consider that the lower richness at the equator (0 to 5°N) was an artefact of sampling bias. However, a recent review has shown that marine species richness does decrease at the equator, where the highest sea temperatures occur, across almost all taxa¹¹. Thus the dip in richness may be due to thermal stress above an annual average SST of 28 °C. For example, bivalve larvae (e.g., Mytilopsis leucophaeata) have higher embryo and larval mortality in temperatures above 28 °C⁸,⁶⁶,⁶⁷, and annual average sea surface temperature was above 28 °C (29 to 30 °C) at the equatorial western Pacific from 1870 to 2005⁶⁸.

Considering the best model fit, both hexagon-based and 5° latitudinal band models showed that ES50 was the better estimate of species richness in both shallow and deep sea. We also concluded that the best model explaining species richness in the shallow and deep sea was when all the environmental factors were considered rather than single variables. In other words, no single variable explained the pattern in species richness for either hexagons or latitudinal bands. Alpha diversity (species richness per hexagon), and gamma diversity (species richness per 5° latitudinal band), showed different regression levels and outcomes in the model outputs. However, the high spatial heterogeneity (hexagons) resulted in poor correlations between species richness and environmental conditions. While clearer correlation trend lines were found when data were aggregated into 5° latitudinal bands, the reduced sample size limited the statistical outcome of the results.

Sampling effort and bias are important factors to consider in interpreting species richness patterns. Indeed, studies on global latitudinal species richness gradients showed that alpha, and to a lesser extent gamma, species richness patterns were affected by sampling effort⁸,¹¹. Our comparison of alpha, gamma, and ES50 (to account for sampling bias) showed that all measures still peaked between 5 to 10°N in both shallow and deep sea (Fig. 3). However, ES50 also showed high values in the Arctic Ocean, not indicated by alpha and gamma species richness. Average ES50, alpha, and gamma species richness decreased with depth below the lower continental slope in the study area, but ES50 richness trend was not decreased as sharp as alpha and gamma species richness (Fig. 6). Some bathymetrical peaks of species richness in the Arctic and lower-slope have been well documented by other studies²⁰,⁶⁹,⁷⁰.

We found that the NW Pacific conforms to the recent global findings of species declining with latitude and depth, such that most species occur in tropical coastal depths⁷. Furthermore, the dip in richness observed by Chaudhary et al. (2016, 2017) is also present from 0–5° latitude. By calculating alpha, gamma, and ES50 against latitude and depth, we showed that using different diversity indices may influence perceived patterns of species richness over large spatial scales.

Supplementary Material

Open In Web Browser

Footnotes

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary information accompanies this paper at 10.1038/s41598-019-45813-9.

Acknowledgements

This paper was part of the “Biogeography of the NW Pacific deep-sea fauna and their possible future invasions into the Arctic Ocean project (Beneficial project)”. Beneficial project (grant number 03F0780A) was funded by Federal Ministry for Education and Research (BMBF: Bundesministerium für Bildung und Forschung) in Germany. We would like to thank Chhaya Chaudhary for her collaboration in primary dataset preparation and Marianna Simões for reformatting Figure S4.

Author Contributions

H.S., M.J.C. and A.B. designed this study, H.S. performed all the statistical analyses, D.W. contributed to statistical analyses, H.S. wrote the first draft, and all authors discussed the results and co-wrote the manuscript.

Data Availability

The dataset will be available upon the request.

Competing Interests

The authors declare no competing interests.

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group

PMC Copyright Notice

The articles available from the PMC site are protected by copyright, even though access is free. Copyright is held by the respective authors or publishers who provide these articles to PMC. Users of PMC are responsible for complying with the terms and conditions defined by the copyright holder.

Users should assume that standard copyright protection applies to articles in PMC, unless an article contains an explicit license statement that gives a user additional reuse or redistribution rights. PMC does not allow automated/bulk downloading of articles that have standard copyright protection.

See the copyright notice on the PMC site, https://www.ncbi.nlm.nih.gov/pmc/about/copyright/, for further details and specific exceptions.

References

Rex MA, Stuart CT, Coyne G, authors. Latitudinal gradients of species richness in the deep-sea benthos of the North Atlantic. Proceedings of the National Academy of Sciences of the United States of America. 2000;97:4082–4085. [CrossRef] [PubMed]

Gray J, author. Marine diversity: the paradigms in patterns of species richness examined. Scientia Marina. 2001;65:41–56. [CrossRef]

Yasuhara M, et al., authors. Cenozoic dynamics of shallow-marine biodiversity in the Western Pacific. Journal of Biogeography. 2017;44:567–578. [CrossRef]

Crame J. Alistair, McGowan Alistair J.; Bell Mark A., authors. Differentiation of high-latitude and polar marine faunas in a greenhouse world. Global Ecology and Biogeography. 2018;27(5):518–537. [CrossRef]

Roy K, Martien KK, authors. Latitudinal distribution of body size in north-eastern Pacific marine bivalves. Journal of Biogeography. 2001;28:485–493. [CrossRef]

Chaudhary C, Saeedi H, Costello MJ, authors. Bimodality of Latitudinal Gradients in Marine Species Richness. Trends in Ecology & Evolution. 2016;31:670–676. [CrossRef] [PubMed]

Costello MJ, Chaudhary C, authors. Marine Biodiversity, Biogeography, Deep-Sea Gradients, and Conservation. Current Biology. 2017;27:2051 [CrossRef] [PubMed]

Saeedi H, Dennis TE, Costello MJ, authors. Bimodal latitudinal species richness and high endemicity of razor clams (Mollusca). Journal of Biogeography. 2017;44:592–604. [CrossRef]

Tittensor Derek P.; Mora Camilo, Jetz Walter, Lotze Heike K.; Ricard Daniel, Berghe Edward Vanden, Worm Boris, authors. Global patterns and predictors of marine biodiversity across taxa. Nature. 2010;466(7310):1098–1101. [CrossRef] [PubMed]

Rabosky DL, et al., authors. An inverse latitudinal gradient in speciation rate for marine fishes. Nature. 2018;559:392–395. [CrossRef] [PubMed]

Chaudhary C, Saeedi H, Costello MJ, authors. Marine Species Richness Is Bimodal with Latitude: A Reply to Fernandez and Marques. Trends in Ecology & Evolution. 2017;32:234–237. [CrossRef] [PubMed]

Saeedi H, Basher Z, Costello MJ, authors. Modelling present and future global distributions of razor clams (Bivalvia: Solenidae). Helgoland Marine Research. 2017;70:23 [CrossRef]

Saeedi, H. & Costello, M. J. In Reference Module in Earth Systems and Environmental Sciences (Elsevier, 2019).

Saeedi, H. & Costello, M. J. A world dataset on the geographic distributions of Solenidae razor clams (Mollusca: Bivalvia). Biodiversity Data Journal7, 10.3897/BDJ.7.e31375 (2019).

Danovaro R, Gambi C, Lampadariou N, Tselepides A, authors. Deep-sea nematode biodiversity in the Mediterranean basin: testing for longitudinal, bathymetric and energetic gradients. Ecography. 2008;31:231–244. [CrossRef]

Brown JH, Svenning JC, authors. Why are there so many species in the tropics. Journal of Biogeography. 2014;41:8–22. [CrossRef] [PubMed]

Valentine JW, Jablonski D, Crame A, authors. A twofold role for global energy gradients in marine biodiversity trends. Journal of Biogeography. 2015;42:997–1005. [CrossRef]

Danovaro R, et al., authors. Deep-Sea Biodiversity in the Mediterranean Sea: The Known, the Unknown, and the Unknowable. PLOS ONE. 2010;5:e11832 [CrossRef] [PubMed]

Rex MA, author. Deep-Sea Species Diversity: Decreased Gastropod Diversity at Abyssal Depths. Science (New York, N.Y.). 1973;181:1051–1053. [CrossRef]

Rex, M. A. & Etter, R. J. Deep-Sea Biodiversity: Pattern and Scale. (Harvard University Press, 2010).

Costello MJ, Cheung A, De Hauwere N, authors. Surface Area and the Seabed Area, Volume, Depth, Slope, and Topographic Variation for the World’s Seas, Oceans, and Countries. Environmental Science & Technology. 2010;44:8821–8828. [CrossRef] [PubMed]

Woolley SNC, et al., authors. Deep-sea diversity patterns are shaped by energy availability. Nature. 2016;533:393-+ [CrossRef] [PubMed]

McClain CR, Allen AP, Tittensor DP, Rex MA, authors. Energetics of life on the deep seafloor. Proceedings of the National Academy of Sciences. 2012;109:15366–15371. [CrossRef]

Yasuhara M, Hunt G, Cronin TM, Okahashi H, authors. Temporal latitudinal-gradient dynamics and tropical instability of deep-sea species diversity. Proceedings of the National Academy of Sciences. 2009;106:21717–21720. [CrossRef]

Myers AA, author. Biogeographic Barriers and the Development of Marine Biodiversity. Estuarine, Coastal and Shelf Science. 1997;44:241–248. [CrossRef]

Davies AJ, Guinotte JM, authors. Global Habitat Suitability for Framework-Forming Cold-Water Corals. PLOS ONE. 2011;6:e18483 [CrossRef] [PubMed]

Powell MG, Beresford VP, Colaianne BA, authors. The latitudinal position of peak marine diversity in living and fossil biotas. Journal of Biogeography. 2012;39:1687–1694. [CrossRef]

Berke Sarah K.; Jablonski David, Krug Andrew Z.; Valentine James W., authors. Origination and Immigration Drive Latitudinal Gradients in Marine Functional Diversity. PLoS ONE. 2014;9(7):e101494 [CrossRef] [PubMed]

Erwin DH, author. Climate as a Driver of Evolutionary Change. Current Biology. 2009;19:R575–R583. [CrossRef] [PubMed]

Renema W, et al., authors. Hopping Hotspots: Global Shifts in Marine Biodiversity. Science (New York, N.Y.). 2008;321:654–657. [CrossRef]

Hillebrand H, author. On the Generality of the Latitudinal Diversity Gradient. The American Naturalist. 2004;163:192–211. [CrossRef] [PubMed]

Fernández M, Astorga A, Navarrete SA, Valdovinos C, Marquet PA, authors. Deconstructing latitudinal species richness patterns in the ocean: does larval development hold the clue? Ecology Letters. 2009;12:601–611. [CrossRef] [PubMed]

Tomašových A, Jablonski D, Berke SK, Krug AZ, Valentine JW, authors. Nonlinear thermal gradients shape broad-scale patterns in geographic range size and can reverse Rapoport’s rule. Global Ecology and Biogeography. 2015;24:157–167. [CrossRef]

Yasuhara M, Hunt G, Dowsett HJ, Robinson MM, Stoll DK, authors. Latitudinal species diversity gradient of marine zooplankton for the last three million years. Ecology Letters. 2012;15:1174–1179. [CrossRef] [PubMed]

Clarke A, Gaston KJ, authors. Climate, energy and diversity. Proceedings. Biological sciences. 2006;273:2257–2266. [CrossRef] [PubMed]

Mittelbach GG, et al., authors. Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography. Ecology Letters. 2007;10:315–331. [CrossRef] [PubMed]

Alistair Crame J., author. Evolution of taxonomic diversity gradients in the marine realm: evidence from the composition of Recent bivalve faunas. Paleobiology. 2000;26(2):188–214. [CrossRef]

Krug AZ, Jablonski D, Valentine JW, authors. Species–genus ratios reflect a global history of diversification and range expansion in marine bivalves. Proceedings of the Royal Society B: Biological Sciences. 2008;275:1117–1123. [CrossRef] [PubMed]

Sanciangco JC, Carpenter KE, Etnoyer PJ, Moretzsohn F, authors. Habitat Availability and Heterogeneity and the Indo-Pacific Warm Pool as Predictors of Marine Species Richness in the Tropical Indo-Pacific. PLOS ONE. 2013;8:e56245 [CrossRef] [PubMed]

Jablonski D, Roy K, Valentine JW, authors. Out of the tropics: evolutionary dynamics of the latitudinal diversity gradient. Science (New York, N.Y.). 2006;314:102–106. [CrossRef]

Asaad I, Lundquist CJ, Erdmann MV, Costello MJ, authors. Delineating priority areas for marine biodiversity conservation in the Coral Triangle. Biological Conservation. 2018;222:198–211. [CrossRef]

WoRMS. World Register of Marine Species, http://www.marinespecies.org/aphia.php?p=stats (2018).

Oksanen, J. F. G. B. et al. vegan: Community Ecology Package, https://CRAN.R-project.org/package=vegan (2017).

Assis J, et al., authors. Bio-ORACLE v2.0: Extending marine data layers for bioclimatic modelling. Global Ecology and Biogeography. 2018;27:277–284. [CrossRef]

Tyberghein L, et al., authors. Bio-ORACLE: a global environmental dataset for marine species distribution modelling. Global Ecology and Biogeography. 2012;21:272–281. [CrossRef]

Basher, Z., Bowden, D. A. & Costello, M. J. Vol. Version 1.0 (Rev.01.2014) http://gmed.auckland.ac.nz (World Wide Web electronic publication, 2014).

Wood, S. N. mgcv: Mixed GAM Computation Vehicle with GCV / AIC / REML Smoothness Estimation, http://cran.r-project.org/web/packages/mgcv (2013).

Burnham KP, Anderson DR, Huyvaert KP, authors. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology. 2011;65:23–35. [CrossRef]

Akaike N, author. The origin of the basal cell potential in frog corneal epithelium. The Journal of physiology. 1971;219:57–75. [CrossRef] [PubMed]

Burnham, K. P. & Anderson, D. R. Model selection and. (2002).

Costello MJ, author. Biodiversity: The Known, Unknown, and Rates of Extinction. Current Biology. 2015;25:R368–R371. [CrossRef] [PubMed]

Bellwood, D. R., Renema, W. & Rosen, B. R. In Biotic Evolution and Environmental Change in Southeast Asia Systematics Association Special Volume Series (eds Brian Rosen et al.) 216–245 (Cambridge University Press, 2012).

Flessa KW, Jablonski D, authors. Biogeography of Recent marine bivalve molluscs and its implications for paleobiogeography and the geography of extinction: a progress report. Historical Biology. 1995;10:25–47. [CrossRef]

Carpenter KE, Springer VG, authors. The center of the center of marine shore fish biodiversity: the Philippine Islands. Environmental Biology of Fishes. 2005;72:467–480. [CrossRef]

Briggs JC, Bowen BW, authors. A realignment of marine biogeographic provinces with particular reference to fish distributions. Journal of Biogeography. 2012;39:12–30. [CrossRef]

Lucifora Luis O.; García Verónica B.; Worm Boris, authors. Global Diversity Hotspots and Conservation Priorities for Sharks. PLoS ONE. 2011;6(5):e19356 [CrossRef] [PubMed]

Neiber Marco T.; Hartke Tamara R.; Stemme Torben, Bergmann Alexandra, Rust Jes, Iliffe Thomas M.; Koenemann Stefan, authors. Global Biodiversity and Phylogenetic Evaluation of Remipedia (Crustacea). PLoS ONE. 2011;6(5):e19627 [CrossRef] [PubMed]

Shenkar Noa, Swalla Billie J., authors. Global Diversity of Ascidiacea. PLoS ONE. 2011;6(6):e20657 [CrossRef] [PubMed]

Fautin DG, Malarky L, Soberon J, authors. Latitudinal Diversity of Sea Anemones (Cnidaria: Actiniaria). Biological Bulletin. 2013;224:89–98. [CrossRef] [PubMed]

Kerswell Ailsa P., author. GLOBAL BIODIVERSITY PATTERNS OF BENTHIC MARINE ALGAE. Ecology. 2006;87(10):2479–2488. [CrossRef] [PubMed]

Asaad I, Lundquist CJ, Erdmann MV, Costello MJ, authors. Ecological criteria to identify areas for biodiversity conservation. Biological Conservation. 2017;213:309–316. [CrossRef]

Selig Elizabeth R.; Turner Will R.; Troëng Sebastian, Wallace Bryan P.; Halpern Benjamin S.; Kaschner Kristin, Lascelles Ben G.; Carpenter Kent E.; Mittermeier Russell A., authors. Global Priorities for Marine Biodiversity Conservation. PLoS ONE. 2014;9(1):e82898 [CrossRef] [PubMed]

Bellwood DR, Hughes TP, authors. Regional-scale assembly rules and biodiversity of coral reefs. Science (New York, N.Y.). 2001;292:1532–1535. [CrossRef]

Förderer M, Rödder D, Langer MR, authors. Patterns of species richness and the center of diversity in modern Indo-Pacific larger foraminifera. Scientific Reports. 2018;8:8189 [CrossRef] [PubMed]

Leprieur, F. et al. Plate tectonics drive tropical reef biodiversity dynamics. Nature Communications7, 11461, 10.1038/ncomms11461https://www.nature.com/articles/ncomms11461#supplementary-information (2016).

Talmage Stephanie C.; Gobler Christopher J., authors. Effects of Elevated Temperature and Carbon Dioxide on the Growth and Survival of Larvae and Juveniles of Three Species of Northwest Atlantic Bivalves. PLoS ONE. 2011;6(10):e26941 [CrossRef] [PubMed]

Verween A, Vincx M, Degraer S, authors. The effect of temperature and salinity on the survival of Mytilopsis leucophaeata larvae (Mollusca, Bivalvia): The search for environmental limits. Journal of Experimental Marine Biology and Ecology. 2007;348:111–120. [CrossRef]

Hansen J, et al., authors. Global temperature change. Proceedings of the National Academy of Sciences of the United States of America. 2006;103:14288–14293. [CrossRef] [PubMed]

Rex MA, Crame JA, Stuart CT, Clarke A, authors. Large-scale biogeographic patterns in marine mollusks: A confluence of history and productivity? Ecology. 2005;86:2288–2297. [CrossRef]

Yasuhara M, et al., authors. Patterns and controlling factors of species diversity in the Arctic Ocean. Journal of Biogeography. 2012;39:2081–2088. [CrossRef]

[Back]

Figure 1.

Study area located in the NW Pacific and adjacent Arctic Ocean (latitude: 0 to 90°N in latitude and longitude: 100 to 180°N) showing 14 Sea Basins with different colours. ArcMap 10.5.1 was used to create this figure (https://support.esri.com/en/products/desktop/arcgis-desktop/arcmap/10-5).

[Back]

[Back]

Figure 2.

Total number of records, species, and ES50 calculated per c. 50 000 km² hexagonal cells in depths above and below 500 m. ArcMap 10.5.1 was used to create this figure (https://support.esri.com/en/products/desktop/arcgis-desktop/arcmap/10-5).

[Back]

[Back]

Figure 3.

The number of sampling records, alpha species richness (number of species per hexagon), and ES50 ± SE calculated per hexagon against latitude for both shallow water (0–500 m) and deep sea (>500 m) species. Solid lines shows order 6 polynomial trend lines regression trends over the latitude.

[Back]

[Back]

Figure 4.

Gamma species richness (number of species per 5° latitudinal band), average ES50 ± SE, average sea surface (SST) and bottom (BT) temperatures (10⁴ ± SE), and ocean area against 5° latitudinal bands.

[Back]

[Back]

Figure 5.

Total number of species (Gamma species richness) against latitude for five selected taxa.

[Back]

[Back]

Figure 6.

The relationship of sampling effort (total number of records), alpha species richness (average number of species), gamma species richness (total number of species), and average ES50 per hexagon per 100 depth intervals. The solid lines show the logarithmic trends.

[Back]