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

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 km2 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 km2). 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.


Species Counts and Environment
Setup Shallow water

Data load-in and visualization
First we're going to load in our data and then trim the data frame down to just the columns we need.

Shallow water
Deep water shallow.numsp.intercept <-gam(num.species ~ 1, data = shallow, family = "nb", meth od = "REML", select = TRUE) gam.check ( GAMs for number of species, shallow water We're going to develop a number of GAMs here. The "intercept" GAM represents the fit of a model that assumes no relationship between the species counts and the environment and no spatial autocorrelation. The "latlon" model represents the fit of a model that fits spatial autocorrelation, but no environmental effects. The "env" model represents the combined effects of all environmental predictors, and the remainder of the models ("temp", "oxygen", etc.) estimate the effects of a single predictor at a time.
shallow.numsp.numrec <-gam(num.species ~ s(num.records), data = shallow, family = "nb", method = "REML", select = TRUE) gam.check ( .037 2.38e-05 *** ## --shallow.numsp.prod <-gam(num.species ~ s(lat, lon, bs = "sos") + s(num.records) + s(primary.productivity), data = shallow, family = "nb", method = "REML", select = T RUE) gam.check ( Here we use AIC to compare models for goodness of fit while penalizing for overparameterization. Models are ranked in order of AIC score, with lowest scoring model (model with the best compromise between fit and complexity) first in the list. The delta AIC represents how much the AIC score for each model differs from the top model, an can be used as an estimate of the relative support for each model. A delta AIC of 2 is considered to be potentially a significantly better fit, with higher delta AICs between models indicating increasingly larger differences in model fit while correcting for the number of parameters. Here we see that the model containing salinity is the best fit, but many other models are almost as good. It is particularly noteworthy that the model containing only the effects of spatial autocorrelation (the "latlon" model) has a delta AIC of < 2 compared to the best-performing model. However, we know that latitude is highly correlated with several of the environmental predictors. That being the case, it is hard to determine which, if any, environmental predictors are affecting the number of species in shallow water. T GAMs for ES50, shallow water These models are identical to the above, except that the response variable is the species richness estimated from the rarefaction curves. These are meant to be estimates of species richness if all sites were sampled equally. Since these measurements are intended to take sampling effort into account, we exclude the "no.records" term for these models.

Model
shallow.es50.intercept <-gam(ES50 ~ 1, data = shallow, family = "nb", method = "RE ML", select = TRUE) gam.check ( shallow.es50.latlon <-gam(ES50 ~ s(lat, lon, bs = "sos") , data = shallow, family = "nb", method = "REML", select = TRUE) gam. Here we find strong support for the effects of the environment on ES50 in shallow water. Given that primary productivity alone does almost as good a job as all environmental predictors at once, we believe that the predictive power here is mostly coming from primary productivity. The relatively poor performance of the spatial autocorrelation-only model (delta AIC of 29.6) suggests that there is significant predictive power in the environmental variables over and above what is expected due to spatial autocorrelation alone.

Deep water
Repeating the above analyses for deep-dwelling species.

Data load-in and visualization
GAMs for number of species, deep water As above, we're going to develop a number of GAMs including one for "intercept", one for "latlon", one for all environmental predictors ("env"), and then one for each individual predictor.

GLMs for number of species, shallow water
We're going to develop a number of GLMs here. The "intercept" glm represents the fit of a model that assumes no relationship between the species counts and the environment and no spatial autocorrelation.
The "numrec" model represents the fit of a model that only contains sampling bias. The remainder of the models ("temp", "oxygen", etc.) estimate the effects of a single predictor at a time, controlling for number of records. GLMs for ES50, shallow water These models are identical to the above, except that the response variable is the species richness estimated from the rarefaction curves. These are meant to be estimates of species richness if all sites were sampled equally. Since these measurements are intended to take sampling effort into account, we exclude the "num.records" term for these models.

GLMs for number of species, deep water
We're going to develop a number of GLMs here. The "intercept" glm represents the fit of a model that assumes no relationship between the species counts and the environment and no spatial autocorrelation. The "latlon" model represents the fit of a model that fits spatial autocorrelation, but no environmental effects. The "env" model represents the combined effects of all environmental predictors, and the remainder of the models ("temp", "oxygen", etc.) estimate the effects of a single predictor at a time.