Environmental Drivers of Mesophotic Echinoderm Assemblages of the Southeastern Pacific Ocean

Mesophotic ecosystems (50–400 m depth) of the southeastern Pacific have rarely been studied because of the logistical challenges in sampling across this remote zone. This study assessed how oxygen concentrations and other environmental predictors explain variation in echinoderm assemblages at these mesophotic systems, where this group is among the predominant fauna. We compiled data on echinoderm taxa at 91 sampling stations, from historical and recent surveys (between 1950 and 2019), covering a longitudinal gradient of approximately 3,700 km along with the Nazca, Salas y Gómez, and Juan Fernández ridges. Uni- and multivariate model-based tools were applied to analyze the patterns of benthic fauna in relation to environmental factors. Our results indicate a significant positive relationship between echinoderm species richness and depth, oxygen, and salinity. Changes in echinoderm community composition were significantly explained by oxygen, longitude, and chlorophyll-a. We observed notable species turnovers at ∼101 and ∼86°W, where assemblages tend to be more variable across stations. This turnover possibly reflects the effects of physical barriers to dispersion (e.g., currents) and habitat changes. Echinoderm assemblages observed around Easter and Desventuradas Islands presented a high occurrence of potentially endemic taxa and distinct species assemblages. This study is the first to assess the structure of mesophotic echinoderm assemblages of the southeastern Pacific Ocean along a large spatial scale. The information reported here could help design appropriate management tools for the vast, recently created, marine protected areas in the southeastern Pacific.


INTRODUCTION
The study of vulnerable marine ecosystems, such as seamounts and oceanic islands, is critical for the conservation and management of the marine ecosystem (Kvile et al., 2014;Watling and Auster, 2017). Seamounts are topographic structures that rise more than 1,000 m above the surrounding seafloor. Their topography affects marine currents, producing local and mesoscale circulation patterns, such as Taylor columns that usually reinforce vertical mixing (Vic et al., 2019), increase nutrient concentrations in the euphotic zone , and promote the growth of primary producers that sustain benthic and pelagic animals (Morato et al., 2010). This increased productivity and the availability of hard substrate on seamounts and oceanic islands allow the establishment of diverse marine communities (Morato et al., 2013). Such communities frequently have high endemism and potentially low resilience levels because their populations are usually restricted to relatively few individuals with limited distributions compared with their continental counterparts (Clark and Rowden, 2009;Clark et al., 2010). According to the International Guidelines for the Management of Deep-sea Fisheries in the High Seas (FAO, 2009), "a marine ecosystem should be classified as vulnerable based on the characteristics that it possesses" and those ecosystem characteristics could be, for example, their uniqueness, rarity, structural complexity or fragility (Paragraph 42, FAO DSF Guidelines). Hence, many seamounts and oceanic islands could be classified as vulnerable marine ecosystems because of their potentially fragile communities and habitat-forming species they support (Watling and Auster, 2017). Besides, the scattered distribution of islands and seamounts along oceanic ridges makes them important "stepping stones" for the distribution of species across ocean basins and, therefore, contributes to shaping biogeographical patterns of oceanic biodiversity (Miller and Gunasekera, 2017).
Echinoderms usually constitute an important group of motile invertebrate grazers, carnivores, and deposit feeders in benthic habitats (Furman and Heck, 2009) and provide essential ecosystem services. Most asteroids are carnivores, and some species have an important role in benthic ecosystems since they can control populations of dominant fouling organisms (Paine, 1966). The impact of these asteroids on the benthic community structure favors the local biodiversity, making them true keystone species (Paine, 1966). Some echinoderms are important ecosystem engineers in many marine systems since they provide shelter and refuge for juvenile fish recruits (Fonseca et al., 2014) and affect substrate characteristics through sediment bioturbation (Uthicke, 2001;Shiell and Knott, 2010) and accumulation of carbonate spines and skeletons after their death (Ogden, 1977). Their ecological importance and high abundance in benthic ecosystems make them good indicators of the marine ecosystem health (Gonzalez-Irusta et al., 2012;Mironov et al., 2014); however, little information on their diversity and ecological role in mesophotic ecosystems is available Pyle and Copus, 2019). This lack of information arises, in part, because the communities inhabiting seamounts and oceanic islands at mesophotic depths are largely understudied compared with shallower and deeper communities. The mesophotic zone is characterized by low-light conditions that, in general, extend from depths of ∼30 m to the limits of photosynthetically dependent biota at ∼150 m (Pyle and Copus, 2019) or even deeper at some ultraoligotrophic places, e.g., to ∼300 m (Hinderstein et al., 2010;Easton et al., 2019). This lack of information is particularly evident for the seamounts and oceanic islands scattered within the southeastern Pacific Ocean (SEP). Although the benthic systems of a few oceanic islands in this region have been studied since the late 1950s, most surveys were conducted at depths between 0 and 50 m (Ziesenhenne, 1963;Fell, 1975) or below 150 m (Parin et al., 1997). Although these studies have expanded our knowledge of the biodiversity of echinoderms in the SEP, few attempts were made to explain benthic faunal patterns in relation to abiotic factors. Parin et al. (1997) described the presence of low oxygen conditions (bottom water dissolved oxygen: 0.15-0.20 mL/L) on shallow seamounts close to the junction of the Salas y Gómez and Nazca ridges, but they did not explore correlations between echinoderm communities and dissolved oxygen or other abiotic factors (e.g., depth and temperature) that can be major drivers of echinoderm community structure. For example, O'Hara and Tittensor (2010) observed that temperature was a stronger driver than depth for ophiuroids on southwest Pacific seamounts; and Iken et al. (2010) found that salinity, sea surface temperature, chlorophyll-a, and primary productivity were strongly correlated with echinoderm assemblages.
Since the easternmost islands and seamounts of the SEP extend to areas under the influence of oxygen minimum zones, especially near the South American continent (Silva and Neshyba, 1979), it is expected that oxygen concentration would be a key factor shaping echinoderm diversity patterns in mesophotic benthic ecosystems in the SEP. In this way, to assess the effects of oxygen and other potential drivers of echinoderm community structure in the SEP, we analyzed the variations in echinoderm species richness and assemblages between 50 and 400 m on seamounts and oceanic islands of the SEP and analyzed the relationship between these patterns and those of select environmental factors (i.e., latitude, longitude, depth, temperature, oxygen, salinity, and chlorophyll-a). To do so, we compiled a dataset that incorporates historical biological data from publications and online databases and data collected during recent oceanographic cruises (2016-2019).

Study Area
We compiled data from the Desventuradas Islands, Easter Island, and seamounts from three ridges within the SEP: Salas y Gómez, Nazca, and Juan Fernández (Figure 1). These three ridges contain more than 144 seamounts (51 within the Chilean EEZ) with summits from 0 to 2,000 m depth representing 54% of the reported seamounts in the SEP (Gálvez-Larach, 2009). The Salas y Gómez Ridge was formed by the Easter-Salas hotspot and extends from ∼112 to ∼80 • W, where it converges with the Nazca Ridge (Ray et al., 2012). Although it is relatively narrow (∼150 km of latitudinal extension), the Salas y Gómez Ridge spans more than 2,900 km along its longitudinal axis. The islands of Rapa Nui (Easter Island) and Motu Motiro Hiva (Salas y Gómez) are the only emerged structures along this ridge. No emerged structures are present along the Nazca Ridge, which also formed from the Easter-Salas hotspot. The Nazca Ridge extends latitudinally from ∼15 to ∼25 • S (∼1,300 km) and longitudinally from ∼76 to ∼85 • W (∼1,200 km). The Desventuradas Islands, which are small islands of volcanic origin, are located east of the junction of the Salas y Gómez and Nazca ridges (Gálvez-Larach, 2009). The Juan Fernández Ridge lies on a ∼800 km (longitudinally) ridge named after the Juan Fernández Archipelago and extends latitudinally at ∼33 • S and longitudinally from ∼76 to ∼81 • W. The Archipelago is composed by three small islands (Robinson Crusoe and Alejandro Selkirk Islands and Santa Clara Islet) located ∼840 km south of the Desventuradas Islands and ∼670 km west of the continental coast of Chile.
We assigned the islands and seamounts according to geographic latitude and longitude to one of six regions (Figure 1): (i) Easter Island, (ii) Salas y Gómez Ridge, (iii) Salas y Gómez Ridge/Nazca Ridge junction, (iv) North Nazca Ridge, (v) Desventuradas Islands, and (vi) Juan Fernández Ridge.

BIOLOGICAL DATASET
The extreme oligotrophy of the waters in the westernmost areas explains the presence of light-dependent biota (i.e., crustose coralline algae) to depths of ∼280 m at Easter Island  and depths > 300 m (potentially to almost 400 m, unpublished data) on a shallow seamount of the Salas y Gómez Ridge. To include these deep mesophotic communities, we expanded the target depth range of this study to 400 m.
Echinoderm data from 91 stations (50-400 m depth) sampled between 1957 and 2019 (data provided in the Supplementary Material) was compiled, of which 36 stations were sampled during recent expeditions: CIMAR22 (Crucero de Investigación Científico-Marina en Áreas Remotas 22), EPIC (East/Central Pacific International Campaign) and Easter Island expeditions (2014)(2015)(2016). During the CIMAR 22 cruise, which is part of a multidisciplinary research project funded by the Chilean government, 25 stations were sampled by Agassiz trawl or ROV within the Chilean EEZ around the archipelagos of Desventuradas and Juan Fernández (Figure 1). The cruise was carried out on the oceanographic vessel AGS 61 Cabo de Hornos (Chilean Navy) during October and November 2016. During the multidisciplinary oceanographic EPIC cruise onboard the RV Mirai (JAMSTEC, Japan), four stations along the Salas y Gómez Ridge were sampled with Agassiz trawl in January and February of 2019. Data from seven stations sampled between 2014 and 2016 at 160-280 m around Easter Island were collected by ROV . These stations include a subsurface peak, located 13 km SW of the island, locally known as Apolo, and a seamount known as Pukao, located 98 km west of the island. Data from the remaining 58 stations were compiled from existing publications (Allison et al., 1967;Parin et al., 1997) and online databases such as the Ocean Biogeographic Information System (OBIS 1 ) and Seamounts Online 2 .

Environmental Data
Physical (i.e., temperature, salinity) and chemical (i.e., oxygen) data were obtained from the global three-dimensional CSIRO Atlas of Regional Seas (CARS) climatology (2009 version) and EPIC cruise. CARS climatology combines all the available oceanographic data from the World Ocean Database over the last 50 years, along with the Array for Real-time Geostrophic Oceanography (ARGO) buoy profiles. CARS is based on an adaptive-length-scale loess mapper to maximize resolution in data-rich regions and the presence of steep topography (Ridgway et al., 2002). This procedure results in a better definition of oceanic structures (e.g., coastal upwelling off Chile) and accuracy of point values than other global ocean climatologies (e.g., NODC's World Ocean Atlas). The horizontal resolution is 0.5 • . As data availability has significantly increased in recent years, the CARS mean values are inevitably biased toward the current ocean state. Data were selected from the closest location (longitude, latitude) to the 91 stations of the cruise data and vertically interpolated on a 1-meter vertical grid using a cubic spline function to match their depth. The CARS data, therefore, provide 12 monthly mean values for each variable corresponding to the calendar months. The results presented in the paper used the annually averaged values; however, we also use monthly mean calendar data (12 mo) to test the robustness of some of our analyses (cf. Supplementary Figure S1B) by performing them using the minimum and maximum values as environmental predictors. Preliminary analysis showed that minimum and maximum CARS data values were highly colinear with the annual averages (Pearson correlation coefficient > 0.99) and, therefore, were not included as predictors in most analyses. The sensitivity to monthly variation in CARS predictors within and across sampling stations was assessed by performing a principal component analysis (PCA, Supplementary Figure S1). A convex hull represented intra-annual variation within each station and the mean value for each station as a point (Supplementary Figure  S1A). The oceanographic data (oxygen, depth, salinity, and sea surface temperature) collected during the EPIC cruise were also included in this PCA. The data was scaled and centered prior to the PCA, which was performed using the R package vegan (Oksanen et al., 2018).
Data on the superficial concentrations of chlorophyll-a (mg m −3 ) was also considered in our analyses. This data was obtained from the Bio-Oracle online database (Tyberghein et al., 2012;Assis et al., 2018), which provides data on different oceanographic variables for applications in studies on marine biogeography. These global maps represent an average of 2000-2014 and were generated from local observations, satellite images, and modeled estimates (Tyberghein et al., 2012;Assis et al., 2018). The data on superficial chlorophyll-a concentrations were extracted from a global raster with a spatial resolution of 5 arcmins (c. 9.2 km) and a uniform landmask. The chlorophyll-a values at each sampling station were extracted considering a spatial buffer of 30 km (i.e., the average of all values within a 30 km radius around the sampling station). This values extraction was made by using tools in the R package raster (Hijmans, 2019).

Data Analysis
All analyses and graphs were performed using the software R and associated packages (R Core Team, 2019) and maps made with QGIS 3.6.2 Version.
We used taxa richness (S) as a proxy of echinoderm diversity (Tecchio et al., 2011). Given the differences in sampling effort and techniques across stations, we applied a data standardization protocol before analyzing taxa richness. To apply a rarefactionbased standardization, we grouped multiple stations within sites when they were distributed on and around the same seamount or oceanic island. Thus, sampling stations were considered as subsamples of sites in a particular area. Sparsely distributed stations were included within the nearest site when the distance to it was smaller than 50 km. Stations of North Nazca were included in the same site since these three isolated stations were ∼50 km apart from each other. Single stations that were more than 50 km apart from any site or station were not considered in the taxa richness analysis. We used tools of the R package iNEXT (Hsieh et al., 2016) to standardize the estimations of taxa richness for each site via seamless rarefaction and extrapolation sampling curves. We also used the iNEXT package to assess sample completeness (sample coverage) across multiple sites. We applied asymptotic diversity estimates for each site to extrapolate taxa richness. Estimates of sites with low sample completeness (<0.3) were not included in the further analyses. This process resulted in the estimates of taxa richness of a total of fifteen sampling sites.
The relationship between taxa richness (number of species) and the environmental variables was assessed by applying generalized linear models (GLMs). We used routines for model selection and averaging in the R package MuMIn (Barton, 2009) to generate a final model containing the best predictor estimates to model echinoderm taxa richness. First, we constructed a global GLM based on Poisson distribution and log-link function that considered latitude, longitude, depth, salinity, dissolved oxygen, superficial chlorophyll-a, sea surface temperature, and water temperature at 300 m depth as predictors. We applied a model selection routine to classify different subsets of the global model based on their performance (based on the corrected Akaike Information Criterion, AICc). The best model fits (i.e., AICc weight > 0.5) were then used to calculate the full average of the model estimates. This process generated a final averaged model with the predictor estimates that better explain taxa richness variation. The general effect of each predictor was tested by using an analysis of deviance and assessed by plotting the partial residuals against each predictor.
We visualized the general trends in the variation of echinoderm assemblages across stations by fitting pure latent variable models. Pure latent variable models offer a model-based approach to unconstrained ordination for visualizing sites and the indicator species characterizing them on a low-dimensional plot (Hui, 2016). The ordinations produced by such models use latent variables as axes representing missing (latent) predictors and the main axes of covariation across species (Hui et al., 2015;Warton et al., 2015). This model-based ordination technique incorporates the statistical properties of the data at hand. It provides an advantage over distance-based methods that only describe the data's multivariate geometry . The R package boral has different routines to fit pure latent variable models and extensions of multivariate generalized linear models (Hui, 2016). We fitted a pure-latent variable model based on binomial (Bernoulli) error distribution to produce a biplot based on the posterior median estimates of the latent variables. The model fit was restricted to two latent variables, and the derived ordination only showed the centroids of species with three or more occurrences. Star characters indicate the data source (see Section "Environmental Data"): *Data from CARS; **Data from Bio-Oracle.
We used the routine manyglm of the R package mvabund  to apply multivariate GLMs to analyze how environmental variables explain changes in species composition across seamounts and oceanic islands. This analysis uses multiple GLMs to model the variation of individual species and then generates a global model depicting environmental variables' effects on the assemblage. First, we fitted a multivariate GLM based on binomial distribution considering all environmental variables used in the global model of the univariate variables (i.e., longitude, latitude, depth, temperature, salinity, and dissolved oxygen and chlorophyll-a). Then we perform a manual model reduction by dropping predictors that did not improve model performance (expressed by AICc) using the function drop1. The significance of each predictor was tested by performing an analysis of deviance on the final multivariate GLM. We visually validated the results of the multivariate GLM by fitting a correlated response model to produce a partial ordination of residuals (residual ordination). This analysis consists of fitting GLMs to each species and including latent variables to account for residual correlation between species, for example, due to unmeasured covariates (Hui, 2016). By comparing the residual ordination with the unconstrained ordination, one can visualize how much of the variation in assemblage structure is explained by the predictors. Thus, the more the residual ordination reflects the patterns in the unconstrained ordination (i.e., based on pure latent variable modeling), the less the predictors used in modeling explain the variation in assemblage structure.

Environmental Variables
A marked longitudinal gradient in water mass properties is observed, with the mean thermocline (see 15 • C isotherm in Figure 2 and Supplementary Figure S2) sloping upward from west to east and outcropping near the continental coast and the minimum dissolved oxygen concentration occurring at ∼250-300 m and decreasing below 2 mL/L eastward from ∼87 • W at 26 • S to ∼81 • W at 33.5 • S. The lowest levels of dissolved oxygen (0.90 ± 0.03 mL/L) occur at stations within the North Nazca Ridge, which comprise the closest stations to the western limit of the OMZ. The average oxygen concentrations of the other zones (3.9 ± 1.2 mL/L) were more than four times higher than at the North Nazca Ridge. The highest oxygen concentrations (4.86 ± 0.16 mL/L) were at stations in the Easter Island region. The summary of the environmental conditions at each sampled area is shown in Table 1.

Taxonomic Composition of the Assemblage
For the target depth range of 50-400 m depth, our dataset included a total of 68 OTUs (operational taxonomic units) of echinoderms belonging to 35 families ( Table 2) and four classes. The classes Echinoidea and Asteroidea were represented by the highest number of species, accounting respectively for 47% (32) and 40% (27) of the 68 OTUs. Ophiuroidea accounted for 9% (6), and Holothuroidea for 4% (3) of the OTUs. No species of the class Crinoidea were observed. The classes Echinoidea and Asteroidea were predominant in the number of OTUs in all analyzed areas ( Table 2).
No species common to all the studied areas were found (Table 2 and Figure 3). The irregular urchin Scrippsechinus fisheri was the most widely distributed, occurring in four of the six studied areas (central Salas y Gómez Ridge, Salas y Gómez and Nazca Ridges junction, Desventuradas Islands, and Juan Fernández Ridge). Since most of the OTUs observed in the North Nazca area were poorly classified, we had to remove them from the analyses, resulting in the total separation of this area regarding the others (Table 2 and Figure 3). Easter Island only shared OTUs with the adjacent central Salas y Gómez Ridge area (Table 2 and Figure 3). The junction of the Salas y Gómez and Nazca Ridges contained the highest number of OTUs (Table 2) and shared 16 OTUs the central Salas y Gómez Ridge.

Species Richness
We did not identify any well-defined geographical pattern in the distribution of the taxa richness of echinoderm assemblages across the SEP (Figure 4A). The estimated values of taxa richness on sites ranged between 2 and 47, with the maximum observed at a station in the Salas y Gómez and Nazca Ridges junction. The GLM model results indicated that dissolved oxygen together with depth and salinity (Figures 4B-D) had a significant positive effect on echinoderm species richness ( Table 3).

Community Composition
The general patterns in the variation of the echinoderm assemblage structure across sampling stations are depicted in the pure latent variable ordination (Figure 5). The most notable    pattern is the clear separation of Easter Island stations from the other areas, indicating distinct echinoderm assemblages (Figures 5A,B). The differences among the remaining stations were obscured by the high variation in assemblages observed in the Salas y Gómez and Nazca ridges' junction. Stations from the North Nazca Ridge tended to differ from the Desventuradas Islands and Juan Fernández Ridge. Still, they all shared some similarities with the Salas y Gómez and Nazca ridges' junction zone (Figures 5A,B). The results of the multivariate GLM indicate that oxygen, longitude, chlorophyll-a explained changes in echinoderm assemblages across stations ( Table 4). The residual ordination derived from the correlated response model still showed some of the structure depicted in the unconstrained ordination (Supplementary Figure S3), suggesting that other unmeasured factors might also be driving the echinoderm assemblage structure.

DISCUSSION
Our results indicate that echinoderm species richness and community composition observed on mesophotic benthic ecosystems distributed across the SEP are shaped by gradients in oxygen, which is in line with our main prediction. Besides, echinoderm assemblages were also significantly affected by other environmental parameters, such as depth, salinity, longitude, and chlorophyll-a. Even though species richness did not present a notable spatial pattern, we detected a clear longitudinal gradient in community composition, explained by gradients in dissolved oxygen and chlorophyll-a. These changes in assemblage structure possibly reflect the sensitivity of echinoderm assemblages to the low oxygen zones located in the easternmost areas of our study area and the presence of potential barriers to dispersion for the westernmost assemblages. Considering that the available information on echinoderms of the SEP is extremely limited, the information reported in our results represents a significant addition to the understanding and management of this region of the Pacific Ocean.

Echinoderm Species Richness
Although we did not detect a clear gradient in echinoderm species richness across the seamounts and islands of the southeast Pacific, we observed that species richness had significant positive relationships with depth. Variables such as oxygen and salinity also affected taxa richness significantly, suggesting that these variables, or a combination of them, could act as environmental constraints for the establishment of diverse echinoderm assemblages in mesophotic ecosystems. A study assessing ophiuroid assemblages on seamounts of the southwest Pacific (O'Hara and Tittensor, 2010) found that environmental variables highly correlated with depth, such as water temperature range, significantly affected ophiuroid species richness across a much wider bathymetric range (i.e., 100-3,000 m depth). The GLM model was based on Poisson error distribution and log-link function.
They found that dissolved oxygen was a significant covariate to explain ophiuroid species richness, although their results indicated a negative relationship (O'Hara and Tittensor, 2010). This finding might suggest that the observed trends in relation to dissolved oxygen in our system might be restricted to euphotic and mesophotic zones, the zones in which the major clines in oceanographic conditions are observed (i.e., light incidence, temperature, density, and pelagic productivity). The positive effect of dissolved oxygen on taxa richness reflects a well-known response of the megafauna to this variable in deep-sea ecosystems (Levin et al., 2001). Dissolved oxygen can be a significant constraint for benthic organisms in marine ecosystems, especially within the OMZ (Levin and Gage, 1998;Breitburg et al., 2018). It can dictate zonation patterns in benthic communities, as reported by Sellanes et al. (2010) for the Chilean continental margin. Furthermore, the OMZ depth and geographic range experiences fluctuations on scales from intraseasonal to centennial that can be reflected in animal abundance and diversity (Sellanes et al., 2010). Apart from the physiological constraints on megafauna within the OMZs, certain groups of echinoderms can thrive on the edges of OMZs because of associated increases in food supply and decreased competition (Steckbauer et al., 2015).
In the SEP, the core of the OMZ is usually associated with the Equatorial Subsurface Water, which is characterized by low oxygen levels and high salinity (Silva and Neshyba, 1979). The southeast Pacific OMZ, as the result of the combined effect of a sluggish subsurface circulation (i.e., low ventilation) and a high export of organic matter from the productive coastal zone (Paulmier and Ruiz-Pino, 2009), is characterized by low dissolved oxygen concentrations (<1 mL/L) at intermediate depths (∼50-1,000 m). Its presence and western boundary are clear in the two main longitudinal sections of our study along which biological Chlorophyll-a 87 203.8 <0.001 The multivariate GLM was based on the binomial distribution, and the analysis of deviance was performed on 999 resampling iterations.
data were examined (Figure 2): at 26 • S (i.e., Salas y Gómez and Nazca Ridges), reaching 80 • W at 300 m depth; and at 33.5 • S (i.e., Juan Fernández Ridge), reaching 75 • W at 260 m depth. These tongues of low-oxygen waters are associated with high salinity; their southward extents are restricted by mixing with the low-salinity and well-ventilated Antarctic Intermediate Waters and Sub-Antarctic Waters (Silva and Neshyba, 1979;Llanillo, 2012). Although it is known that some echinoderms on the Chilean continental shelf present adaptations to hypoxic events (Steckbauer et al., 2015), insufficient information is available for the species observed in our study area. Besides, because of the scarcity of oxygen data in the region, the degree to which the OMZ influences echinoderm richness in the North Nazca Ridge, Desventuradas Islands, and junction of the Salas y Gómez and Nazca ridges is still uncertain. Thus, further studies are needed to provide insights into how hypoxia events limit echinoderm species richness in this region.
Considering that many deep marine systems are energetically dependent on the productivity in the euphotic zones, the spatial variation of pelagic productivity might significantly affect the productivity and diversity of deep benthic systems (Levin and Gage, 1998). Because our study area is characterized by a major gradient in pelagic productivity along the longitudinal gradient, some of the variation observed in echinoderm taxa richness might also be reflecting a response to increased food supply derived from primary production (Hobson et al., 1995;Lauerman and Kaufmann, 1998). Indeed, our dataset covers a longitudinal gradient in nutrient concentrations from the eutrophic waters associated with the Humboldt Current (Desventuradas Islands, Juan Fernández Ridge, and north Nazca Ridge) westward to the oligotrophic waters of the southern Pacific gyre (Easter Island and Salas y Gómez Ridge) (Thiel et al., 2007;Fernández and Hormazábal, 2014). Variation in productivity within the epipelagic zone driven by nutrient variability would affect deep-sea communities that depend on this allochthonous food source (Billett et al., 2001(Billett et al., , 2010. As the food supply for deep communities decreases, the system's capacity to support viable populations of multiple species also decreases, leading to communities with low species richness (Levin et al., 2001). The fact we failed to detect a significant effect of chlorophylla suggests that water masses, traceable by salinity and oxygen, could play a role in moderating the longitudinal effect of reduced surface productivity in overlying waters.
Although we observed that echinoderm species richness is relatively similar across areas, the number of reported species has been steadily increasing as new studies are published. For example, 14 new echinoderms records have been reported for Easter Island within the last 3 years , and five new species were observed during the two previous cruises (i.e., CIMAR 22 and EPIC) in the area of the Salas y Gómez Ridge and the Nazca Ridge junction (Mecho, personal obs.). Thus, we expect that the trends observed for species richness in the present study might be altered as new information is obtained.

Echinoderm Assemblage Composition
We observed a notable longitudinal change in the echinoderm assemblage composition westward of ∼86 • W, previously suggested by Parin et al. (1997) as a biogeographic break, and an abrupt change observed between stations at ∼101 • W and those near Easter Island at ∼107 • W. Echinoderm community composition was significantly affected by longitude, oxygen, and chlorophyll-a ( Table 4). Considering these findings, we suggest that a combination of environmental and external factors, such as barriers to dispersion (i.e., currents and the distance across shallow seamounts) and changes in habitat types (i.e., from sandy to rocky bottoms), might be driving species turnover. Potential barriers to dispersion within the Easter Island region have been suggested based on studies finding extremely low levels of connectivity of metapopulations of commercially important species (e.g., the cowry Monetaria caputdraconis, the spiny lobster Panulirus pascuensis and the sea chub Kyphosus elegans) between Easter Island and the Salas y Gómez Islet (located ∼389 km east of Easter Island) that support the presence of a potential hydrological barrier east of Easter Island (Meerhoff et al., 2018a,b). Seamount summit depths may also result in potential barriers between the Easter Island and the Central Salas y Gómez Ridge and Salas y Gómez and Nazca Ridges junction because the seamounts between these studied regions are deeper and thus have few summits with mesophotic benthic systems (i.e., depth < 300 m) that could act as stepping stones for some megafaunal species (Kim and Wessel, 2011).
Apart from geographic isolation, oceanographic processes associated with the South Pacific gyre area, which is characterized by higher temperatures and lower nutrient concentrations compared to zones near the Humboldt Current, may act as a barrier between Easter Island and the easternmost regions of the Salas y Gómez Ridge (Fell, 1974;Boyko, 2003). Thus, the contrasting echinoderm assemblages and high level of potentially endemic taxa observed along the Salas y Gómez Ridge  supports the hypothesis that Easter Island represents an independent ecoregion (the so-called Rapanuiian Province), at least at mesophotic depths (Schilder, 1965;Rehder, 1980;Parin et al., 1997;Boyko, 2003). Additional surveys and sample collections are necessary to determine the transition zones between these ecoregions and the geographic ranges and taxonomic identities of potentially endemic echinoderm taxa.
Species on Desventuradas Islands and Juan Fernández also present high levels of endemism (e.g., 77% of fish in the area are endemic; Friedlander et al., 2016), which is potentially reflected in our echinoderm samples. The species of these areas are probably adapted to colder water conditions and lower oxygen levels or specific hypoxic events (when the OMZ expands or when a westward propagating eddy with low-oxygen in its core passes by) (i.e., Scrippsechinus fisheri). They probably are also physically isolated from communities of South American coastal waters, given the relatively limited off-shore connectivity and short dispersal distances (Garavelli et al., 2014), which explains the differences observed between the communities on the oceanic islands compared to the fauna on the continental shelf (Friedlander et al., 2016). On the other hand, the net equatorward flux of the Humboldt Current System could be driving the connectivity between Juan Fernández Archipelago and Desventuradas Islands, as supported by three shared species of echinoderms (e.g., Parvulastra calcarata, Scrippsechinus fisheri, Brissopsis sp.). However, this number is lower than shared species between the other areas with a similar distance (i.e., Desventuradas Islands and Salas y Gómez/Nazca ridges, eight shared species).
Our analyses also depicted patterns in echinoderm assemblages that were coincident with gradients in oxygen, temperature, and depth (see Supplementary Figure S1) for Salas y Gómez Ridge and Desventuradas Islands. The latter was associated with shallower depths and higher temperature and oxygen (higher oxygen, but closer to the oxycline) so that small variation in depth can produce strong variations in oxygen. Stations from the Salas y Gómez Ridge and the westernmost part of the Salas y Gómez/Nazca ridges had deeper depths and lower oxygen levels than the Desventuradas Islands (see Supplementary Figure S1). For shallow-water communities, Desventuradas Islands have been considered part of the Juan Fernandez Ecoregion (Spalding et al., 2007). However, our data identified only 3 echinoderm species that are shared between both areas at similar depths and oxygen concentrations (Figure 2), suggesting that a reconsideration of these assignments could be done, at least for mesophotic communities.

Data Limitations
Our dataset comprises data from 12 different oceanographic expeditions that encompass a temporal window of 55 years. These expeditions rarely sampled the same zones of the SEP and used different sampling technologies. Given that different benthic habitats (e.g., rocky, pebble, gravel, sandy, and muddy bottoms) were unequally sampled, we admit that a considerable proportion of the observed variation in echinoderm community structure (especially within areas) might reflect this unaccounted habitat heterogeneity. Unfortunately, we could not incorporate bottom characteristics as a predictor in our analyses since this data was not available for a large proportion of our dataset (i.e., for those stations sampled with trawls; see Supplementary Table  S1). The unbalanced nature of our dataset also can introduce some bias in our results given the natural temporal variability that some ecological systems can present. Mesophotic assemblages can change markedly with local environmental conditions over decadal time scales (Smith et al., 2019). For instance, oceanic circulation changes associated with the El Niño-Southern Oscillation can affect pelagic productivity and marine food webs (Thiel et al., 2007). These changes can affect mesophotic assemblages both by thermal stress (Smith et al., 2019) and by altering the survivor rates of the meroplanktonic larvae as well as their post-settlement survivorship (Meerhoff et al., 2018a).
Although the caveats mentioned above suggest precaution when interpreting our results, we highlight that this pioneering study aims to present general variation trends of this rarely studied ecosystem as a first approach to elucidate the drivers of echinoderm communities. Thus, we expect that the patterns described here might motivate the development of standardized sampling efforts to assess how mesophotic echinoderms respond to gradients in oxygen and other variables in the SEP.

Implications for Hypothetical Climate Change Scenarios
Considering the potential sensitivity of echinoderm assemblage structure to dissolved oxygen suggested by our findings, we can hypothesize that climate alterations may influence benthic biota in the SEP's seamounts and islands system. Recent studies predict a significant decrease in dissolved oxygen under global warming conditions , and the current trend of ocean deoxygenation may lead to expansion of the geographic extent of OMZs and eventually lead to more extensive and persistent anoxic areas within OMZs (Stramma et al., 2008;Schmidtko et al., 2017). This extreme form of hypoxia or "low oxygen" can form dead zones on marine bottoms (Breitburg et al., 2018). Given the oxygen profiles described here, we would expect that stations near the Desventuradas Islands, North Nazca Ridge, the junction of the Salas y Gomez and Nazca ridges, and the Juan Fernandez Ridge could be further influenced by OMZ waters if it expands westward. Because some stations of the North Nazca Ridge are located only a few hundred kilometers westward from the current average location of the OMZ (Figure 2), it is expected that their assemblages would be affected in the short term (e.g., years to decades) by a reduction in the number of species or by changes in their component species as they experience increasingly more deoxygenated conditions. However, such a scenario does not consider the effects associated with changes in connectivity patterns due to changes in the circulation. For instance, the predicted increase in meridional winds off central Chile (Goubanova et al., 2011;Rykaczewski et al., 2015) is expected to intensify the net equatorward transport of the Humboldt Current System, potentially favoring the connectivity between Juan Fernandez and Desventuradas. Both factors (ocean deoxygenation and change in connectivity) may have compensating or cumulative effects on echinoderm assemblage structure.
Our study provides science-based information necessary for assessing the effectiveness of recently established marine protected areas within the Chilean EEZ, including the Motu Motiro Hiva Marine Park (∼150,000 km 2 ) created in 2010, the Nazca-Desventuradas Marine Park (∼300,000 km 2 ) created in 2016, the Rapa Nui Multiple Use Marine and Coastal Protected Area (∼740,000 km 2 ) created in 2018, and the Mar de Juan Fernández Marine Park (∼262,000 km 2 ) created in 2018. For example, our data reveal the extent of connections between areas, the environmental factors that affect these benthic communities, and therefore provide data for inferring the potential impacts on echinoderm communities associated with predicted environmental changes in those areas. Our data also suggest a potential revision of the Juan Fernandez Ecoregion mesophotic communities.

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/s.

AUTHOR CONTRIBUTIONS
AM, EE, and JS secured funding for participation in CIMAR22, EPIC, and Easter Island cruises and conducted fieldwork. AM developed the manuscript proposal and processed and identified specimens, and compiled the bibliographical data. BD and SG processed oceanographic data and edited environmental models. JG and AM processed biological data and conducted the statistical analysis. All authors contributed to the concept, writing, and editing of the manuscript.

ACKNOWLEDGMENTS
We would like to thank the Chilean Navy and the crew of the AGS Cabo de Hornos as well as the crew of the RV Mirai from JAMSTEC. We also thank Enrique Hey, Enrique Pate, Ricardo Hito, and Poky Tane Haoa from Rapa Nui for providing infrastructure and assistance during surveys around the island. Special thanks to Sergio Rapu and the Rapa Nui Heritage Foundation for providing the facilities for our on-island laboratory. We would also like to thank Dr. Tim O'Hara for his help with the ophiuroid classification. We would like to thank science teams and crew of the Cabo de Hornos and Mirai.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2021.574780/full#supplementary-material Residual ordination describing the variation in echinoderm assemblage data after controlling the effect of the environmental predictors (oxygen, salinity, depth, temperature, chlorophyll-a, latitude, and longitude). The more the patterns of (a) and (b) differ, the more the variation in echinoderm assemblages are explained by the environmental predictors.
Supplementary Table 1 | Detailed species occurrence data, including poorly classified OTUs that were excluded from the analyses.
Supplementary Table 3 | "enviro": averages of each environmental variable (salinity, chlorophyll-a, sea surface temperature, and temperature at 300 m depth) for each station. This data was used as predictors in the analyses.
Supplementary Table 4 | "dat_pca": CARS data including seasonal variation at each sampling station. This data was used to construct the PCA in Supplementary Figure 1.