Latitudinal diversity of planktonic copepods in the Eastern Pacific: overcoming sampling biases and predicting patterns

Rivera, Reinaldo; Escribano, Ruben; González, Carolina E.; Pérez-Aragón, Manuela

doi:10.3389/fevo.2024.1305916

ORIGINAL RESEARCH article

Front. Ecol. Evol., 24 April 2024

Sec. Biogeography and Macroecology

Volume 12 - 2024 | https://doi.org/10.3389/fevo.2024.1305916

Latitudinal diversity of planktonic copepods in the Eastern Pacific: overcoming sampling biases and predicting patterns

Reinaldo Rivera^1*

Ruben Escribano^1,2

Carolina E. González¹

Manuela Pérez-Aragón¹

¹Millennium Institute of Oceanography (IMO), University of Concepción, Concepción, Chile
²Department of Oceanography, Faculty of Natural and Oceanographic Sciences, University of Concepción, Concepción, Chile

Gradients of latitudinal diversity are one of the biogeographic features calling the most attention in ecology and macroecology. However, in pelagic communities of the marine environment, geographic trends and patterns are poorly known. We evaluated the latitudinal variation in species richness of marine planktonic copepods in the Eastern Pacific using spatial statistical models and approaches that mitigate and account for biases in occurrence data. A Boosted Regression Tree (BRT) and regression-Kriging based models allowed us to estimate and predict alpha diversity in poorly sampled regions, whereas beta diversity patterns were assessed using generalized dissimilarity analysis (GDM). Species richness showed a bimodal pattern, with a maximum of 291 species in the Northern Hemisphere and Tropical Eastern Pacific Ocean. Particulate Organic Carbon, salinity (max), spatial autocovariate, range of salinity and temperature, and Mixed Layer Depth, explained 85.2% of the latitudinal variability of copepods. Beta diversity was structured into four macrozones associated with the main water masses of the North and South Pacific.Our analytical approaches can overcome the limitations of data gaps, predicting greater diversity in subtropical and coastal areas, while providing insights into key drivers modulating spatial diversity patterns.

Introduction

The latitudinal diversity gradient (LDG) is one of the most striking features in ecology and biogeography. Contrary to the classic notion of an increase in diversity from the poles towards the equator (Rohde, 1992), studies in marine systems have shown that diversity is higher at mid-latitudes and decreases towards the equator (Rombouts et al., 2009; Barton et al., 2010; Chaudhary et al., 2016; Saeedi et al., 2017; Rivadeneira and Poore, 2020; Lin et al., 2021), which constitutes a pattern that has remained constant for millions of years (Yasuhara et al., 2012).

Over a large latitudinal gradient, the highest diversity of marine organisms has been reported in subtropical areas of the Northern Hemisphere, and a tropical-subtropical plateau in the Southern Hemisphere (Rutherford et al., 1999; Rombouts et al., 2009; Fautin et al., 2013; Chaudhary et al., 2016; O’Brien et al., 2016; Saeedi et al., 2017; Saeedi et al., 2019b; Rivadeneira and Poore, 2020; Chaudhary et al., 2021; Lin et al., 2021), together with a decrease in species richness near the equator (Chaudhary et al., 2016; Menegotto and Rangel, 2018). Although this pattern has been attributed to sampling biases (Menegotto and Rangel, 2018), empirical evidence supports this latitudinal trend (e.g., Chaudhary et al., 2016, 2017; Saeedi et al., 2017; Chaudhary and Costello, 2023, 2019b), reported for mollusks (Saeedi et al., 2017, 2019b), foraminifera (Brayard et al., 2005), coccolithophores (O’Brien et al., 2016), ophiuroids (Woolley et al., 2016), crustaceans (Rivadeneira and Poore, 2020), polychaete worms (Pamungkas and Glasby, 2021), brown macroalgae (Fragkopoulou et al., 2022), and higher trophic level organisms, such as fishes, sharks, squids, and cetaceans (Rutherford et al., 1999; Tittensor et al., 2010). This spatial pattern among taxonomically distant species suggests the existence of a common underlying mechanism that supports their similar distributions (Worm et al., 2005; Rombouts et al., 2011), which has been mainly attributed to the energy or temperature hypothesis (Hawkins et al., 2003; Rombouts et al., 2011).

The geographic patterns of diversity in marine environments have long been described (see Tittensor et al., 2010; Beaugrand et al., 2013; Miller et al., 2018); however, the ecological and evolutionary mechanisms of these patterns are not well understood (Roy et al., 1998; Gagné et al., 2020; Melo-Merino et al., 2020; Moreno et al., 2021). Most hypotheses attempting to explain contemporary patterns of diversity are based on the link between the abiotic environment and species diversity (Pianka, 1966; Currie, 1991; Rutherford et al., 1999; Kerr, 2001). Nevertheless, there is a lack of consensus regarding the mechanisms that cause such relationships (Stokes, 2018). Indeed, a particular set of environmental ecological hypotheses (e.g., energy, productivity, and environmental heterogeneity) has been suggested to explain the processes underlying the distribution of diversity in the oceans (e.g., Tittensor et al., 2010; Worm and Tittensor, 2018; Brandão et al., 2021). Among these factors, temperature and productivity may exhibit the highest explanatory power (Rutherford et al., 1999; Tittensor et al., 2010). Thus, they have been proposed as the main driver of diversity in multiple marine taxa (Irigoien et al., 2004; Brayard et al., 2005; Tittensor et al., 2010; Saeedi et al., 2019a; Lin et al., 2021).

In the same context, studies on planktonic copepods’ diversity and geographic variability in relation to environmental parameters have been performed (Rutherford et al., 1999; Rombouts et al., 2009, 2011; Hooff and Peterson, 2006; Medellín-Mora et al., 2016; González et al., 2020), showing a strong temperature dependence, resulting from the effect of this variable on the vital rates and ecological responses of copepods. Nevertheless, despite their ecological importance and the relatively high amount of data available on their distribution, open ocean studies on their biogeographic patterns are scarce and underrepresented in global biodiversity assessments, with studies mostly concentrated on local systems of the Pacific (e.g., Williamson and McGowan, 2010; Medellín-Mora et al., 2016; González et al., 2020, 2021). Furthermore, the mechanisms or factors that explain and modulate diversity on a wide scale are insufficiently understood, and too few studies have been conducted using explicit spatial analysis to understand and predict biogeographic patterns (Wilson et al., 2016; Brandão et al., 2021; Rivera et al., 2023).

When studying biogeographic patterns, the total diversity in a region (gamma diversity) can be decomposed into alpha and beta diversity. Alpha diversity refers to diversity within a local community, whereas beta diversity refers to variability between different communities or sites. It is relevant to note that the factors underlying alpha, or beta diversity are different (Leprieur et al., 2011), including historical and environmental processes. Therefore, it is important to explicitly know the factors that modulate the spatial structure of copepods along broad environmental gradients when evaluating ecological hypotheses (e.g., energy, heterogeneity, productivity) that would explain these diversity trends.

Beta-diversity is a less studied attribute of copepod diversity over a broad geographic scale. Beta diversity represents the extent of change in community composition or the degree of community differentiation in relation to environmental gradients (Whittaker, 1960). The importance of beta diversity lies in the fact that it not only accounts for the relationship between local and regional diversity, but also informs about the degree of differentiation among biological communities (Baselga, 2010). Therefore, knowing the beta diversity of copepods allows us to understand and evaluate the ecological drivers that modulate it, as well as and predict patterns in non-surveyed regions. At the ocean scale, biological sampling is often sparse and biased to certain regions; however, environmental data, such as satellite images, are often more readily available and cover a wide geographic space (Ferrier et al., 2007; McArthur et al., 2010), constituting useful inputs for estimating and predicting the geographic patterns of beta diversity through modeling. It is well-known that large areas of the open ocean remain unexplored or greatly under-sampled; thus prediction models may be considered highly valuable tools to overcome sampling gaps with low sampling coverage (Dornelas et al., 2018).

Planktonic copepods are a key basal component of the trophic web in marine ecosystems, linking primary production to all higher trophic levels. They have relatively short life cycles and respond quickly to oceanographic variations, such as fluctuations in temperature, oxygenation, stratification, primary production, and upwelling (Hays et al., 2005; Hooff and Peterson, 2006). Therefore, planktonic copepods constitute a suitable study model to assess biodiversity patterns in the ocean for animals with restricted migration and free-drifting assemblies. An important aspect in the study of spatial biodiversity patterns in copepods is that they can be conditioned by the quantity and gaps in the information (Hughes et al., 2021). Regarding the quantity, in recent decades, the availability of data on species occurrences in global databases has greatly improved, such as the Ocean Biodiversity Information System (OBIS) (Klein et al., 2019) together with the Geographic Information Systems (GIS) platforms, allowing us to evaluate spatial patterns and infer the underlying ecological processes. This progress also allows for the testing of multiple hypotheses explaining the geographical distribution of the richness and composition of planktonic organisms (Tittensor et al., 2010; Righetti et al., 2019). Regarding information bias, occurrences from field and online databases have temporal and geographic biases (Beck et al., 2014; Isaac and Pocock, 2015; Dornelas et al., 2018), which can mask the patterns of diversity (Hughes et al., 2021). In some cases, occurrences are directed toward certain groups of economic or ecological importance (taxonomic bias and societal preferences) (Titley et al., 2017; Troudet et al., 2017), proximity to research centers, or easily accessible areas (see Hortal et al., 2008; Ball-Damerow et al., 2019). Nevertheless, this is not an impediment to study large-scale diversity patterns given the development of approaches that allow the evaluation and consideration of biased information (Rocchini et al., 2011; Ruete, 2015; Monsarrat et al., 2019; Boyd et al., 2021; Tessarolo et al., 2021; Zizka et al., 2021; D’Antraccoli et al., 2022). These methodologies allow the estimation of diversity based on the wide spatial coverage of environmental predictors (Hortal and Lobo, 2011; Alves et al., 2020), allowing the elucidation of the modulating factors of copepod diversity.

The objectives of this study were: (i) evaluate the relationships between the latitudinal gradient of copepod diversity and multiple environmental variables, (ii) predict copepod richness in regions where sampling is currently deficient, and (iii) estimate and predict beta diversity patterns and their underlying environmental factors. Our analytical approaches aim to overcome biases and limitations derived from data gaps, while providing insights into key drivers modulating spatial patterns of diversity and composition for planktonic free-drifting organisms over a large scale in the ocean.

Materials and methods

Study area and occurrences

The study area corresponds to the Eastern Pacific (66.5°N to 60°S, and 180°W to 67°W). Georeferenced occurrences of copepods were obtained from the COPEPOD database (https://www.st.nmfs.noaa.gov/copepod) (O’Brien, 2014) and OBIS database (http://www.iobis.org) using the ‘robis’ package (Provoost and Bosch, 2022) and OBIS mapper application (https://mapper.obis.org). All online databases were accessed on January 10, 2022 (see Supplementary Table S1). Initial zooplankton observations gathered 81.417 occurrences. The search was limited to the 0–4000 m depth range, given the higher availability of occurrences for this layer in the water column, which constitutes the average depth of the oceans (3.682 m). Occurrences without geographic coordinates, coordinate duplicates, equal to zero, or located within continents were excluded.

The taxonomic validity of 692 species names was assessed by checking all names of the World Register of Marine Species (WoRMS, http://www.marinespecies.org) through the match_taxa function of the ‘robis’ package (Provoost and Bosch, 2022). Taxonomic assignment was performed at the species level, eliminating the occurrences at higher taxonomic levels. Only the species name classification was considered for occurrences classified as subspecies. Unaccepted species names have been corrected and synonyms and misspellings have been reconciled. After the data-screening and filtering procedures, 37.158 valid records of 525 species were retained for downstream analyses (see Supplementary Table S1).

Diversity of copepods

Species richness (SR) is the total number of different species present in a community, or specific geographic area. Unlike SR, taxonomic distinctiveness is a measure of biological diversity that is used to evaluate the uniqueness of species in a community or ecosystem. We thus calculated the average taxonomic distinctness (Δ+) and variation in taxonomic distinctness (Λ+) to estimate the average and expected deviation of diversity for the study area when compared with a species inventory. These indices show the average and expected deviations for the study area when compared to an inventory of species (aggregation matrix). For the region, Razouls et al. (2023) indicates 849 species (Zones 20, 24, 25 and 26) (https://copepodes.obs-banyuls.fr/en/searching.php). With the list of registered species, an aggregation matrix was built at the level of class, order, family, genus and species. The average ( $Δ^{+}$ ) (Equation 1) and variation in taxonomic diversity (Λ+) (Equation 2) are defined according to the following equations:

\begin{array}{l} Δ^{+} = 2 \frac{\sum \sum i < j w i j}{N (N - 1)} & (1) \end{array}

\begin{array}{l} Λ^{+} = 2 \frac{\sum \sum i < j {(ω i j - Δ^{+})}^{2}}{N (N - 1)} & (2) \end{array}

Where the double sum is applied to all species i, j. N is the sample number, and ωij is the assigned taxonomic path length of branch between i and j species (Clarke and Warwick, 1998). Λ+ is the variance of the taxonomic distances (ωij) between each pair of species i and j from their mean value (Δ+). The aggregation matrix was organized, and a classification tree was built from a random selection (n = 1000) of species from the general species inventory, thus comparing the indices and establishing a confidence interval of 95%. The TAXDTEST routine available in PRIMER 7 software was used (Clarke and Gorley, 2015).

To assess the spatial variability of the copepod diversity in the eastern Pacific, we created a map of the study area using a hexagonal cell. Each hexagonal cell covered an area of approximately 3 degrees. The species richness and rarefied species richness (ES50) were calculated for each hexagonal cell. We calculated the ES50 among 50 random samples, repeatedly sampled to standardize the data, and accounted for the sampling effort (Gotelli and Colwell, 2011; Saeedi et al., 2019a). It repeatedly resampled 50 randomly chosen records from all available records and calculated the average number of species per 50 records (Oksanen et al., 2022; Saeedi et al., 2022). It is important to highlight that if the data are not standardized by reference to the sampling effort, it is not possible to distinguish whether the patterns are biased or not (Fernandez and Marques, 2017). The ES50 estimates the total number of species in a community, including those which may have not been observed in the sample. This is based on the abundance of rare species present in the sample to infer the presence of non-observed species. It differs from rarefaction in that the latter standardizes the comparison of species richness between different samples, especially when the samples have different sizes. ES50 was calculated using the ‘vegan’ package (Oksanen et al., 2022).

To assess the modality in the number of species and the ES50, Hartigan’s dip statistic (HDS) (Hartigan and Hartigan, 1985) was used through the ‘diptest’ package (Maechler, 2021). Unimodality is assumed when p > 0.05 (i.e., not significant) and not unimodal if p < 0.05 (bi- or multimodal) (Hartigan and Hartigan, 1985).

Environmental database

Twenty-three environmental variables were selected to define candidate predictors for modeling the distribution of richness, ES50, and beta diversity. Chlorophyll-a concentration (chla-a) (maximum, minimum, and average), dissolved oxygenconcentration (maximum, minimum, and average), mixed layer depth (MLD), nitrate (mean), net primary productivity (NPP) (maximum, minimum, and average), photosynthetically active radiation (PAR), pH, phosphate concentration, particulate organic carbon (POC), sea surface salinity (SSS) (maximum, minimum, average, and range) and sea surface temperature (SST) (maximum, minimum, average, and range). These variables are recognized for their relationship with the physiology of copepods and for constraining their distribution (Escribano and Hidalgo, 2000). The oceanographic variables were obtained from Bio-ORACLE v 2.2 (Assis et al., 2018) at a spatial resolution of 5 arcminutes (0.08° or 9.2 km at the equator) (https://www.bio-oracle.org), spanning from 2000 to 2014. The variables PAR and POC were downloaded from the Global Marine Environment Dataset (GMED Version 2.0) https://gmed.auckland.ac.nz/index.html. The MLD was obtained from the Copernicus database (https://www.copernicus.eu) at a resolution of 0.25 degree. All the environmental variables were aggregated at a spatial resolution of 0.08 degree. To obtain the average value and range of the environmental variables for each 3° cell, zonal statistics were used. This method allows calculating descriptive statistics of the data that are associated with specific zones and how these vary within different geographic units (i.e., 3° cells).

Variable selection

Multicollinearity between environmental predictors was tested using the Variance Inflation Factor (VIF) and Spearman’s rank correlation coefficient. Highly correlated variables were considered if VIF > 3 and ρ > 0.85 (Dormann et al., 2013). This procedure ensures the removal of collinear descriptors and reduces overfitting (Guillaumot et al., 2020). Based on the low correlation and available empirical evidence (e.g., Escribano and Hidalgo, 2000; Rombouts et al., 2009), eight predictors were retained and used to develop the Regression-Kriging and Generalized Dissimilarity Modelling (described below): MLD, NPP (max), NPP (min), pH, POC, SSS (max), SSS (range) and SST (range) (Supplementary Table S2; Supplementary Figure S1). VIF analysis was performed using the ‘usdm’ package (Naimi et al., 2014), and correlation analyses were performed using the ‘raster’ (Hijmans, 2023) and ‘corrplot’ (Wei and Simko, 2021) packages.

Relationships between species richness and environmental variables

To evaluate the relationships between oceanographic variables and the number of species, as well as ES50, we used a Boosted Regression Tree (BRT) (Elith et al., 2008). Exploratory analysis highlighted that species richness and ES50 data have non-linear relationships with the explanatory variables; therefore, we used BRT because it deals with non-linear and non-monotonic relationships between response and explanatory variables (Elith et al., 2008). The BRT parameters were selected to optimize the model using the ‘caret’ (Kuhn, 2008) and ‘dismo’ packages (Hijmans et al., 2021). The tree complexity of the model was fixed at 2, and the learning rate was 0.01. Although BRT is a flexible technique and has shown high performance relative to other algorithms (see Elith et al., 2008), it does not explicitly consider spatial autocorrelation. Therefore, we used the residual autocovariate (RAC) technique to address spatial autocorrelation. Spatial autocorrelation was accounted for by adding an autocovariate term to the final selected BRT. The RAC represents the influence of neighboring observations on the response variable at a specific location (Escalle et al., 2016). First, the BRT model was computed using uncorrelated variables (Supplementary Figure S1 and Supplementary Table S2). The residuals calculated for each grid cell were used to compute the autocovariate using focal calculation. The residual autocovariate was considered an explanatory variable in the BRT model (Crase et al., 2012). Second, the best fitting model was selected as the model with the lowest Akaike Information Criterion (Burnham and Anderson, 2004) and the highest deviance was explained. The Moran index was calculated to evaluate the autocorrelation of the residuals of BRT without RAC (simple BRT) and BRT+RAC models using the ‘spdep’ package (Bivand, 2022). The models were estimated using a Poisson distribution, and model selection was performed using forward and backward approaches considering the total explained deviance. The BRT models were performed using the ‘gbm’ package (Greenwell et al., 2022).

Estimation of species richness

To predict the geographic gradient of species richness, regression-kriging (R-K) was performed (see Alves et al., 2020). This method considers the richness-environment relationship and its autocorrelation structure to predict the richness gradients (Alves et al., 2020). First, well-sampled cells were selected to calculate the semi-variogram. A fitted species accumulation curve was calculated using a rational function to describe the association between species richness and sampling effort for each cell (Lobo et al., 2018). Following the proposal of Alves et al. (2020), three sampling effort metrics were combined: a) completeness (i.e., percentage of species richness estimated by the accumulation curve), b) the ratio between the number of occurrences within the cell and the observed species richness, and c) the slope between the accumulated number of species and the number of occurrences (Lobo et al., 2018). First, the values to establish well-sampled cells were selected according to the values proposed by Alves et al. (2020): completeness > 95%, ratio > 15, slope < 0.02, and minimum number of occurrences (10). The number of well-sampled cells was estimated by using the ‘KnowBR’ package (Lobo et al., 2018). Secondly, the semi-variance values of the variogram were fitted with a spherical function because it presented the lowest root mean squared error (RMSE) (see Supplementary Table S3). Generalized Linear Models (GLM) were fitted only from well-sampled cells to identify the significant predictors. A model selection approach was performed to retain the subset of variables that best explained species richness, and these variables were predicted by R-K (Figure S2). Finally, species richness was predicted using ‘automap’ (Hiemstra et al., 2009) and ‘gstat’ packages (Pebesma, 2004; Pebesma and Graeler, 2021).

Geographic pattern of beta diversity and clustering of species composition

To model and evaluate the spatial pattern of beta diversity, we used Generalized Dissimilarity Modelling (GDM). This method relates biological distance to ecological distance, and estimates and predicts spatial patterns of turnover in community composition across large areas (Ferrier et al., 2007; Brown et al., 2014). GDM is an extension of the GLM, allowing the modeling of community dissimilarity against a set of environmental predictors (Ferrier et al., 2007). GDM estimates the magnitude and rate of turnover along environmental gradients, representing the dissimilarity between pairs of sampling units as a function of environmental differences and geographical distances (Fitzpatrick et al., 2013). GDM was fitted to the presence of species at each site, and compositional dissimilarity was predicted at unsampled localities based on environmental data. Eight environmental covariates were used: MLD, NPP (max), NPP (min), pH, POC, SSS (max), SSS (range), and SST (range) (Supplementary Figure S2; Supplementary Table S2). The possible influence of spatial processes was evaluated by including the geographic distance between pairs of locations as a predictor variable. The importance of the variables was assessed by a permutation test (n = 100). The analyzes were performed with the ‘gdm’ package (Fitzpatrick et al., 2020).

To determine regions with similar species composition (GDM analysis), we performed cluster analysis using the unsupervised classification algorithm K-means (Lloyd’s algorithm) to partition the data into a set of K clusters. To estimate the optimal number of clusters, we used the Silhouette method with the ‘factoextra’ package (Kassambara and Mundt, 2020), and a map of this classification was created using the ‘ecbtools’ package (Williamson, 2016).

Finally, we followed best-practice standards in species distribution modeling, providing a full description of the modeling steps following the ODMAP (Overview, Data, Model, Assessment and Prediction) protocol (Zurell et al., 2020) (Supplementary Table S4). All statistical analyzes were performed in the R software (R Core Team, 2023) and the environmental variables were processed in ArcGIS 10.4.1 (ESRI, 2016).

Results

Species richness

Based on the inventory of 692 species of copepods described in the Eastern Pacific and the list of 525 species recorded in this study, the mean (Δ+) and variation taxonomic distinctness (Λ+) showed that the values of taxonomic distinctness were mostly within the 95% confidence level, indicating that the observed copepod diversity contained the expected level of diversity for the region (Funnel plot, Figure 1). The Λ + had a mean of 69 and an expected variance of 61.5 (Figure 1).

Figure 1

Figure 1 (A) Average taxonomic distinction index (Δ+) and (B) variation in taxonomic distinctness’ (Λ+) for the copepod community in the Eastern Pacific. The segmented line in red shows the mean, whereas the continuous lines are the distribution of probability at 95%.

Relationship of number of species, ES50 and environmental variables

The BRT+RAC model had the best fit to explain the number of species and retained six variables. The percentage of total deviance explained was approximately 85.2%. POC accounted for most of the relative influence (21.4%), followed by SSS, spatial autocovariate, SSS (range), SST (range), and MLD (Table 1). Because spatial autocorrelation was detected in the BRT model without autocovariates (Table 1), we chose the BRT+RAC model given its lack of spatial autocorrelation (Moran’s index = 0.004; p value = 0.485) (Table 1).

Table 1

Table 1 Summary of the BRT Poisson and BRT + residual autocovariate (RAC) used to model number of species and estimated species richness (ES50) of copepods.

The BRT+RAC model had the best fit to explain the ES50 and retained six variables. The percentage of total deviance explained was approximately 81.3%. The autocovariate accounted for most of the relative influence (25.7%), followed by SST (range) (18.2%), SSS (max), MLD, POC, and SSS (range) (Table 1 and Supplementary Figure S3).

Observed and predicted latitudinal diversity gradients

The observed richness showed a heterogeneous pattern in the Eastern Pacific, with the highest values observed (291 species) off the coast of Alaska, Ecuador, northern Chile, and around 30°S–179°W (Figure 2A). ES50 showed a spatial pattern similar to that of the number of species, revealing the greatest richness at intermediate latitudes, with higher diversity in coastal areas in the Northern Hemisphere and in oceanic areas in the Southern Hemisphere (Figure 2B). The relationship of species richness and ES50 with latitude was significantly different from a unimodal distribution, indicating bimodality in both trends (Hartigan’s dip test, D = 0.048 and and D = 0.045; p-value < 0.01 and < 0.05, Figures 3A, B, respectively). Well-sampled cells were observed primarily at the Bering Sea, Gulf of Alaska, Gulf of Mexico, Eastern Tropical Pacific Ocean (ETPO), and the southern hemisphere between 40–58°S and 180°W (Supplementary Figures S4 and S5). The variogram model that best fit the pattern of species richness among the well-sampled cells was the Spherical function (sill = 1489; range = 588; nugget = 62; kappa = 1.7; Supplementary Table S3). The diversity predicted by R-K revealed higher predicted richness, with a maximum of 362 species (Figure 2C). The highest species richness was concentrated in continental margins, mainly in tropical and subtropical areas, and decreased in temperate areas. The difference between the observed and estimated richness showed extensive poorly sampled regions (Supplementary Figure S6), and higher residual variation was detected in the areas of low observed richness (Figure 2D). The species richness estimated by R-K showed a bimodal pattern (Hartigan’s dip test, D = 0.059; p-value < 0.05) (see Supplementary Figure S7).

Figure 2

Figure 2 (A) Spatial distribution of number of species, (B) estimated species richness (ES50), (C) richness estimated from Regression-Kriging, and (D) residuals. Isolines represent richness predicted by R-K. Datum=WGS84. The map was constructed with ArcGIS 10.4.1 (ESRI, 2016). Map uses Behrmann equal area projection.

Figure 3

Figure 3 (A) Bimodality of number of species and (B) estimated species richness (ES50) of Eastern Pacific copepods. Black dots = empirical species richness; blue line = estimated trend; gray shade represents 95% point wise confidence interval. The figure was constructed with R software (R Core Team, 2023).

Spatial patterns of beta diversity

GDM explained 53.8% of the deviance in beta diversity. The predictors with the best fit to explain this pattern were geographical distance, SST (range), POC, SSS (max), NPP (min), MLD, and SST (range) (Supplementary Figure S8). Spatial prediction of the variation in copepod composition revealed regions with similar copepod communities distributed along a latitudinal gradient (Figure 4A). The composition differs between the temperate and equatorial zones, as well as between oceanic and continental areas (Figure 4A). Cluster analysis revealed four areas with similar species composition. The areas were numbered latitudinally from C1 to C4 for better interpretation (see Figure 4B). There is an exclusive group for the Northern Hemisphere in temperate and cold zones (C1) covering the Bering Sea, Gulf of Alaska, and part of the California Current, and another exclusive group for the Southern Hemisphere (C4) in the Southern Ocean (Figure 4B). The C2 cluster is shared between both hemispheres and is mainly oceanic, whereas the C3 cluster is oceanic-coastal with a large latitudinal distribution, encompassing the Humboldt Current System and the Pacific Central-American coast (Figure 4B).

Figure 4

Figure 4 (A) Predicted spatial variation in copepods species composition. Locations with similar colors are expected to contain similar communities. (B) Map illustrating the clustering of GDM. A four-class classification of the region derived from predicted dissimilarities (C1 to C4). Symbol colors represent each plot’s membership in a specific assemblage. The figure was constructed with ArcGIS 10.4.1 (ESRI, 2016). Map uses Behrmann equal area projection.

Discussion

The geographic tendency of the specific richness of copepods in the Eastern Pacific is concordant with the bimodal pattern previously described for other taxa (e.g., Chaudhary et al., 2016; Saeedi et al., 2017; Rivadeneira and Poore, 2020; Pamungkas and Glasby, 2021), showing a significant asymmetric and bimodal pattern with higher diversity in the Northern Hemisphere (Figures 2A, B, and Supplementary Figure S7). From an ecological perspective, the bimodal pattern reported for copepods has the potential to be explained by SST, bathymetry, and NPP. Regarding temperature, this could exclude species intolerant to high values that are typical in tropical areas (Chaudhary and Costello, 2023). Regarding bathymetry, a wider continental margin in the Northern Hemisphere would provide a larger habitat, which is concordant with the species-area hypothesis (Lomolino, 2000). This is complemented by high primary productivity, which is characteristic of upwelling zones in the North Pacific (Zaytsev et al., 2003) and in the South Pacific off the coast of Chile and Peru (Escribano et al., 2012). Other mechanisms that have not been explicitly evaluated to date, such as physiological limitations, biotic interactions, and trophic resource availability, could also explain the dip in species richness in the equatorial zones (Chaudhary et al., 2016; Chaudhary and Costello, 2023). Therefore, the bimodality and asymmetry of species richness indicate that there is not only a general factor that determines diversity patterns on a broad scale. These tendencies could be investigated in regional studies to evaluate the effect of biotic interactions (e.g., co-occurrence) as a potential causal mechanism for bimodality.

According to the BRT+RAC model, the most important variables that explained and modulated the number of species were POC, SSS (max), and autocovariate (Table 1). Although there is no single predictor of species richness, there are fundamental variables, such as temperature, that modulate the large-scale patterns of pelagic diversity (Rombouts et al., 2009; Tittensor et al., 2010). Based on to Tittensor et al. (2010), SST alone is one of the most efficient predictors of coastal and oceanic species distribution on a global scale. Our results revealed that variability of temperature, and not its average value, was an efficient predictor of richness along the latitudinal gradient. Therefore, high temperature variability, together with the presence of warm waters in the equator, could constitute the main cause of the observed bimodal gradient because equatorial regions are too warm for some species (Chaudhary et al., 2021). Another physical mechanism that could explain the difference between areas of high diversity in the subtropics and less diversity in the equatorial zone would be the presence of two opposing currents (Equatorial Counter Current and South Equatorial Current). These two currents have distinct hydrological features: the warm pool to the west and the Pacific Equatorial Divergence, regions separated by well-defined fronts in salinity, pC02, and macronutrients (Le Borgne et al., 2002). These physical characteristics would generate a barrier such as the East Pacific Barrier reported for tropical fish (Gaither et al., 2016), and plankton (Bowen et al., 2016), which would limit the spread of organisms with a consequent break in diversity.

Regarding SSS as a modulator of species richness, it has a direct influence on feeding and egg production (Dutz and Christensen, 2018), growth rate, and population composition (Magouz et al., 2021), which are factors that ecologically regulate the community structure of the systems. Another modulator of copepod diversity is NPP, which has a positive relationship with richness, and a reported trend for marine and terrestrial organisms (Mittelbach et al., 2001; Chase and Leibold, 2002; Woodd-Walker et al., 2002; Fukami and Morin, 2003; Rombouts et al., 2009). This relationship can be explained by the fact that productivity is a key factor in the composition of copepod diets (Vargas et al., 2006, 2010). This is particularly important in coastal upwelling systems, given that high-quality food resources and lower temperature variability exert an advantageous effect on copepods, despite constituting an advective environment that may generate important changes in their life cycles (Peterson, 1998), biomass, and abundance (Escribano et al., 2012; Medellín-Mora et al., 2016) with potential consequences for the ecosystem food web.

Another important factor spatially structuring the observed and estimated diversity of copepods is the MLD, which indicates the relevant (effective) temperature of the upper-ocean copepods at a macroscale (Rombouts et al., 2009, 2011; Rajakaruna and Lewis, 2018). Indeed, a strong negative relationship was observed between diversity and MLD (Supplementary Figure S3) indicating that diversity was high in stratified waters, which may be due to that a stable physical structure in the near-surface ocean provides higher vertical niche availability than can support a high number of species (Rutherford et al., 1999). Indeed, copepod diversity was higher in areas with stable seasonal ocean temperatures, trends that have also been reported for the Atlantic Ocean (Woodd-Walker et al., 2002). This result corresponds with the observation that oligotrophic regions, with low seasonal variability have high copepod diversity (Medellín-Mora et al., 2021).

The high contribution of POC to the spatial variability of copepod richness is related to the fact that zooplankton in general, are important mediators of the flux of organic material into the deep ocean through contributions to both the active and passive carbon fluxes (the biological pump) (Ducklow et al., 2001), and their community composition has a significant influence on particle export, repackaging of POC with depth, carbon attenuation, and export through the mesopelagic zone (Steinberg et al., 2008; Wilson et al., 2008). In this regard, Omand et al. (2015) pointed out that disaggregation of large POC by zooplankton and the detrainment of POC-rich water from the mixed layer via the so-called mixed-layer pump (Gardner et al., 1995) constitutes one of the mechanisms that explain the export of small particles (0.2 to 20 μm) to depths up to 1000 m (Dall’Olmo and Mork, 2014). The discrepancies in POC, SST and SSS values in the Pacific Ocean are influenced by a complex interaction of various factors (Gebbie and Huybers, 2012). These factors include different water masses present in the region and the circulation patterns of the thermohaline current (Gebbie and Huybers, 2012). These water masses originate from different geographical areas and exhibit a variety of temperatures, salinity levels and ages, which interact and intermingle. The thermohaline current plays a fundamental role in this process, as it facilitates the transport of warm and saline waters towards higher latitudes, while directing cold and less saline waters towards lower latitudes. This causes significant mixing in the subtropical and tropical regions of the ocean, and therefore could constitute a macroscale structuring of the unequal distribution of species richness along the latitudinal gradient (Emery, 2001).

Our species richness prediction allowed us to identify areas of high diversity on the coasts of the Gulf of Alaska, California Current, Eastern Tropical Pacific Ocean (ETPO), and central-northern Chile (Figure 2B). The reported pattern can represent the differential sampling effort made in both hemispheres, with gaps and biases in the spatial distribution of the occurrences (Hickisch et al., 2019; Maitner et al., 2023). To overcome these limitations, the R-K method based on well-sampled cells (see Supplementary Figure S4) was used to fill these gaps (Lobo et al., 2018). Robustness in the prediction of the variation on richness by R-K at a broad geographic scale is supported by considering the autocorrelation intrinsically present in species richness (Legendre, 1993; Dormann, 2007; Snickars et al., 2014; Alves et al., 2020), allowing for a better understanding of the drivers of its geographic variation (Diniz et al., 2003; Dormann, 2007). Similarly, the latitudinal gradient of beta diversity showed changes from the coastal to the open ocean, as well as between hemispheres (see Figure 4A). For the East Pacific, the four clusters (Figure 4B) represent unique communities, suggesting high idiosyncratic and dispersal limitations in the East Pacific. The spatial tendency of the composition revealed an area exclusive to the Northern Hemisphere and one in the Southern Hemisphere, and two shared clusters in the temperate and southern regions (Figure 4B). This indicates that although there is a bimodal gradient of richness, the specific composition has a defined structure that obeys the dynamics of the current systems in both hemispheres (e.g., California Current, Humboldt Current System, Pacific Equatorial Countercurrent and North/South Pacific Subtropical Gyres) with a unique fauna with high levels of endemism as a result of basin-scale circulation patterns and mesoscale eddies, which come from the coastal upwelling zone, generating plankton transport (Medellín-Mora et al., 2021).

GDM revealed that the turnover in community composition changed at a greater rate with increasing geographical distance, SSS (range), and POC, showing a broad latitudinal gradient in turnover from Alaska to Tierra del Fuego (Figure 4A). The relationship was approximately linear for SST (range), indicating a weak effect on compositional turnover along these gradients (Supplementary Figure S9), for which processes related to oceanographic dynamics, such as POC circulation, productivity, and MLD, were the main modulators of replacement patterns. The productivity gradients are strongly influenced by local winds, bathymetry, and upwelling patterns (Sobarzo et al., 2007; Medellín-Mora et al., 2021), which generates heterogeneity in the water column chemistry. These environmental gradients can generate shifts in species composition and diversity patterns, as has been reported in some groups, such as phytoplankton on a global scale (Righetti et al., 2019), as well as in other regions, such as the South Australian Gulfs and semi-enclosed bays in Victoria (Leaper et al., 2011). An aspect of interest in future research would be to evaluate the patterns of community composition in upwelling and non-upwelling regions to study regional processes that may affect species turnover at the mesoscale.

To our knowledge, quantitative analysis and prediction of beta diversity have not yet been conducted in copepods of the East Pacific (but see Chaudhary and Costello, 2023 for a different approach). Knowledge of these geographic trends contributes to a better understanding of species shifts and similarities across broad environmental gradients (Buckley and Jetz, 2008; Jankowski et al., 2009). The geographic patterns of the composition of copepods studied through GDM constitute the first explicit analysis of a large dataset to geographically delimit regions confirmed by changes in copepod communities.

A critical aspect in the study of spatial patterns of diversity is the limitation that may be generated by spatial and temporal biases in the occurrence of species (Boakes et al., 2010; Bowler et al., 2022). The availability of occurrences on a global scale through OBIS is valuable for the study of marine macroecological patterns; that together with an appropriate curation and utilization of suitable statistical tools that consider autocorrelation (e.g., spatial regression, geostatistics) allows the identification of diversity trends on a broad geographic scale (see Alves et al., 2020; Cavalcante et al., 2022). Indeed, spatial prediction techniques, such as R-K and GDM, allow the generation of a robust prediction in regions that are less studied, enabling the overcoming of shortfalls in knowledge that limit our understanding of diversity patterns and gradients (see Hortal and Lobo, 2011; Alves et al., 2020).

Finally, knowledge on the main ecological factors driving diversity patterns (richness and composition) of copepods constitutes an important opportunity for understanding ecological processes that operate on a broad geographic scale and allows the identification of areas where greater sampling efforts are required.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

RR: Conceptualization, Data curation, Formal analysis, Methodology, Writing – original draft, Writing – review & editing. RE: Conceptualization, Writing – original draft, Writing – review & editing. CG: Formal analysis, Writing – review & editing. MP-A: Conceptualization, Formal analysis, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was funded by Grant the Millennium Institute of Oceanography (Grant ICN12_019N.

Acknowledgments

RR thanks Magdalena Sanhueza for her constructive and critical comments on the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2024.1305916/full#supplementary-material.

References

Alves D. M. C. C., Eduardo A. A., da Silva Oliveira E. V., Villalobos F., Dobrovolski R., Pereira T. C., et al. (2020). Unveiling geographical gradients of species richness from scant occurrence data. Glob. Ecol. Biogeogr. 29, 748–759. doi: 10.1111/geb.13055