Patterns in Bathyal Demersal Biodiversity and Community Composition Around Archipelagos in the Tropical Eastern Pacific

Giddens, Jonatha; Goodell, Whitney; Friedlander, Alan; Salinas-de-León, Pelayo; Shepard, Charles; Henning, Brad; Berkenpas, Eric; Sala, Enric; Turchik, Alan

doi:10.3389/fmars.2019.00388

ORIGINAL RESEARCH article

Front. Mar. Sci., 09 July 2019

Sec. Deep-Sea Environments and Ecology

Volume 6 - 2019 | https://doi.org/10.3389/fmars.2019.00388

This article is part of the Research Topic Pacific Deep-Sea Discoveries: Geological and Biological Exploration, Patterns, and Processes View all 24 articles

Patterns in Bathyal Demersal Biodiversity and Community Composition Around Archipelagos in the Tropical Eastern Pacific

$\r\nJonatha Giddens,*$ Jonatha Giddens^1,2*

Whitney Goodell^1,2,3

Alan Friedlander^2,3

Pelayo Salinas-de-León^3,4

Charles Shepard¹

Brad Henning¹

Eric Berkenpas¹

Enric Sala³

Alan Turchik¹

¹Exploration Technology Lab, National Geographic Society, Washington, DC, United States
²Fisheries Ecology Research Laboratory, Hawaii Institute of Marine Biology, Kaneohe, HI, United States
³Pristine Seas, National Geographic Society, Washington, DC, United States
⁴Charles Darwin Research Station, Puerto Ayora, Ecuador

Deep-sea fauna of the Tropical Eastern Pacific (TEP) have remained largely undescribed because exploration has historically been a challenge in these remote locations. Consequently, little is currently known about deep-ocean biodiversity in the TEP. An enriched understanding of biogeographic patterns and the factors that influence them is crucial to the success of efforts to conserve deep-sea communities. The objectives of this research were to provide the first assessment of biodiversity and community composition of bathyal demersal (200 – 3,000 m) fauna across multiple ocean/island archipelagos in the TEP. To better understand biodiversity and its drivers in the region, we tested the hypotheses that biodiversity and faunal assemblage composition vary with location and environmental drivers (e.g., bathymetry, oxygen concentration). Using Deep-Ocean Dropcams, we collected video footage from Clipperton, Galápagos, Revillagigedo, and Malpelo oceanic islands, in depths ranging from 200 to 3,000 m. Taxa richness, Shannon diversity, and Simpson diversity were calculated for each location. In general, biodiversity values were high closer to the equator and closer to the continent, and low in the locations farther north and farther west. Fish family biodiversity was highest in the Galápagos, followed by Malpelo. In contrast, Clipperton had significantly lower fish family richness compared to all other locations. For invertebrate orders, Revillagigedo was significantly lower in biodiversity compared to the Galápagos and Malpelo. Broad-scale Bathymetric Position Index (BPI) and dissolved oxygen had a positive effect on fish diversity metrics. Both fish and invertebrate assemblages were significantly different among locations. Around each of the island groups are vastly different current patterns, which likely accounts for some of the variation we observed in biodiversity and community composition. Due to the longitudinal East Pacific Rise, which limits population connectivity by creating a spatial barrier to gene flow and dispersal among marine organisms, the TEP is largely isolated from the wider Pacific. The islands and archipelagos likely represent important stepping-stones and corridors for species connectivity across the TEP and the wider Pacific. Therefore, the preservation of these intact ecosystems is crucial for maintaining biodiversity both within the TEP as well as across the wider Pacific region.

Introduction

The deep-sea is a severely under-sampled ecosystem and much remains to be understood about deep-sea biogeography (reviewed in Wilson and Kaufmann, 1987; Rogers, 1994). An enriched understanding of biogeographic patterns and the factors that influence them is crucial to the success of efforts to conserve deep-sea communities, especially in remote islands and archipelagos in the Tropical Eastern Pacific (TEP).

The TEP is separated from the central and western Pacific by the East Pacific Rise, a divergent tectonic plate boundary on the Pacific floor running north to south that coincides with the boundaries of a unique biogeographic province with a high proportion of endemic species (Ekman, 1953; Allen, 2008; Figure 1). The TEP can be divided into three biogeographic provinces based on shore fish (0 – 100 m) assemblage similarities (Robertson and Cramer, 2009). The continental coast contains two provinces: the Cortez (Gulf of California and lower Pacific Baja) and the Panamic (southward province). These two distinct areas are identified by a peak in abundance of local endemics and overall species richness. The third province within the TEP consists of five isolated ocean islands and archipelagos: the Galápagos, Revillagigedo, Clipperton, Cocos, and Malpelo; collectively termed the “Ocean Island Province” (Hubbs and Briggs, 1974; Robertson and Cramer, 2009).

FIGURE 1

Figure 1. Four study locations in the Tropical Eastern Pacific (TEP), and the associated currents influencing each site, with mean sea surface temperature (World Ocean Atlas) (adapted from Fiedler and Talley, 2006; Kessler, 2006), areas of severe hypoxia (<0.2 ml l⁻¹; adapted from Diaz and Rosenberg, 1995), axis of the East Pacific Rise, and biogeographic provinces (adapted from Robertson and Cramer, 2009).

Small ocean islands, such as those in the TEP, typically support faunal assemblages of lower species richness compared to adjacent continental areas. The fauna of these oceanic islands also have different functional-group composition, including both more broad-ranging species and more localized endemics. Thus, this area forms a hotspot for both the origin and maintenance of biodiversity in the sea (Robertson and Cramer, 2009).

The smaller size of island assemblages compared with neighboring mainland areas makes them more sensitive to extinctions (Dulvy et al., 2003; Hobbs et al., 2011). Due to the longitudinal East Pacific Rise, which limits population connectivity by creating a spatial barrier to gene flow and dispersal among marine organisms, the TEP is largely isolated from the wider Pacific (Lessios and Robertson, 2006). The islands and archipelagos likely represent important stepping-stones and corridors for species connectivity across the TEP and the wider Pacific (Glynn et al., 1996).

Previous research has assessed marine biodiversity at these islands, but has focused largely on the shallow-water ecosystems (Allen and Robertson, 1994; Glynn et al., 1996; Jost and Andréfouët, 2006; Alzate et al., 2014; Fourriére et al., 2014). Overall, the TEP is noted for its high level of endemism (72%) and relatively diverse shore-fish assemblages (Zapata and Ross Robertson, 2007). Additionally, the eutrophic mesopelagic zone (200 m – 1000 m), hosts a large proportion of endemic deep-water fish species adapted to the low oxygen in these waters (Priede, 2017).

However, the bathyal demersal (200 m – 3000 m) fauna of the TEP have remained largely undescribed because exploration has historically been a challenge in these remote locations (but see Auster et al., 2016; Lubetkin et al., 2018; Sánchez-Jiménez et al., 2018). Slopes and seamounts around Isla del Coco (Costa Rica) have been surveyed to depths of 450 m using manned submersibles (Cortés and Blum, 2008), recent explorations at the Galápagos used Remotely Operated Vehicles and manned submersibles to depths of 4000 m (Carey et al., 2016). Besides these localized studies, little is currently known about biodiversity in deep waters across multiple oceanic islands in the TEP.

The objectives of this research were to characterize biodiversity and community composition of bathyal (200 m – 3000 m) demersal fauna at four ocean/island archipelagos to better understand deep-ocean biogeography and its drivers in the TEP. We tested the hypotheses that biodiversity and faunal assemblage composition vary with location and environmental variables such as bathymetry and oxygen concentration. Our study represents the first characterization of deep-sea biodiversity across multiple islands and archipelagos in the region.

Materials and Methods

Study Region

The TEP is characterized oceanographically by extensive stratification with a sharp permanent pycnocline that prevents local ventilation of subsurface waters, a shallow thermocline (<80 m), and very low oxygen (<50 μmol⋅kg⁻¹ in the largest extent of the oxygen minimum zone at ∼ 10 degrees north) (Karstensen et al., 2008; Robertson and Cramer, 2009; Fiedler and Lavín, 2017; Priede, 2017). The boundaries of this biogeographic region are located northward to Magdalena Bay, Baja California (∼25°N) and southward to the Gulf of Guayaquil’s southern shore (∼4°S). These northern and southern boundaries are defined not by hard physical barriers, but by sharp water temperature gradients between tropical and temperate conditions (Figure 1).

The TEP is a region of complex oceanographic patterns. These waters lie between the two large subtropical gyres of the North and South Pacific, and at the convergence of two eastern boundary currents, coming from the north (California Current) and the south (Peru, or Humboldt, Current) (Kessler, 2006). Revillagigedo sits in the North Equatorial Current (NEC) and at the terminus of the colder waters of the California Current along the eastern boundary from the north (Figure 1). Clipperton, located ∼1,000 km to the south of Revillagigedo, sits in the East Pacific Warm Pool, a region of warm surface waters between the NEC and the North Equatorial Countercurrent (NECC) (Wang and Enfield, 2001). Moving south and eastward, the Galápagos Islands straddle the equator, and are influenced by the cold tongue of the South Equatorial Current, and the Peru (or Humboldt) Current, which brings cold waters along the eastern boundary from the south. Malpelo Island, furthest to the east and five degrees north of the equator, sits in the eastern curve of the northern gyre (NEC and NECC).

Data Collection

National Geographic’s Exploration Technology Lab developed Deep-Ocean Dropcams to observe deep-sea life in situ by capturing high quality imagery of the sea floor (Turchik et al., 2015). Deep-Ocean Dropcams are autonomous (free-falling) lander instruments, with high definition cameras (Sony Handycam HDR-XR520V 12 megapixel) encased in a 43-cm diameter borosilicate glass sphere that are rated to 10,000 m depth. We also deployed a Dropcam Mini, encased in a 33-cm diameter borosilicate glass sphere and rated to 7,000 m. This Dropcam Mini housed a Sony Handycam FDR-AX33 4K Ultra-High Definition video with a 20.6 megapixel still image capability. Viewing area per frame for both cameras was between 2 and 6 m², depending on the steepness of the slope where the Dropcam landed. Cameras were baited with ∼1 kg of previously frozen fish (bonito, mackerel, and tuna) and deployed for 6 to 9 h. Over the deployment period, cameras were programmed to record in 30–60 min intervals.

Biological Data

Using the Dropcams, we collected video footage from four ocean island/archipelago sites in the TEP: Clipperton, Galápagos, Revillagigedo, and Malpelo, in depths ranging from 200 m to ∼3000 m (n = 9–34 drops per site). From the video footage, fauna were classified to the lowest possible taxonomic level. Invertebrates were typically identified to taxonomic levels of order or class, whereas fish were identified to family level in most cases. The relative abundance of each identifiable taxon was calculated as the maximum number of individuals per frame (MaxN). For presence/absence analyses, these observations were converted to incidence per deployment.

Geospatial Data

Environmental predictor data were derived from publicly available geospatial databases. Variables examined included environmental/oceanographic data [e.g., Particulate Organic Carbon (POC) flux, and dissolved oxygen], geographic data (e.g., distance to geomorphic features), and bathymetric/habitat feature data [e.g., slope-of-slope (second derivative of slope), bathymetric position index (BPI)] (Supplementary Table S1). Data were extracted for each Dropcam deployment, based on the deployment’s lat/long position. For some predictor variable datasets, data were extracted directly from the source data pixel underlying the deployment point. For others, data were derived from source data via geoprocessing steps, such as calculating a maximum value of a variable within a given radius of the deployment point. All geoprocessing was executed in ArcMap 10.2 (ESRI, 2011). Here, we describe only the variables included in the final model, but see Supplementary Table S1 for more details on all geospatial variables examined.

Oxygen availability is a limiting factor for deep-sea faunal distributions (Levin, 2003; Stramma et al., 2010). Data for dissolved oxygen (annual mean) for defined depth bins were derived from the National Oceanographic Data Center (NODC) World Ocean Atlas (WOA05) and resampled by Watling et al. (2013) to 1 km × 1 km grids for defined depth bins. Dissolved oxygen values were extracted from the corresponding depth layer from Watling et al. (2013) associated with each deployment’s depth (Supplementary Table S1).

BPI, an index which highlights areas of topographic anomalies (Jones and Brewer, 2012), is a continuous metric of the elevation of a pixel relative to other pixels in a landscape (Lundblad et al., 2006), used to identify benthic features (e.g., ridges and valleys). Positive BPI values indicate features that are higher than the surrounding area, such as peaks and ridges, while negative cell values indicate features such as troughs. BPI values at and near zero are either flat areas or areas of constant slope, while sharp peaks or pits are represented by higher absolute values (Lundblad et al., 2006; Wright et al., 2012). BPI was calculated in ArcMap 10.2 using the Benthic Terrain Modeler extension (Wright et al., 2012). BPI can be calculated at broad or fine scales, with the scale of detected features dependent on the scale of BPI. Both broad- and fine-scale BPI have been shown to be associated with variations in marine faunal communities (Dolan et al., 2008; Jones and Brewer, 2012; Robert et al., 2015), including those of the deep-ocean (Leitner et al., 2017). Chosen BPI scales reflected the resolution of the available bathymetric and geomorphic data (100 m and ∼1 km, respectively) and the scale of the geomorphic features assessed (Supplementary Table S1). Fine scale BPI used an inner radius of 1 cell (100 m) and an outer radius of 10 cells (1000 m). Broad-scale BPI used an inner radius of 25 cells (2500 m) and an outer radius of 50 cells (5000 m).

Substrate type for each deployment was classified into either hard or soft bottom habitat. Hard bottom habitats were generally rocky outcrops or pavement and consisted of <20% sand, whereas soft-bottom habitats were areas with greater than ≥80% sand (Table 1).

TABLE 1

Table 1. Number of drops per site, including number per benthic habitat type (hard vs. soft), and number per depth zone (Mesobenthic vs. Bathybenthic).

Statistical Methods

Biodiversity Indices

Hill numbers (or the effective number of species) have been increasingly used to quantify the species/taxonomic diversity of an assemblage. The sample-size- and coverage-based integrations of rarefaction (interpolation) and extrapolation (prediction) of Hill numbers represent a unified standardization method for quantifying and comparing species diversity across multiple assemblages. To assess biodiversity, taxon observations were converted to incidence matrices. Three Hill numbers were calculated for each site: species richness, Shannon diversity, and Simpson diversity (Hill, 1973; Chao et al., 2014). We extrapolated biodiversity and assessed sample completeness for each location in the R package iNEXT Version 2.0.19 (Hsieh et al., 2016). Differences in biodiversity among locations, habitats (hard vs. soft), and depth zones (Table 1) were tested using linear analysis of variance (ANOVA) with a Tukey HSD post hoc test in the R programing environment, Version 1.0.136 (R Core Team, 2016).

Generalized Linear Mixed Effects Models

We used Generalized Linear Mixed Effects Models (GLMMs) to investigate environmental drivers of biodiversity among locations, as well as faunal assemblage composition among sites (islands/archipelagos) and sea-floor types (hard vs. soft bottom habitat).

To test for environmental drivers of biodiversity, the response variables (richness, Shannon diversity, Simpson diversity) were modeled as a Poisson distribution. For all hierarchical mixed models, location (island/archipelago) was included as a random effect to account for the nested sampling structure of the dataset. Model parameter significance was assessed using a Likelihood Ratio Test (LRT) with the statistical package LmerTest (Bates et al., 2014). All data analyses were performed in the R programing environment, Version 1.0.136 (R Core Team, 2016).

For community composition models, biological data were subset for taxa that occurred in more than three samples (drops). The probability of taxa occurrence was assessed using a GLMM with random slopes and intercepts for each taxon, allowing the taxon-specific relationship to vary by location or sea-floor type (hard vs. soft bottom), and modeled as a binomial distribution in the statistical package lme4 (Bates et al., 2014). The mixed model allowed individual taxon to respond differentially to location or sea-floor type, and the change in community composition was quantified in terms of how much the slopes of the response differ across taxa.

Multivariate Analysis

To describe the patterns of variation in community structure (patterns of distribution of abundance of fish families and invertebrate taxa) and their relationship with environmental data (Supplementary Table S1) we used indirect gradient analysis. Non-linear models were most appropriate for our data because a preliminary detrended correspondence analysis showed long gradient lengths (>2 SD) (ter Braak and Šmilauer, 2002). We performed Canonical Correspondence Analysis (CCA) on square root-transformed numerical abundance of fish families and invertebrate taxa using the ordination program CANOCO for Windows version 5.0 (ter Braak, 1994). Fish families and invertebrate taxa that occurred in <5% of the samples were excluded from the analyses. Environmental data were centered and standardized for analyses.

Results

Fish and Invertebrate Assemblages of the Bathyal Zone

From a total of 70 deployments of the Deep-Ocean Dropcams, 43 different fish families and 42 different invertebrate orders were identified across four survey locations in the TEP Ocean Island Province. Hard and soft habitats were surveyed roughly evenly among the four locations (Table 1). Within the bathyal zone (200 m–3,000 m), mesobenthic (200 m–1,500 m), and bathybenthic (<1,500 m–3,000 m) were also sampled roughly evenly (Table 1).

For fishes, comparisons between habitat types (hard vs. soft bottom) showed no difference in richness (F_1,78 = 0.11, P = 0.74), Shannon (F_1,68 = 0.01, P = 0.93), or Simpson (F_1,68 = 0.01, P = 0.94) diversity, therefore samples were pooled among habitats for further analyses. For invertebrates, soft bottom habitats had significantly lower diversity than hard bottom habitats for all three diversity indices: (F_5,64 = 9.66, P = 0.001), Shannon (F_5,64 = 12.42, P = 0.002), and Simpson (F_5,64 = 11.86, P = 0.024).

Biodiversity Among Locations

Fishes

For all sites, three diversity indices were in agreement, showing consistent trends across sites (Figure 2). Asymptotic diversities indicated that the Galápagos were highest in richness [index value 25 (SE 0.27)], followed by Malpelo [index value 17 (SE 0.72)]. The Revillagigedo [index value 10 (SE 1.69)], and Clipperton [index value 11 (SE 3.09)] had similarly, low diversities and overlapping 95% confidence intervals on the rarefaction/extrapolation curve (Figure 2). Assessment of sample completeness for locations were: Clipperton (84.3%), Galápagos (99.1%), Malpelo (100%), and Revillagigedo (94.9%).

FIGURE 2

Figure 2. Rarefaction/extrapolation curves for asymptotic fish family diversity indices: Richness, Shannon, and Simpson calculated by incidence frequency matrices. Sample-size-based rarefaction (solid line segment) and extrapolation (dotted line segments) sampling curves with 95% confidence intervals (shaded areas) for four locations in the TEP, separately by diversity family. The solid dots/triangles represent the reference samples.

Richness (F_3,76 = 6.91, P < 0.001), Shannon (F_3,66 = 8.81, P < 0.001), and Simpson (F_3,66 = 10.68, P < 0.001) diversity were statistically different among locations based on linear analysis of variance. Clipperton had significantly lower fish family richness compared to all other locations (P < 0.05). Shannon and Simpson diversity were significantly lower at Clipperton compared to Galápagos and Malpelo (P < 0.001). Although not significant, differences in Simpson diversity between Clipperton and Revillagigedo were suggestive (P = 0.06). Galápagos, Malpelo, and Revillagigedo overlapped in 95% confidence intervals around mean richness and diversity estimates. Geographically, locations closer to the equator and closer to the continent had higher richness and diversity, while assemblages in locations farther north and farther west were lower in species richness and diversity.

Invertebrates

Similar to fish biodiversity methods, rarefaction/extrapolation curves were constructed for invertebrate diversity indices (Figure 3). For all of these sites, the three diversity indices were in agreement, showing consistent trends across sites. Asymptotic diversities indicated that the Revillagigedo [index value 13.5 (SE 10.9)] had lower richness than the other three sites, with Clipperton having 24.3 (SE 12.2), Galápagos 25.6 (SE 9.8), and Malpelo 43.0 (SE 23.3) for richness asymptotic estimates. Assessment of sample completeness for locations were: Clipperton (84.6%), Galápagos (95.1%), Malpelo (89.1%), and Revillagigedo (87.3%).

FIGURE 3

Figure 3. Rarefaction/extrapolation curves for asymptotic invertebrate order diversity indices: Richness, Shannon, and Simpson calculated by incidence frequency matrix. Sample-size-based rarefaction (solid line segment) and extrapolation (dotted line segments) sampling curves with 95% confidence intervals (shaded areas) for four locations in the TEP, separately by diversity order. The solid dots/triangles represent the reference samples.

For invertebrates, richness (F_3,66 = 10.68, P < 0.001), Shannon (F_3,66 = 15.6, P < 0.001), and Simpson (F_5,64 = 11.86, P = 0.02) diversities were statistically different among locations based on linear analysis of variance. The Revillagigedos were significantly lower compared to the Galápagos and Malpelo for richness (F_3,66 = 10.68, P < 0.001), Shannon (F_3,66 = 10.68, P < 0.001), and Simpson diversity (F_5,64 = 11.86, P < 0.001).

Environmental Drivers of Biodiversity

Fishes

Twelve environmental variables were included in the global GLMM model (Supplementary Table S1). The reduced model, selected by AICc, included broad scale BPI and dissolved oxygen. In all fish models, both BPI and dissolved oxygen had a positive effect on fish diversity metrics. The pattern was stronger with Richness and Shannon diversity, showing significant positive effects. Richness (β = 0.151, se = 0.08, P = 0.05) and Shannon diversity (β = 0.12, se = 0.06, P = 0.05) were both marginally significant [richness: LRT: χ²(2) = 3.56, P = 0.05; Shannon: LRT: χ²(2) = 3.78, P = 0.05], but only explained a small portion of the variability in these data (Richness R² = 0.17; Shannon R² = 0.30) (Figures 4A,B, respectively). Simpson diversity had a positive effect of BPI, though not statistically significant (β = 0.04, se = 0.03, P = 0.13) [LRT: χ²(2) = 2.31, P = 0.12]. Dissolved oxygen positively effected diversity metrics, though not significant for all three metrics tested; (richness: β = 0.16, se = 0.09, P = 0.08; Shannon: β = 0.07, se = 0.06, P = 0.25; and Simpson: β = 0.03, se = 0.08, P = 0.27).

FIGURE 4

Figure 4. GLMM fixed effects from fish family Richness (A) and Shannon diversity (B), for scaled and centered dissolved oxygen (left panels) and Bathymetric Position Index (BPI) (right panels). Blue bands are 95% confidence intervals.

Invertebrates

Invertebrate biodiversity models showed no significant effect of BPI and dissolved oxygen concentration. Hierarchical models with location as a random effect significantly explained variation in the invertebrate taxa data for Shanon and Simpson diversity, but not for richness [Shannon diversity LRT: χ²(5) = 6.17, P = 0.01; Simpson diversity LRT: χ²(5) = 8.69, P < 0.01; richness LRT: χ²(4) = 1.21, P = 0.27].

Assemblage Composition Among Locations

Fishes

Models showed significantly different assemblages [LRT: χ²(14) = 56.2, P < 0.001] among locations. The change in community composition was quantified in terms of how much the slopes of the response differ across families. The magnitude of variation was 5.35 (se = 0.71) on a logit scale, indicating that the change in community composition was substantial in terms of which fish families were present, depending on location. The variation in the location effect across families was large as standard deviation ranged from 0.5 to 1.8 (on the logit scale) depending on location (Figure 5).

FIGURE 5

Figure 5. GLMM random effects (modeled slopes) for fish family specific probability of occurrence by location.

Many fish family occurrences were highly variable. Their modeled slopes overlapped with zero in all but four - six families per location (Figure 5). Grenadiers (Macrouridae) were common at three locations (with less at the Revillagigedos), with the highest association at Clipperton. Chimaeras (Chimaeridae) were common at Clipperton and the Revillagigedos, but were not frequently encountered at the Galápagos or Malpelo. Catsharks (Scyliorhinidae) were the most common family encountered at the Galápagos, but less common at Malpelo. Morids (Moridae) were less common at Galápagos and Malpelo (Figure 5).

Similar binomial GLMMs were used to test for family specific effects of habitat (hard vs. soft). While these did show a slight shift in community composition these were not statistically significant [LRT: χ²(5) = 4.90, P = 0.08] (Figure 6).

FIGURE 6

Figure 6. GLMM random effects (modeled slopes) for fish family specific probability of occurrence by habitat (hard vs. soft bottom).

Fish family assemblage structure could be explained by several environmental variables. The first two axes of the CCA plot explained 27% of the explained variation and 66% of the explained fitted variation (Figure 7 and Table 2). The main factors influencing this ordination were location (e.g., Malpelo, Galápagos, and Revillagigedo), distance to seamounts, and distance to the continent. Chimaeridae and Moridae were correlated with depth and bathybenthic strata. Liparidae, Echinorhinidae, and Etmopteridae were most closely correlated with seamounts. Rajiidae and Scyliorhinidae were correlated with soft bottoms and closer to the continent.

FIGURE 7

Figure 7. Canonical Correspondence Analysis (CCA) of fish family numerical abundance. Vectors indicate direction and strength (length) of relationships between taxa and explanatory geospatial variables defined in Supplementary Table S1 [e.g., location; distance to seamounts; continent; or slope; depth; bathybenthic strata (bathybenthic vs. mesobenthic); and habitat type (hard vs. soft); particulate organic carbon (POC); dissolved oxygen (DO); and temperature].

TABLE 2

Table 2. Results of Canonical Correspondence Analysis (CCA) on square root-transformed fish abundance data.

Invertebrates

Models showed significantly different assemblages [LRT: χ²(9) = 69.9, P < 0.001] among locations. The change in community composition was quantified in terms of how much the slopes of the response differ across orders. The magnitude of variation was 5.35 (se = 2.35) on a logit scale, indicating that the change in community composition was substantial in terms of which invertebrate orders were present, depending on location. The variation in the location effect across orders was large. The standard deviation ranged from 1.3 to 3.6 (on the logit scale) depending on location (Figure 8).

FIGURE 8

Figure 8. GLMM random effects (modeled slopes) for invertebrate order-specific probability of occurrence by location.

Each site was unique in terms of invertebrate order probability of occurrence (Figure 8). The Revillagigedos showed the least variability in estimates, but no strong affinities, whereas Malpelo showed a large variability, with sea stars (Paxillosida) strongly negatively associated with the site, and amphipods (Amphipoda), and krill (Euphausiacea) strongly positively associated with the site.

Similar binomial GLMMs were used to test for invertebrate order-specific effects of habitat (hard vs. soft). These did not show a shift in community composition [LRT: χ²(2) = 1.36, P = 0.50] (Figure 9). The signal and variation were greater in hard bottom habitats (standard deviation 0.79 on the logit scale) compared to soft bottom habitat (standard deviation 0.24 on the logit scale).

FIGURE 9

Figure 9. GLMM random effects (modeled slopes) for invertebrate order-specific probability of occurrence by habitat (hard vs. soft bottom).

Invertebrate community structure was correlated with a number of environmental variables. The first two axes of the CCA plot explained 26% of the explained variation and 65% of the explained fitted variation (Figure 10 and Table 3). The main factors influencing this ordination were location (e.g., Clipperton, Malpelo, Galápagos, and Revillagigedo), distance to seamounts and distance to the continent. POC and slope also explained modest amounts of variation in invertebrate community structure. Echinodermata, Asteroidea, and Sceptrulophora were correlated with seamounts. Alcyonacea was most closed correlated with proximity to the continent, while Ophiurida and Amphipoda were correlated with POC.

FIGURE 10

Figure 10. Canonical Correspondence Analysis (CCA) of invertebrate taxa numerical abundance. Vectors indicate direction and strength (length) of relationships between taxa and explanatory geospatial variables defined in Supplementary Table S1 (e.g., location; distance to seamounts; and continent; and POC).

TABLE 3

Table 3. Results of Canonical Correspondence Analysis (CCA) on square root-transformed invertebrate abundance data.

Discussion

Patterns in Biodiversity Across Islands/Archipelagos

The overall geographic pattern of biodiversity revealed higher values closer to the equator and closer to the continent, and lower values in the locations farther north and farther west. This finding is consistent with the biodiversity pattern predicted based on marine species ranges across latitudinal gradients (Costello and Chaudhary, 2017). Fish family biodiversity was highest in the Galápagos, followed by Malpelo. In contrast, Clipperton had significantly lower fish family richness compared to all other locations. For invertebrate orders, the Revillagigedos were significantly lower in biodiversity compared to the Galápagos and Malpelo.

The strong location effect could be explained, in part, by the vastly different current regimes around each island group, and the unique hydrographic and environmental characteristics of the associated water masses (Figure 1). Gradients in environmental parameters, such as temperature, salinity, nutrients, and oxygen, vary drastically across this region and co-occur within larger oceanographic processes, such as current regimes [see review by Fiedler and Talley (2006)]. These interrelated variables affect patterns in biodiversity through influencing the physiology, biology, and interactions of marine organisms (Levin et al., 2001). Similarly, co-varying environmental parameters have been implicated in causing depth zonation in the deep-sea (Rex, 1976; Carney, 2005; Yeh and Drazen, 2009) and latitudinal gradients in biodiversity (Fortes and Absalão, 2010; Costello and Chaudhary, 2017). While depth zonation has been well documented, current regimes will require further assessment as large-scale drivers of biodiversity within the region (Levin et al., 2001).

Benthic habitat types (hard vs. soft bottom) had little influence on fish diversity in our study. This finding is not surprising given that highly mobile bathyal demersal fishes are expected to have large ranges, irrespective of habitat type (Ekman, 1953; Pineda, 1993; Stevens, 1996; Smith and Brown, 2002). Of the 284 described fish species with circumglobal distributions, more than 80% of these species are represented by deep-sea fishes (Gaither et al., 2016). In contrast, for the low-mobility macro invertebrates assessed in this study, soft bottom habitats harbored significantly lower diversity compared to hard bottom habitats for all three diversity indices. This finding is consistent with the well-documented positive association of low-mobility benthic fauna diversity and measures of habitat heterogeneity (Etter and Grassle, 1992; Leduc et al., 2012).

Faunal Assemblage Composition Across Islands/Archipelagos

Our results show significant differences in community composition among locations for fish and invertebrate assemblages. Variations in the spatial distribution in fish families were mostly explained by location, distance to seamounts, distance to the continent, and depth. Temperature, and broad-scale complexity (BPI) were also important in explaining the distribution of fish family abundance within the region. Similarly, invertebrate community structure was explained by location, distances to seamounts and to the continent, as well as POC and bathymetric slope.

Habitat effects (hard vs. soft) on fish family composition were suggestive, though not statistically significant in our GLMMs, a finding that is similar to comparisons of abyssal scavenging community composition in the Pacific (Leitner et al., 2017). In contrast, studies of bathyal fish-habitat associations on seamounts in the Atlantic found significant differences in community composition between substrate types (Porteiro et al., 2013). The loose habitat association detected here could be expected given that the assemblage of scavengers/predators observed are highly mobile. In contrast, for invertebrates, the shifts in composition with habitat were significant. The order-specific signal and variation in habitat effects were greater in hard compared to soft bottom habitats, likely due to increased habitat heterogeneity (Etter and Grassle, 1992). The positive associations could be related to a greater abundance and diversity of potential epibenthic prey items for less-mobile invertebrate taxa (McClain and Barry, 2010).

Overall, measures of seascape configuration (e.g., distance to seamounts, distance to the continent) were stronger explanatory variables structuring community composition in our study, compared to gradients in environmental parameters (e.g., temperature, salinity) (Supplementary Table S1). These co-varying environmental parameters respond to large-scale oceanographic processes, such as currents, and are affected by seascape configuration (Levin et al., 2001; McClain and Barry, 2010).

Environmental Predictors of Deep-Sea Biodiversity

We measured a significant relationship between bathyal fish biodiversity and bathymetric predictors, namely broad-scale (2,500 m) BPI. As a topographic variable, BPI is used to identify benthic features (e.g., ridges and valleys), with positive values indicating peaks in topography while negative values indicate troughs. Fish diversity increased with BPI, indicating a positive association with ridges and peaks as opposed to valleys on the seafloor. Our findings are consistent with previous studies in the Pacific (Leitner et al., 2017), which also detected changes in diversity along topographic gradients. Broad-scale BPI provides a metric for bathymetric characteristics at a spatial scale appropriate for highly mobile species, such as those that dominate the fish communities assessed in this study. Lower-mobility invertebrate communities did not show a significant relationship with either fine- or broad-scale BPI. Likely these communities are influenced by topographic metrics at a scale finer than the BPI scale (minimum inner radius 100 m) afforded by the resolution of our bathymetric data (Dolan et al., 2008; Guinan et al., 2009; Turner et al., 2018).

Dissolved oxygen concentration emerged as a positive significant predictor of fish diversity. Annual mean dissolved oxygen values derived from depth bins associated with deployment sites indicate that many of the biological communities sampled would have been under severely hypoxic conditions (<0.2 ml l⁻¹) (Figure 1). More specifically, deployments at intermediate depths were located in oxygen minima. In the TEP, oxygen minima lie below the surface layer where high productivity peaks oxygen consumption, and above depths replenished by advection of oxygenated waters subducted at higher latitudes (Fiedler and Talley, 2006). These oxygen minimum zones (OMZs) typically exhibit low macrofaunal and megafaunal species richness, with communities highly dominated by only a few taxa (Levin et al., 2001). This pattern was apparent for fish communities in our study, and has been observed elsewhere in the eastern Pacific (Levin et al., 1991). This pattern of increased community diversity with dissolved oxygen concentration has also been observed for invertebrate macrofauna (Diaz and Rosenberg, 1995) and megafauna (Levin, 2003), though was not detected in our study. As large taxa are more sensitive than small taxa to hypoxia (Levin et al., 2001), we expect that with more extensive sampling, and finer resolution taxonomic classification, patterns of association between dissolved oxygen concentration and invertebrate diversity may be revealed. These may not be linear, but rather unimodal, with peaks in diversity at intermediate levels of oxygen concentration (Levin et al., 2001).

We did not see positive influence of POC flux (a proxy for food availability) on biodiversity, a surprise given that primary productivity is implicated in generating latitudinal diversity gradients (Rex et al., 1993). However, the relationship between biodiversity, horizontal variation in vertical POC flux, and the many interrelated factors across water masses and current regimes that co-occur with latitude, requires additional sampling (Lampitt and Antia, 1997; Lutz et al., 2007; McClain and Hardy, 2010). For example, food availability to the seafloor declines exponentially with depth (Lutz et al., 2007; Smith et al., 2008) therefore POC effects could be dampened by the strong influence of bathymetry on surrounding environments. In abyssal depths of the eastern Pacific, peaks and ridges, rather than the overlying surface productivity were implicated in influencing abyssal scavenging community composition and diversity at local scales, owing to either decreased food availability in bathymetric lows, or to enhanced flow and food flux to suspension feeders in bathymetric highs (Leitner et al., 2017).

Limitations and Data Gaps

Many environmental predictor datasets in our study remain at a coarse scale, potentially shrouding relationships between predictors and biodiversity. For example, the oxygen data applied in our analysis were for coarse vertical zones (i.e., 800–3500 m; Watling et al., 2013), therefore the insignificant relationship with invertebrate orders could be an artifact of data resolution. More extensive sampling, and ideally in situ environmental measurements, would greatly enhance our ability to determine the individual effects of environmental parameters on patterns in biodiversity. Nevertheless, this study represents the first assessment of bathyal demersal biodiversity across multiple Ocean Island Province sites in the TEP, and while sampling coverage limits our ability to detect the individual influence of co-occurring hydrographic and environmental variables, significant differences among locations in fish and invertebrate diversity were apparent. These may be attributed, within the limitations of the study, to the different current regimes around each island/archipelago.

Conclusion

Biodiversity and community composition of bathyal demersal fauna differ across the TEP associated with large-scale oceanographic processes. We have provided insight into differences in biological communities among four island groups, and identified several environmental drivers of these variations, including dissolved oxygen concentration and habitat type for some fauna. This expands the previously limited knowledge of biodiversity patterns and its drivers in the TEP. Because the TEP is largely isolated from the wider Pacific (Lessios and Robertson, 2006), the islands and archipelagos of the Ocean Island Province likely represent important stepping-stones and corridors for species connectivity across the TEP and the wider Pacific (Glynn et al., 1996). Therefore, the preservation of these intact ecosystems is important for maintaining biodiversity both within the TEP as well as across the wider Pacific region.

Data Availability

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

JG, WG, and AF wrote the manuscript, performed the statistical analyses, processed the imagery, and identified the species. PS identified the species. BH, CS, and AT built the Dropcams and collected the data on expeditions. EB, AT, and ES organized the research program. All authors assisted in manuscript editing.

Funding

Funding for this study was provided by the National Geographic Society, Washington, DC, United States.

Conflict of Interest Statement

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.

Acknowledgments

We thank the Pristine Seas team and boat crews who assisted in data collection. Jessica Elfadl of the Exploration Technology Lab assisted tremendously with expedition logistics and equipment. We thank the experts who assisted with species identification, including Chris Kelley, Ross Robertson, Bruce Mundy, Beatriz Mejía-Mercado, and Dave Ebert. Graham Wilhelm was instrumental in setting up the initial mechanical design of the Dropcams. We thank Graham and all members of the National Geographic Exploration Technology Lab. The Galápagos Islands research was authorized by the Galápagos National Park Directorate under the permit PC-26-15. This publication is contribution number 2275 of the Charles Darwin Foundation for the Galápagos Islands. The insightful comments of the two reviewers and the editor greatly improved the quality of this manuscript.

Supplementary Material

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

References

Allen, G. R. (2008). Conservation hotspots of biodiversity and endemism for Indo-Pacific coral reef fishes. Aquat. Conserv. 18, 541–556. doi: 10.1002/aqc.880