Environmental Gradients and Spatial Patterns of Calanoid Copepods in the Southeast Pacific

Understanding the mechanisms maintaining the biodiversity of plankton communities in marine ecosystems subject to a strongly variable ocean has become a critical issue for modern oceanography. Here, we used data on distribution of calanoid copepods in the upper layer of the ocean (0–500 m), a widely distributed taxonomic group in the pelagic realm, to assess the effects of changing oceanographic conditions on their diversity patterns and family and species richness. Copepods abundance and occurrence were evaluated from 2002 to 2015 covering the region extended between the coastal upwelling zone (CUP-Z) and the offshore region of Chile at subtropical and temperate areas. We used spatial analyses of community structure descriptors, such as abundance and diversity (family and species richness), multivariate analysis and General Additive Models (GAMS) in order to study the effect of surface temperature and its gradients, mixed layer depth, salinity and Chlorophyll-a on copepod diversity. Seventeen families were identified comprising 151 species distributed in 3 predefined zones in the onshore-offshore gradient covering the coastal upwelling, the coastal transition and the oligotrophic zones, whereas over the alongshore gradient, same families were majorly linked to the northern and southern portions of the sampled area (20–40°S). Families and species were significantly structured over the zonal gradient, revealing the dominant habitat for each of the families. Spatial patterns revealed the presence of transitional zones comprised by mixed taxa. Over the alongshore gradient this transition zone was linked to the subtropical convergence (30°S). The spatial variation of sea surface temperature (SST) revealed strong environmental zonation of temperature gradients across onshore-offshore and alongshore dimensions. Mean SST combined with mean mixed layer depth explained about 40% and about 29% of variation in family and species richness, respectively over the onshore-offshore axis. We thus conclude that the environmental zonation imposed by SST and its spatial gradients, considered as ecological barriers, is the key driver for maintaining diversity of copepods in the southeast Pacific.


INTRODUCTION
The marine pelagic environment presents a high diversity of organisms as part of a complex trophic food web (Lenz, 2012). One of the main groups within the pelagic community is the zooplankton, mostly comprised by species with short life cycles (<1 year), and which is considered a key trophic link in the marine food web and. Zooplankton exhibit a high taxonomic diversity highly sensitive to environmental changes (Peijnenburg and Goetze, 2013;Pino-Pinuer et al., 2014), in particular their changes in abundance and distribution can be used as indicators of ecosystem response to climate driven variation (Beaugrand, 2003;Richardson, 2008).
Actually, the knowledge and understanding of oceanographic drivers influencing patterns of the large-scale distribution of zooplankton over large-scale domains, such as ocean basins is poor, limiting our capacity to predict changes in diversity and spatial distribution under a changing ocean. Biodiversity and species distribution of zooplankton are known to be affected by changes in oceanographic conditions, such as temperature, oxygenation, salinity and stratification, as it has been demonstrated in other regions of the world (Beaugrand et al., 2002;Beaugrand, 2003;Richardson and Schoeman, 2004;Peterson et al., 2006;Aronés et al., 2009;Gewin, 2010;Seibel, 2011), revealing the strong dependence of these organisms on oceanographic variables due to their limited migration capacity and so reflecting the strong effects of hydrographic patterns on their distribution and that of their prey (McClain and Barry, 2010). In this context, the Southeast Pacific (SEP) shows a great heterogeneity in hydrographic conditions from coastal to open sea gradient evidenced by a progressive increase of depth of the mixing layer, dissolved oxygen, temperature and surface salinity (Palma and Silva, 2006;González et al., 2019). Also, within the same gradient it is possible to observe drastic changes in biological productivity, which is high in the coastal zone (Humboldt Current) and extremely low in the ultraoligotrophic central gyre (Bonnet et al., 2008;Morel et al., 2010;Kletou and Hall-Spencer, 2012;Mann and Lazier, 2013). Across the alongshore gradient (72-110 • S), strong differences in temperature, salinity and levels of oxygenation of surface waters can occur. These changes are also accompanied by variable seasonal regimes from low to high latitudes (Strub et al., 1998).
The strong hydrographic variability observed in the SEP over both zonal and meridional gradients allows us to distinguish distinct ecoregions, where the environmental characteristics are seasonal, although relatively stable and homogeneous with respect to the adjacent ones, and such properties are also reflected in distinct communities between these ecoregions. For instance, along coastal and oceanic areas of the SEP five ecoregions have been suggested: Central Peru, Humboldtian, Central Chile, Araucanian, Juan Fernandez, and Desventuradas (Spalding et al., 2007). Although, a recent work based on ecological and biogeochemical division of zooplankton, proposed the existence of four regions over the zonal gradient, predefined in terms of surface chlorophyll-a as a proxy for productivity (González et al., 2019). Same work suggested that changes in the community structure of zooplankton in the SEP could be linked to the factors or processes that modulate the availability of nutrients at lower trophic levels, such as N-limitation and N-fixation (González et al., 2019). Nevertheless, it has now become clear that a changing temperature of the ocean is the major driver for altered diversity and distribution patterns of marine organisms (Antão et al., 2020). Changes in temperature and sea level have modified the geographical range of some copepod species (Gee, 1991;Crame, 1993;Molfino, 1994), also affecting their evolutionary history (Cronin and Schneider, 1990;Beaugrand et al., 2002). On a smaller scale, ecological factors influencing the physiology of the species can modify the spatial and temporal distribution of species, depending on their tolerance and optimal ranges of their populations, regulating diversity patterns. Such spatial and temporal variation, from microscale to large scale, has contributed synergistically to promote speciation events, and hence shaping up the actual patterns of biodiversity of pelagic organisms. However, the mechanisms that regulate biogeographical patterns need a better understanding to gain tools for predicting the variation in the species distribution under the effect of natural or anthropogenic impacts. Under this view, the pelagic ecosystem of the Southeast Pacific may constitute a suitable model to enhance our understanding of underlying mechanisms modulating zooplankton distribution patterns and other higher trophic levels. With these ultimate goals in mind, in this work we assessed the coastal-oceanic (zonal) and alongshore (meridional) distribution of pelagic Calanoid families and their species over pre-defined zones to test the hypothesis that ecological zonation is determined by environmental gradients of which the temperature gradient can play a unique role. We thus aimed at understanding the underlying mechanisms influencing the large-scale patterns of species distribution within ocean basins.

Study Area and Data Bases
The Southeast Pacific region off Chile has been the focus of various zooplankton studies and samplings over the last two decades carried out for our research team, although most of them in the coastal zone, also known as the southern portion of the Humboldt Current system covering the area between 18 • S and 40 • S. We have also conducted a few oceanographic cruises in the offshore region, even reaching the area near Easter Island (∼110 • W). Figure 1 illustrates our study area and the sites from which zooplankton samples have been obtained and used in this study and information of the samples sources and data are summarized in Supplementary Table S1. González et al. (2019) distinguished four ecological zones in this region, as based on the surface chlorophyll-a (Chla-a) concentration, from the highly productive upwelling area toward the central south Pacific gyre. This study found significant changes in zooplankton community structure among these zones, defined as: the CUP-Z, the coastal transition zone (CTZ), the oligotrophic zone (OLIGO-Z) and the ultra-oligotrophic zone (U-OLIGO-Z) which are illustrated in Figure 1. Here, we used these pre-defined zones to allocate all FIGURE 1 | The southeast Pacific region where the samplings were made between the years 2002-2015. The sampling stations are indicated by a red circle. The study area was divided in three zones according to surface chlorophyll-a values: coastal upwelling zone (CUP-Z), mesotrophic coastal transition zone (CTZ), and oligotrophic zone (OLIGO-Z) (González et al., 2019). The bathymetry of each zone is indicated by its color and isobars. the sampling sites, although the OLIGO-Z and the U-OLIGO-Z were considered only as OLIGO-Z, because of too few samples to separate them.
The sampling period for our study was extended from 2002 to 2015 with seasonal variable periods and with different sampling gears (see Supplementary Table S1), although in all cases the target community has been mesozooplankton (200-2000 µm) which is the zooplankton majorly composed by copepods in the ocean. Although most samples have been performed in the upper 200 m, some of them were obtained down to 500 m over OLIGO-Z zone. In the coastal zone a site is part of a coastal time series, so that it has a much better temporal resolution upon monthly sampling, in comparison with all the others for which in some cases a single sampling was done. Most of these data have been published as indicated in Supplementary Table S1 and the copepod data base has been uploaded in the South Pacific node of the Ocean Biogeographic Information Systems (OBIS) 1 .

Data Analyses
Data on species abundance from all sampling sites were used as number of individual m −3 . Species have been properly identified and only a few uncertain species are considered at the genera level. For each sampling point the sum of abundance for all species and families was made. However, since the source of samples varied with season, depth strata and sampling 1 https://obis.org/ gear, the effects of such factors on species abundances and composition were tested with variance analysis (ANOVA), firstly tested for normality and variance homogeneity by Shapiro-Wilks and Levene tests, respectively (see Supplementary Table S2). ANOVA test showed no seasonal or sampling gear effects, but significant differences between depth strata were found, so that data and analyses were treated separately into shallow (0-200 m) and deep (0-500 m) strata.
Changes in copepod families were evaluated for the different zones using community descriptors, such as family richness (R) which is the number of species per family, numerical abundance (N) and weighed frequency of occurrence (WFO) defined as: where NF is the number of stations in which the I family was present, NS is total number of sampled stations and REF is the number of species present in the I family. Also, we modified the equation of weighted mean depth suggested by Andersen et al. (2004), as to estimate the average weighted family abundance found across the meridional and zonal gradients of the SEP. For this, we estimated the weighted abundance across the zonal (ZWA) and meridional (MWA) gradients of the copepod families as: where Ni is standardized sum abundance to ind.m −3 for family i, Zi is the distance in km between a predefined starting point and end point of a given study area i and Ri is half the difference between these predefined boundaries of their i zone. Satellite data as on sea surface temperature (SST), mixed layer depth (MLD), sea surface salinity (SSS) and surface Chl-a were used as independent variables during the study period (August 2002until October 2015. Mean values of these variables were calculated for each sampling station from monthly satellite data. Also, the spatial gradients of SST, both zonal and alongshore were estimated for these mean values for the entire study area. These spatial gradients represent the either alongshore (meridional) or onshore-offshore (zonal) difference in temperature between two consecutive SST means at 1 degree (110 km) spatial resolution. MLD and SSS were obtained from the ARGO floats 2 with a spatial resolution of 3 degrees, while surface Chl-a was downloaded from the MODIS satellite 3 with a spatial resolution of 0.8 degrees. SST was obtained from OAFlux project 4 with a spatial resolution of 1 degree.

Statistical Analysis
We assessed similarities in distribution of Calanoid copepods families and species only for the shallow stratum (0-200 m), within and between zones, in terms of abundance (individuals m −3 ) across a zonal transect. This zonal gradient could not be assessed for the deep stratum (0-500 m), because most samples were only within oceanic region. The analysis of the alongshore distribution of families and species was also restricted to north-south comparisons with respect to the subtropical convergence (ca. 30 • S), because of low coverage of the whole gradient and unbalanced design, due to a large number of samples for the southern area 36 • S and few of them in the northern portion of the SEP. The structure of the community was evaluated by multivariate analysis performed with PRIMER v.7 (Anderson, 2008). For this, we grouped the different samples over the zonal gradient in according to their location, as CUP-Z, CTZ, and OLIGO-Z. For evaluating similarity in composition and abundance between and within zones, we applied a multidimensional scaling (MDS) and analysis of similarities (ANOSIM) with log-transformed data [Log (x + 1)] for families and presence/absence for species. We used presence and absence data for species, because in some cases only genera were reported, also to reduce extremely high variance of abundances due to presence of dominant and rare species. The matrices of similarities were built for families from the Bray-Curtis index and for species with Jaccard distance. The differences or similarities between the zones were evaluated with SIMPER (percentages of similarity). Also, R Statistic analysis was applied to test differences between groups.
Biota-environmental (BIOENV), Distance-based Linear Model (DISTLM), and Generalized Additive Models (GAM s ) were used to correlate the copepod community with environmental. For this, the abundance of species and families were used as response variables for DISTLM, while for GAM s and BIOENV richness values were used as dependent variable in order to determine variables to be included in the models, a correlation analysis using Spearman's test and variance inflation factor (VIF), among the hydrographic variables was first carried out. SSS was found to have strong correlations (r > 0.60 and VIF > 5) with SST, Chl-a and MLD, and it was therefore removed from the analysis in order to minimize collinearity among independent variables. In the BIOENV analysis, similarity matrices were built between log-normalized environmental variables and both family and species richness using Euclidean distance, delivering as result the Sperman rank correlation between the two variables. An automated, stepwise procedure were used both DISTLM and GAM S analyzes to identify influential variables on families and species richness, and the final model for each diversity measure were obtained according to the Akaike information criterion (AIC). Smaller values of the AIC indicate a better model and yield into the "best" combination of environmental predictor variables that explain the largest amount of variation in the response variable. For DISTLM analysis, similarity matrices were built between the abundance of species and families with the log-normalized environmental variables using Bray-Curtis distance, delivering as result a "Pseudo-F statistic" for each environmental variable. Both the BIOENV and DISTLM analyzes were performed with PRIMER v.7 (Anderson, 2008). Finally, for the GAM S analysis percent deviance explained (DE) and determination coefficient (R 2 ) were calculated for each model to examine the overall goodness-of-fit. The models were constructed using a Gaussian distribution with an identity-link function and smoothing function for oceanographic factors through the "mgcv" package (Wood, 2017) of R (Team, 2000). Supplementary  Table S3: Species list. From the whole data base, 17 Calanoid families were found for a total number of 151 species (see Supplementary Table S3). Information on the specific zone where the species were found is shown in Supplementary Material (Supplementary Tables S4, S5 for zonal and meridional locations, respectively). The family Aetideidae was the most represented with 19 species, followed by Scolecitrichidae with 17 species and then Clausocalanidae with 16 species. The least represented families were Rhyncalanidae and Temoridae with 2 species. Paracalanidae, Acartidae, and Calanidae were the most abundant and common species, mainly present within the CUP-Z. Rare families, in terms of occurrence and abundance, were Temoridae and Calocalanidae, both present only in the CUP-Z and OLIGO-Z, respectively, with low abundances ( Table 1).

Copepod families and species are listed in
When examining the community descriptors, abundance (N), species richness (R) and weighed frequency of occurrence (WFO) for each predefined zone for the 0-200 depth strata (Table 1), it was found that much greater abundances occurred in the CUP-Z, while in the OLIGO-Z abundances were commonly 1 order of magnitude lower when compared same families. The number

Zones
Family Richness (R), mean abundance (N) with its standard deviation (±SD) and weighed frequency of occurrence (WFO) were estimated as divided by zone. Oceanographic zones are defined as: coastal upwelling zone (CUP-Z), mesotrophic coastal transition zone (CTZ) and oligotrophic zone (OLIGO-Z). For the three study areas, sampling was done between 0 and 200 m (Shallow), except for the OLIGO-Z zone, which was also done at 0-500 m (Deep).
of species per family (family richness) tended to decrease in the OLIGO-Z in some families containing most species. This was the case for Aetideidae (14 species) and Scolecitrichidae (11 species). In most cases families had more species in the CUP-Z. However, when comparing the two depth strata within the OLIGO-Z zone, we observed greater abundance and richness at the deep stratum. This was observed in Calanidae (12 species), Aetideidae (12 species) and Paracalanidae (12 species) ( Table 1).

Onshore-Offshore Distribution
For the zonal distribution analysis, we used the zonal weighed abundance (ZWA) of the families at shallow depths (0-200 m) to assess their dominant habitat across the onshore-offshore gradient. This distributional pattern is shown in Figure 2. It was found that many families can occupy the entire zonal range, although a few appeared restricted to the coastal area. It was also found that ZWA values were aggregated within the first 2000 km from the coast, and these varied considerably among families, although with a tendency toward the coastal zone. In any case, ZWA's could help us to allocate dominant habitat of families, such that Acartidae, Augaptilidae, Candaciidae, and Lucicutidae which seemed to have more oceanic habitats compared to the abundant Eucalanidae, Rhincalanidae, and Temoridae, whose habitats seem centered close to shore (Figure 2). Some dominant families in the CTZ were Aetideidae, Heterohabbidae, Metridinidae, Phaennidae, and Scolecitrichidae.
In terms of absolute abundance across the onshore-offshore gradient, Figure 3 represents how families are distributed over this axis. Calanidae and Paracalanidae were the families most widely distributed (with more homogeneous distribution) across this zonal gradient and along with Acartidae and Centropagidae they were the most abundant within the CUP-Z. By contrast, Temoridae, Aetideidae, and Phaennidae were represented with much more restricted distributions (Figure 3).
Changes in community structure at the family level across the onshore-offshore gradient became evident after NMDS analysis using the Bray-Curtis similarity index of distance. This NMDS analysis ( Figure 4A) shows that a large cluster can be arranged at 20% of similarity with most samples and significant grouping may occur at 40% of similarity where the samples of the CUP-Z zone are grouped with the others. Also, ANOSIM analysis showed significant differences in the abundance of families between the zones OLIGO-Z with CUP-Z (p-value = 0.001; Rglobal = 0.75) and CTZ (p-value 0.001; Rglobal = 0.69), but not between the zones CUP-Z and CTZ (p-value 0.14; Rglobal = 0.13).
A similar analysis to that of families was applied to the community structure at the species level. In this case however, using presence or absence since abundance was too heterogeneous for comparisons. The Jaccard index of distance yielded a NMDS shown in Figure 4B. NMDS analysis revealed that at low level of similarity (10%) three groups could be observed in which, as in families, the samples from the CUP-Z zone became grouped with those from other zones. However, ANOSIM analysis at the species level showed a result similar as that observed at the family level. Significant differences were found between zones OLIGO-Z with CUP-Z (p-value 0.001; Rglobal = 0.51) and CTZ (p-value 0.001; Rglobal = 1.00). On the other hand, presence and absence of species was different between the shallow (0-200 m) and deep (0-500 m) layers

Alongshore Distribution
When looking at the alongshore distribution of copepod families at the shallow depth stratum (0-200 m), we used a meridional weighed abundance (MWA), designed in according to our sampling area. This allowed us to allocate the dominant habitat of each family over the latitudinal axis within the CUP-Z. This alongshore pattern is shown in Figure 5. It was found that most families were widely distributed over the region, except for Rhincalanidae mostly inhabiting the southern area. Two major groups can be distinguished located either in the southern or northern portions of the study area. The southern area was characterized by seven families, Augaptilidae, Candaciidae, Centropagidae, Eucalanidae, Scolecitrichidae, and Temoridae which were the most abundant, whereas in the northern portion 10 families Calanidae, Clausocalanidae, Euchaetidae, Heterorhabditidae, Lucicutidae, Metridinidae, and Paracalanidae were well represented (Figure 5). A shade plot (Figure 6) shows that Paracalanidae and Calanidae were widely distributed families along the upwelling zone, also representing the most abundant copepods in the region. In a lesser extent, Acartidae and Metrinidae also appeared as widely distributed, whereas other families exhibited more restricted distributions. The northern area, two families Eucalanidae (4 species) and Centropagidae (2 species) showed the highest abundances, whereas in the southern portion three families Paracalanidae (6 species), Calanidae (6 species) and Acartidae (2) presented the highest abundance values.

Environmental Correlates
In order to examine environmental correlates of the copepod community with habitat conditions, we used the grand mean values of the study period for the different oceanographic variables (independent variables). The spatial patterns of these variables (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) are shown in Figure 7. In the CUP-Z a warmer SST in the range of 16-21 • C prevails in the northern portion, while a range of 12-16 • C prevails in the southern portion of this CUP-Z. The CTZ exhibits a similar northernsouthern pattern and SST range as that of the CUP-Z, whereas the OLIGO-Z shows a gradual increase in SST from about 21 • C up to 27 • C at the northeast portion of the study area ( Figure 7A). Sea surface salinity (SSS) exhibited a high variation, both over the zonal and alongshore gradients. In the CUP-Z a higher SSS in the range of 34.5-35 predominated in the northern zone, compared to a range of 33-34 observed in the southern region. The CTZ showed a similar SSS northern-southern pattern as that exhibited by the CUP-Z (Figure 7B). Contrasting with SST and SSS, the mixed layer depth MLD showed a strong variation only over the zonal gradient. In the CUP-Z the shallowest depths of MLD (<20 m) were observed, while toward the OLIGO-Z zone it shows a gradual deepening mainly in the northern portion (>80 m) ( Figure 7C). Finally, Chl-a exhibited an abrupt increase of its concentration only over the zonal dimension. The highest Chl-a concentrations were observed in the CUP-Z zone (>1 mg m −3 ), while toward the OLIGO-Z zone the concentration was extremely low (<0.05 mg m −3 ) (Figure 7D).
Spatial patterns of SST gradients over the zonal and alongshore axes are shown in Figure 8. The Southeast Pacific showed strong temperature gradients within and across all predefined regions for the study period. The zonal gradient exhibited a very particular spatial pattern, characterized by bands nearly parallel to the coast of about 700-900 km width. These features showed stronger gradients (>0.4 • C), in the region between 80 • and 70 • W, corresponding mainly to the CTZ off the coastal upwelling area ( Figure 8A). The alongshore gradient presented a zone of lower variability (<0.2 • C) between 20-30 • , while toward the southern portion an intensification of this gradient was observed and the differentiation between these two patterns became evident at 30 • S corresponding to the subtropical convergence ( Figure 8B).
The multivariate analysis between the abiotic and biotic variables (families and species) was performed with the family and species matrices. In both cases similar results were found. The BIOENV analysis between oceanographic variables and the abundances of families and species showed a high correlation with SST (families: r = 0.5 and species: r = 0.6) ( Supplementary  Table S6). Indeed, DISTLM analysis indicated that the three variables SST, MLD and Chl-a correlate significantly with the abundance of both species and families (p-value < 0.05), although SST explained the highest percentage of variability in the community matrix (families: 21% and species: 20%). Furthermore, the sequential test with the stepwise procedure selected the SST and MLD as the variables that best explain the variation in the species matrix, for the selection of the best model, the lowest value of the AIC was used. More detailed results from this analysis are shown in Supplementary Tables S7, S8. FIGURE 4 | Non-metric multidimensional scaling diagram (NMDS) of families and species within the order Calanoida found in the 0-200 m layer, across an onshore-offshore gradient, divided in zone which were defined as: coastal upwelling zone (CUP-Z), mesotrophic coastal transition zone (CTZ) and oligotrophic zone (OLIGO-Z) of the southeast Pacific. For families NMDS (A) was calculated using its abundance through the Bray-Curtis distance index from the transformation of data with log (x + 1). For species NMDS (B) was calculated using the Jaccard distance index from the transformation of the abundance data as presence/absence of the species.
When looking for a predictive model for copepod diversity, it was found that most oceanographic variables had a significant effect on family and species richness at the shallow depth stratum, as resulting from application of a Generalized Additive Models (GAMs). Over the zonal gradient we found from the AIC value, as the most parsimonious criteria for best-fit GAMs, that a predictive model incorporated SST and MLD to explain both the family and species richness ( Table 2). Relationships between  family richness with SST and MLD explained together a 39.80% of the variance ( Table 2). These effects were highly significant (P-value < 0.01) and showed a positive relationship with the SST, and a negative relationship with MLD (Figures 9A,B). At species level, it was observed that variation in richness was explained in 29.10% by both SST and MLD (Table 2). However, only SST exhibited a significant and positive effect on richness (Figures 10A,B).

DISCUSSION
Biological diversity of zooplankton can be modulated by environmental gradients which modify the spatial patterns of species and communities (Lawton, 1999;Hillebrand and Azovsky, 2001;Rex et al., 2001). Although the partition of the habitat over the horizontal plane of zooplankton communities has been reported in some works (e.g., Heinrich, 1973; FIGURE 8 | Spatial gradients of sea surface temperature (SST) for the onshore-offshore (Zonal) (A) and the alongshore (meridional) (B) gradients in the southeast Pacific for the period 2002-2015. Gradient resolution is for 1 • from mean monthly values of SST. Gonzalez and Marín, 1998;Linacre and Palma, 2004), there have been difficulties for identifying and understanding the mechanisms that modulate differential distributions. In this context, our analyses revealed a high structuring of the order Calanoida, linked to spatial zonation, both in zonal and meridional gradient. However, some sources of variance may introduce biases affecting the observed species distributions. For example, sampling depth can modify the observed horizontal distribution of copepods (e.g., Judkins, 1980), and even though we separated samplings done Means values of sea surface temperature (SST), Chlorophyll-a (Chl-a) and mixed layer (MLD) were used as predicting variables. Akaike information criterion (AIC), determination coefficient (R 2 ) and deviance explained (DE) are shown.
over different depth strata, variability introduced by vertical distribution, including vertical migrations become difficult to account in our horizontal analyses. Also, most of the samplings were made in the CUP-Z targeting the upper 200 m, which is the zone where most species of copepods aggregate due to vertical restrictions imposed by presence of a shallow oxygen minimum zone (Escribano et al., 2009). Another additional source of variation could be caused by using different sampling gear. In our case, we selected the same mesh size for the zooplankton nets, but different opening diameter of nets and variable towing speed may have occurred. For instance, a comparison of two type of gear, the MOCNESS and WP2 nets, showed that estimates of biomasses differed between these nets, although such differences were mostly accounted by the largest fractions of zooplankton represented by euphausiids and amphipods, and there was no evidence for a significant effect of net opening diameter in the estimates of abundance of smaller organisms, such as copepods (Gjøsaeter et al., 2000). In our analyses we did not find significant differences in copepods abundances estimated by two different sampling gears, Tucker trawl and Multinet and it seems that both devices can properly sample mesozooplankton which was our target group. This result agrees with a previous study in which different vertical, oblique, and multiple opening/closing net systems resulted in similar estimates of zooplankton with comparable mesh-sized nets and a sufficiently high ratio of mesh to mouth opening (Skjoldal et al., 2013). Samples analyzed in this study were obtained in different years and seasons. The coastal transition and oligotrophic zones were sampled during spring-summer, unlike the coastal region (upwelling zone) which was sampled throughout the year, such that it was then possible to assess the effect of seasonality on the species composition and families. Even although, for the coastal zone it has been reported an effect of seasonal variability on the composition of species (Hidalgo and Escribano, 2001;Escribano et al., 2007), we did not detect such effect when evaluating composition of families. This lack of seasonality effect could be explained by continued presence of copepod families, although eventually represented by different species, or different stages of same species, so that showing a lack of variability in the total abundance and diversity at different times of the year. Continuous reproduction and multigenerational populations can also help explain species occurrence year-round (Escribano and McLaren, 1999;Hidalgo and Escribano, 2007).
With respect to zonation over the zonal gradient, as found in the composition of families and species of the order Calanoida, we found highly significant changes from the coastal upwelling to the oligotrophic zones. Other works in the same area also reported changes in species composition; example for euphausiids (Riquelme-Bugueño et al., 2012) and (Palma and Silva, 2006). One of the main explanations for this pattern is attributed to the presence of different water masses, although variable sources of food and nutrient over this gradient has also been suggested, providing distinctive characteristics to each of the study zones (Silva et al., 2009;González et al., 2019). The costal upwelling zone is mainly influenced by Equatorial Subsurface Water, characterized by a high concentration of nutrients and source of fresh carbon from the atmosphere that is rapidly transferred to copepods (Morales et al., 1996;Gruber et al., 1999;Daneri et al., 2000;Silva and Valdenegro, 2003;Mompeán et al., 2013).However, mesotrophic coastal transition and oligotrophic zones are composed by Sub Antarctic Water and Subtropical Water, respectively, which present lower concentrations of nutrients and more regenerated carbon sources compared to the coastal zone (Reid, 1973;Nakatsuka et al., 1992). All these factors along with presence of the anticyclonic gyre in oceanic region can isolates species assemblages from the adjacent areas, and so promoting a semi-closed ecosystem regulated by in situ processes (McGowan, 1971(McGowan, , 1977. Over a longer time-scale samples of marine microfossils collected during the Pleistocene indicate that the longitudinal zonation observed in Southern Pacific persisted allowing to observe the structuring patterns at the level of family and species of calanoid copepods (Riedel and Funnell, 1964).
A high degree of structuring can also be observed over the alongshore gradient, with spatial agreement between analyses done at the levels of both species and families, showing the presence of three groups distributed in the north, a transition region, and the south. Other works done about the same gradient have also reported the presence of these biogeographic breaks in other organisms such as kelps, and bivalves mollusks (Broitman et al., 2001;Cárdenas et al., 2009;González et al., 2012). The boundaries of these biogeographic regions are not fixed and can vary seasonally and between years (Escribano et al., 2003). Despite that the limits between the different regions are not clear, it has been reported that the transition region shows a high correlation with a gradient of surface temperature (Broitman et al., 2001), irregular patterns of circulation (Hormazabal et al., 2004), and a sharp narrowing of the continental shelf (Strub et al., 1998). Furthermore, the presence and interaction of different water masses, such as: equatorial subsurface water and intermediate Antarctic water, can be the source for equatorial species in the north region and polar/sub-Antarctic in the south, possibly generating a transition zone between 30 • S corresponding to tropical (Escribano et al., 2003;Hidalgo et al., 2010).
In this study we report a strong dependence of richness with temperature both at family and species levels. In this respect, it is known that temperature can modify the rates of various biological processes in copepods, such as growth, productivity and mortality, but also influencing physical conditions of the water column, such as stratification and availability of nutrients and indirectly affecting the ecology of copepods (Reese et al., 2005;Kiørboe and Hirst, 2008). Along the zonal gradient it was observed a notorious increase in mean temperature from the coastal zone to the open ocean associated with increasing richness of copepod families. In fact, mean temperature was strongly correlated (ca. 29-39%) with richness both of families and species levels and a positive correlation with temperatures values > 17 • C due to the presence of Subtropical Water (Silva et al., 2009). In several regions of the world a positive correlation between temperature and richness of aquatic zooplankton species has been reported, accompanied by a decrease in abundance (Matsubara, 1993;Castro et al., 2007;Hessen et al., 2007) and with marine organisms at a global scale as well (Antão et al., 2020). Regarding the alongshore gradient we observed that temperature showed a significant effect on species richness and this can be associated with rapid and abrupt change of temperature between 20 and 30 • S in the coastal zone due to a spatial variation in the seasonal upwelling regime between northern Chile and southern Chile (Escribano et al., 2004), and presence of different water masses (Silva et al., 2009). The onshore-offshore (zonal) gradient of SST indeed revealed the existence of discrete habitats imposing constraints for population connectivity of epipelagic plankton. These physical barriers can play a significant role for evolutionary processes and gene flow across this zonal gradient, as shown by a recent work revealing ongoing speciation in the cosmopolitan copepod Pleuromamma abdominalis over this gradient (González et al., 2020). Such processes can certainly be extended to the community level and so influencing both family and species richness.
In addition to temperature, many other factors may influence zooplankton diversity patterns (Angel, 1997 for review). For instance, in the Benguela biome it has been suggested that variation in primary production (PP) and the coupling between primary producer and secondary producer may be the main controllers of the copepod community structure (Woodd-Walker et al., 2002). In our work, Chl-a as a proxy of food or ocean productivity, did not show any correlation with copepods diversity and species distribution. In this regard, the spatial scale is an important issue to consider and it may be that local (patchiness) effects of food resources can indeed influence species distribution and composition. However, over a large, basin scale physical barriers, such as those formed by temperature gradients, and the temperature within them are key drivers for diversity patterns. Depth of the mixed layer also deserves further attention. In the coastal zone vertical distribution, commonly restricted by strong gradients in temperature, oxygen, light or other factors, the available habitat becomes spatially limited, and deepening of the mixed layer which also correlates positively with temperature may help expanding the vertical and so promoting increased diversity (Tittensor et al., 2010;Mora et al., 2011).
In any case, temperature and its spatial arrangements have been reported as a major correlate of zooplankton diversity in other studies as well (Rombouts et al., 2009;Beaugrand and Kirby, 2010), implying that (Hessen et al., 2007) ongoing warming of the ocean may greatly impact copepod diversity patterns by either altering their comfort ranges of temperature or by modifying the spatial gradients.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study have been deposited in the Southeast Regional Node of OBIS [https://ron.udec.cl/].

AUTHOR CONTRIBUTIONS
All co-authors have contributed to the work. CG contributed to data analysis and writing. RE integrated the all results and developed the direction of the manuscript. JM-M and RE participated in discussion and data analyses and commented on the manuscript. All the authors contributed to the article and approved the submitted version.

FUNDING
This work was funded by Grants CIMAR-21 and CIMAR-22 of CONA-Chile and the Millennium Institute of Oceanography (ICN12_019-IMO). Additional funding was provided by Grant FONDECYT 188-1682. CG work was supported by the CONICYT Scholarship No 21160714.