Samples were collected between February 2 and March 21 during the austral summer of 2018 off the northern coast of Chile, on board the German FS Sonne (cruise So261). The stations were located above the Atacama Trench (stations 2 and 4), between 21° 47’ S and 23° 21’ S ( Figure 1A ). Hydrographic data, including temperature, salinity, oxygen concentration, and pressure, were obtained using a Seabird SBE-9 plus CTD equipped with an oxygen sensor, coupled to an oceanographic rosette with 10–24 L Niskin bottles. Once on deck, the Niskin bottles were emptied and pre-filtered using a 100 µm nylon mesh. Water samples were obtained at several discrete depths from the two sites for analysis of particulate organic carbon (POC), particulate organic nitrogen (PON), carbon and nitrogen (C/N) ratios, and isotope signatures of carbon and nitrogen (δ¹³C and δ¹⁵N (‰)). Between 0.5-60 L of seawater, depending on the POM concentration, were gently filtered through pre-combusted (450°C, 24 h) GF/F filters (Whatman, 47 mm diameter), which were immediately dried at 60°C for 24 h. The filters were subsequently wrapped in aluminum foil and stored in a sealed container with silica gel to avoid moisture until analysis in the Laboratory of Biogeochemistry and Applied Stable Isotopes (Pontifical Catholic University of Chile). The filters were placed overnight in a desiccator saturated with HCl fumes to remove inorganic carbonate, packed in tin capsules, and fed via flash combustion into a Thermo Scientific Flash 2000 CHN elemental analyzer (EA) coupled with a Thermo DeltaV Advantage isotope ratio mass spectrometer (IRMS). Analytical errors in the elemental determinations were 3 μg and 7 μg for nitrogen and carbon, respectively, according to linear regression against the standard (acetanilide). The precision of the stable isotope analyses, as determined using four internal standards (acetanilide, atropine, caffeine, and glutamic acid) in triplicate, was 0.33‰ for δ¹³C and 0.21‰ for δ¹⁵N.

Zooplankton samples for metabarcoding analysis were collected at two stations (details are provided in Supplementary Table S1 ). Sampling was conducted between 0 and 5000 m using a multiple opening/closing net and an environmental sensing system (MOCNESS-10) with a 10 m² mouth opening and five 333-μm mesh nets, equipped with pressure and flowmeter sensors. Five nets were closed every 1000 m. Zooplankton samples were split in half with a Motoda splitter (Motoda, 1959) and sieved using a 200-μm mesh sieve to remove excess water. Half of the samples were preserved in 99% ethanol, which was replaced after 24 hours. The samples were stored at -20°C for subsequent DNA analysis in the laboratory.

The extraction method followed that described by Blanco-Bercial (2020). In short, this is a modification of the E.Z.N.A.^® Mollusc DNA Kit (Omega Bio-Tek, Norcross, GA, United States) adapted for large volumes, using Sodium Dodecyl Sulfate (SDS) lysis buffer instead of the ML1 buffer from OMEGA. Two pseudoreplicates were obtained from each sample from each stratum. DNA was quantified for each sample with a Qubit 4.0^® fluorometer, reporting concentrations of 6.8 and 73.8 ng μL^-1.

One PCR reaction and sequencing were performed for each pseudoreplicate per sample for the two different genetic regions at the AUSTRAL-omics laboratory of the Universidad Austral de Chile. Amplification of the mitochondrial COI gene was performed using the primers mlCOIintF and jgHCO2198 (Geller et al., 2013; Leray et al., 2013), whereas primers 1389F and 1510R (Amaral-Zettler et al., 2009) were used for the nuclear V9 region of the 18S gene. Amplifications were performed independently for each amplicon in a final volume of 12 μL, which included 6μL of AmpliTaq Gold enzyme (0.3 U μL^-1), 2 μL of ultrapure water (HyClone), and 1 μL each of forward and reverse primers at a concentration of 2.5 μM each, and 2 μL of template. The amplification protocols for the two amplicons are shown in Supplementary Table S2 . PCR success was confirmed on 1.5% agarose gel in 1x TAE Buffer. The PCR products were purified using magnetic beads in 10 mM Tris buffer.

Libraries were created according to a protocol based on “Illumina 16S Metagenomic Sequencing Library Preparation” (Illumina, 2015) for both markers. In this step, the index primers and the purified amplicon were used in a second PCR reaction with a master mix composed of 25 μL of 2x KAPA HiFi HotStart ReadyMix, 5 μL of each Index Illumina primer, 5 µL of purified PCR product, and 10 µL of PCR-grade water (HyClone), for a total volume of 50 µL. The PCR protocol followed to amplify the libraries included one cycle of initial denaturation at 95°C for 3 min, followed by eight cycles of denaturation at 95°C for 30 s, alignment at 55°C for 30 s, elongation at 72°C for 30 s, and one cycle of final elongation at 72°C for 5 min. The indexed PCR product was purified using AMPure XP magnetic beads, with a final elution volume of 25 μL. Subsequently, the libraries were checked by capillary electrophoresis using an Agilent Fragment Analyzer with the DNF-910 Kit. The libraries were quantified by fluorometry using Qubit 4.0^®. Bidirectional sequencing was performed, including negative controls and samples, on an Illumina MiSeq using the Reagent Kit v3 (600-cycles), spiked with a minimum of 20% PhiX. The sequences obtained are available on the NCBI website under the bioproject: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA908823.

The demultiplexed COI and 18S gene samples were analyzed using MOTHUR ver. 1.43 (Schloss et al., 2009) on the computing cluster of the Instituto Milenio de Oceanografía (IMO) (https://www.imo-chile.com). The full script with annotations is available at https://github.com/carolinagonzaleze/SONNE_above_Atacama_Trench. Forward and reverse sequences were compared nucleotide by nucleotide, and mismatched positions were kept ambiguous if the quality difference between both strands was less than 10. If a base was compared with a gap in the other fragment, it was maintained only if Q > 30. Subsequently, all contigs with ambiguous nucleotides, more than 10 homopolymers, and fewer than 250 base pairs for COI and 115 for 18S were removed. The sequences obtained from the 18S gene were aligned to the v9 region of the SILVA 18S database (Quast et al., 2012). The sequences were trimmed to the length of the V9 region and all incomplete sequences were removed. For both genes, chimeras were detected using VSEARCH (Rognes et al., 2016) and MOTHUR. Errors associated with high-throughput sequencing were removed by UNOISE (Edgar and Flyvbjerg, 2015) using MOTHUR (diffs = 1) for 18S, obtaining amplicon sequence variants (ASVs) that were 100% similar. The COI gene sequences were clustered to 95% similarity, which is the most common threshold for COI species ( Supplementary Figure S1 ) (Bucklin et al., 2021, Bucklin et al., 2010). Potentially erroneous operational taxonomic units (OTUs) were further removed with the “LULU” package in R (Frøslev et al., 2017; R Core Team A, 2020).

Taxonomic assignment for 18S V9 ASVs was performed using the Silva 128 ribosomal database with a naïve Bayesian classifier algorithm implemented in MOTHUR. For COI OTUs, the assignment was done using BLASTN (Zhang et al., 2004) against the GenBank nucleotide database (Benson et al., 2005), comparing the first ten matches and selecting assignments with e-values lower than 10^-6 and identity values greater than 97% (species) or 85% (family or order). Only sequences assigned to metazoans were selected for further analysis of both markers.

Pseudoreplicates from each sample were summed for each station, allowing for greater sequencing depth of the study genes. All samples were standardized using the minimum number of reads per sample:220,984 for ASVs and 109,358 for OTUs. These standardized values are hereafter referred to as sequence abundances. The resulting ASVs/OTUs with total abundances below two were eliminated (global singletons and doubletons). To analyze the effect of the number of reads per sample on diversity, rarefaction curves were performed in the iNEXT (iNterpolation and EXTrapolation) package in R (Hsieh et al., 2016; R Core Team A, 2020), with confidence values obtained after 1000 bootstraps. Diversity indices, such as the total number of taxa (S), Shannon diversity index (H’), Simpson’s inverse (1-λ), and Pielou’s evenness index (J’), were calculated using PRIMER ver. 7 (Clarke and Gorley, 2015). Additionally, Phylogenetic diversity (PD) was calculated from a tree generated using ClearCut, as implemented in MOTHUR (Faith, 1992).

Similarity matrices between samples were constructed using the Bray–Curtis distance after Hellinger transformation (Legendre and Gallagher, 2001) and plotted using Principal Coordinate Analysis (PCoA). Clustering between samples from different strata and stations was performed using a similarity profile routine (SIMPROF test) (Clarke and Gorley, 2015). Differences in the community structures of the generated clusters were compared using the analysis of similarities (ANOSIM) procedure with 10000 permutations for significance testing in PRIMER ver. 7 (Clarke and Gorley, 2015).

Considering the strong gradients in near-surface waters and for statistical purposes, environmental data (including POC, PON, C/N, and their isotopic signatures) from the upper 1000 m were assumed as well represented by mid-layer values (500 m), whereas average values were estimated for deeper layers that had much weaker gradients. A summary statistic for all the environmental data can be found in Table S3 of the Supplementary Material .

The relationship between the abundance of ASVs/OTUs and the oceanographic and biogeochemical conditions of the water column was evaluated using Distance-based Linear Modeling (DistLM) in PRIMER ver.7 (Clarke and Gorley, 2015). Consequently, this analysis provided a pseudo-pseudo-F statistic” and a probability value for each environmental variable after 10000 permutations. To define the environmental variables to be included in the final model, Spearman’s rank correlation was applied to remove highly correlated variables (Spearman’s R > 0.80, p < 0.05) such as POC and PON.

Finally, generalized additive models (GAMs) were constructed to evaluate the effect of oceanographic variability on the diversity indices of OTUs for the most abundant taxonomic groups per sampling site. An automated stepwise procedure was used in the GAMs through the Akaike information criterion (AIC), providing the percentage Deviance Explained (DE), coefficient of determination (R²), and significance value (p) for each model. The models were built using a Gaussian distribution and an identity link function through the mgcv package in R (Wood, 2017; R Core Team A, 2020).

All samples were subjected to DNA extraction and amplification. A total of 4,860,867 reads were obtained for the 18S gene and 2,665,637 for the COI gene, of which 3,696,340 and 1,619,482 passed the quality control and chimera cleanup, respectively. Unique sequences were grouped into ASVs (18S) and OTUs (COI), resulting in 47,871 ASVs for 18S and 2624 OTUs for COI, although only 1681 ASVs and 1631 OTUs could be assigned to metazoans. Rarefaction curves based on the number of ASVs and OTUs showed saturation in most cases, indicating high sampling coverage. Only two samples did not appear to reach saturation (Station 2: 4000-5000 (18S) and 1000–2000 m (COI)) ( Figure 2 ). The raw data are presented in an Excel file ( Supplementary Dataset S1 ), including the abundance per sample, taxonomic assignment, and representative sequences for each ASV and OTU.

Vertical profiles of the diversity indices ( Figure 3 and Supplementary Table S4 ) for the 18S gene exhibited the lowest values in the upper layer at Station 2. The maximum diversity values were also observed at the same station, but in the deepest layer. The COI gene showed the lowest values in the deep layers at Station 4, whereas the maximum values were observed in the 1000–2000 m layer at Station 2. When comparing the stations between 0 m and 4000 m, S, H’, 1-λ, J’, and PD exhibited the same pattern. The diversity indices of the most abundant groups exhibited strong variability (COI gene) ( Supplementary Table S5 ). The diversity was the highest in the upper layer for Peracarida and Halocyprida, whereas Calanoida and Siphonophorae showed the highest values in the middle, bathypelagic layers (1000–3000 m). The other groups exhibited great variability between stations; however, in most cases, the highest values of the S and PD indices were centered between 0 and 2000 m. The 18S gene indicated that the major zooplankton groups were calanoid copepods (>50% relative abundance) at both sampling stations along the water column ( Figure 4 ). Other dominant groups averaged across the water column and stations included Euphausiacea (>10%), Halocyprida (ca. 7%), Doliolida (ca. 5%), Peracarida (ca. 4.9%), Siphonophorae (ca. 3%), other Eucarida (ca. 2.4%), and Salpida (ca. 2.2%). However, the distribution of these groups varied greatly among the strata and sampling stations. For example, the intermediate layer (2000–3000 m) at Station 4 showed a peak in the relative abundance (~ 48%) of Euphausiacea compared with the other groups.

At the species level, provided by the COI data ( Figure 5 ), the number of species varied across groups, with Calanoida having the highest number of species (64), followed by Halocypridae (21 species), whereas fewer than ten species belonged to the remaining groups. In general, the most dominant species showed ubiquitous distribution based on both depth and station, such as Subeucalanus monachus (Calanoida), Euphausia mucronata (Euphausiacea), Erenna sirena (Siphonophorae), Lanceola sayana (Peracarida), and Bentheogennema corbariae (other Eucarida). However, within the Halocypridae group, the distribution of the predominant species Conchoecetta giesbrechti was restricted to the uppermost layer (0–1000 m) at both stations.

The taxonomic assignment based on COI resulted in 177 OTUs that could be identified at the species level. Comparing the number of assignments with the total number of OTUs, it was concluded that there were many unknown (non-assigned) species (~ 80%) throughout the water column at both stations ( Figure 6A ). This pattern of unknown vs. identified species was homogeneous across the depth strata and stations. However, when considering known versus unknown species, based on the relative abundance of OTUs, the number of known species assignments made up 50% of the total number of reads ( Figure 6B ). The highest percentage of unknown species (>78%) was found in the 2000–3000 m (Station 2) and 1000–2000 m (Station 4) strata.

Clustering analysis of 18S and COI genes for the entire diversity of different samples indicated that the two stations did not differ from each other, showing a similar pattern with depth, with five and four clusters for each gene, respectively (SIMPROF test; p < 0.05) ( Supplementary Figure S2 ). ANOSIM of the clusters generated for both genes indicated significant segregation (18S: R = 0.98, p < 0.01; COI: R = 0.93, p < 0.01) ( Supplementary Table S6 ). Temperature, oxygen concentration, and salinity were the main environmental factors that correlated with the PCoA ordination ( Figure 7 ). A similar procedure was performed on the main taxonomic groups identified for both genes, showing two–seven clusters, although some groups indicated an absence of structuring depending on the gene analyzed ( Supplementary Table S6 and Figures S3–S10 ). The Calanoida, Halocyprida, and Peracarida groups showed significant community structuring, concordant with both the genes (R > 0.70; p ≤ 0.01). The structure reflected mainly separate clusters for the layers 0–1000 m and 1000–2000 m, and for the rest of the bathypelagic samples (>2000 m). For the other groups, only one gene indicated structuring, and the patterns differed among them.

A strong gradient of oceanographic conditions was observed at the two study stations, particularly within the upper 500 m layer ( Figure 1B ). The temperature gradually decreased with depth, with an abrupt thermocline at approximately 50 m. Salinity followed a different pattern, exhibiting two abrupt decreases, one in the first 50 m and the other between 500 and 1000 m, associated with Subantarctic Water (SAAW) and Antarctic Intermediate Water (AAIW) respectively (see Supplementary Figure S11 ). The oxygen concentration decreased abruptly below 25 m depth, reaching hypoxic conditions (<50 µM) from 60 to 100 m and anoxic conditions from ~100 to 400 m due to the presence of Equatorial Subsurface Water (ESSW). The two sampling sites showed departures in the hydrographic profiles between stations at 500–1000 m and 2500–3500 m, where lower salinity and more oxygenated conditions were found for Station 4, associated with AAIW and Pacific Deep Water (PDW) mixtures, respectively. Regarding the sources of C and N, the quality and quantity of potential food sources, indicated by POC, PON, C/N ratios, δ¹³C, and δ¹⁵N, were highly variable between the stations and vertical layers ( Supplementary Figure S12 ), revealing highly heterogeneous food conditions in the upper 1000 m compared to the layers below, the latter having more stable food availability in terms of quality and quantity, although much lower than in the upper 1000 m.

Temperature, salinity, oxygen concentration, pressure, and biogeochemical (δ¹⁵N, δ¹³C, C/N, POC, and PON) variables showed high and significant correlations in some cases (Spearman’s R > 0.80, p < 0.05); therefore, POC and PON were analyzed separately. Multivariate analysis between environmental variables and community structure with DistLM (Distance-based Linear Modeling) showed that the patterns differed depending on the gene and group analyzed (see Supplementary Table S7 ). When considering the entire community, temperature and oxygen were significantly correlated with the community structure (p < 0.05), explaining more than 23% of the total deviation in both genes. In some groups, e.g., Euphausiacea (18S) and other Eucarida (COI), C/N explained up to 71% of the community variability (p < 0.05). In contrast, Halocyprida and Peracarida showed a significant correlation with salinity (>25%; p < 0.05) for both genes, whereas Calanoida only presented a significant correlation with the 18S gene, exhibiting the best fit with temperature (deviance > 34%; p = 0.02). The biogeochemical variables POC and PON also had a significant influence on community structure, but with a lower contribution to the total variance than the uncorrelated variables (16–35%, p < 0.05), except for Euphausiacea (up to 75%) in both genes ( Supplementary Table S8 ). Furthermore, when searching for nonlinear predictive models (GAMs) for species diversity, oxygen concentration and temperature were once again the variables with the best fit in the models ( Supplementary Table S9 ). The entire diversity evaluated using the S and PD indices showed that temperature and oxygen concentration significantly explained up to 71.8% of their variance (p < 0.05). The main pattern was similar among the groups, although some other variables, such as salinity and pressure, showed significant relationships of >48% (p < 0.05). Meanwhile, POC and PON only influenced the diversity of specific groups, such as Euphausiacea and other Eucarida, explaining up to 52% of the variability ( Supplementary Table S10 ).

The high zooplankton species diversity found at mid-depths in the subtropical region of the SEP, comparable to or even higher than that found in the upper 1000 m, agrees with the pattern of maximum biodiversity at intermediate depths, as described in the benthic system (Rex, 1981). However, a large part of it was either unknown (>50%) or could not be assigned to any known species found in current genetic diversity databases. This zooplankton diversity is highly structured along the vertical axis of the water column, with a depth pattern parallel to the well-known vertical pattern that describes the benthic system (Rogers, 2015). Concerning the potential drivers controlling zooplankton diversity and its vertical pattern, water temperature and dissolved oxygen appeared to be the most important factors. The role of oxygen in shaping the zooplankton community is not surprising, considering the hydrography of the region, with an acute oxygen minimum zone in the upper 1000 m. Additionally, carbon and nitrogen sources appeared to play a role, although their influence was weaker than those of temperature and oxygen.

Variability in zooplankton diversity with depth has been documented in the tropical and subtropical oceans. Maximum diversity is observed at mid-depths for some specific taxonomic groups, although there is a historical limitation in the sampled depths (max. 3000 m) (Angel, 2003; Yamaguchi et al., 2015; Sommer et al., 2017).This pattern cannot be generalized for the whole community because the maximum diversity for each group may be observed at different depths (Kosobokova and Hopcroft, 2010; Kosobokova et al., 2011). This might be due to environmental preferences, local adaptations, and even lineage evolution (if anything, over evolutionary scales), which may influence the shape of the diversity profiles for each group (Somero, 1992; Seibel and Drazen, 2007; Van der Spoel and Heyman, 2013). For example, Calanoida and Siphonophorae show higher diversity in the intermediate layers (1000–2000 m) (Yamaguchi et al., 2015; Hirai et al., 2021; Ramírez-Flandes et al., 2022), which is possibly associated with intense predatory–driven coevolution (Mackie et al., 1988; Robison, 2004). In contrast, other Eucarida, Euphausiacea, and Halocyprida exhibit higher diversity in the upper layer, possibly due to better adaptation to primary productivity in the photic zone and a close relationship with the shallow oxygen minimum zone (OMZ) typical of this region (Antezana, 2002; Mujica et al., 2022). A few studies have assessed the vertical trend of species richness of the entire zooplankton community, using metabarcoding analysis, and revealed either a peak of diversity in the upper mesopelagic depth (200-300 m) of the North Pacific central gyre (Sommer et al., 2017), or a bimodal vertical pattern with a first peak in the upper 200 m and a second one below 1000 m found in the Indian Ocean (Cheng et al., 2022), and two other regions with a marked oxygen minimum zone in the mesopelagic depths. However, both studies reached maximum depths of 1500 m and 3000 m, respectively. In our study region, analysis of environmental DNA targeting the whole metazoan community suggested a gradual increase in diversity with depth down to 8000 m (Ramírez-Flandes et al., 2022), although this pattern was uncertain for the zooplankton group.

The relative abundance of species observed in this study were consistent with previous research conducted using morphological and genetic approaches, extending down to the abyssopelagic zone (Kosobokova and Hirche, 2000; Vereshchaka et al., 2017; de Puelles et al., 2019; Ramírez-Flandes et al., 2022). In the study region off northern Chile, Calanoida and Euphausiacea dominated the zooplankton communities (Heinrich, 1973; Escribano and Hidalgo, 2000; Antezana, 2002), with Subeucalanus monachus and Euphausia mucronata adapting to hypoxic conditions in the upper 500 m (Binet, 1983; Escribano et al., 2000; Antezana, 2002; Valdés et al., 2007). Halocyprida and Doliolida, which have a detritivorous or sometimes filter-feeding style, might be associated with the high productivity reported in the region that reaches the bathypelagic environment (Nigro et al., 2016; Frischer et al., 2021), playing an important role in vertical carbon transfer to deep waters globally and off the Chilean coast (Martens, 1981; Stone and Steinberg, 2016; Sutherland and Thompson, 2022). In other regions, with lower levels of productivity, such as the Indian Ocean, the zooplankton composition over the vertical axis down to 3000 m was dominated by copepods, hydrozoa, and malacostraca (Cheng et al., 2022), while in the subtropical North Pacific gyre copepods, ostracods and mollusks prevailed in the upper 1500 m (Sommer et al., 2017).

The accumulation curves obtained using the two molecular markers indicated that study sampling may have captured nearly the entire community diversity. However, the percentage of ASVs/OTUs (<62.16%) assigned to metazoans was lower than that reported in other regions using environmental DNA in the abyssopelagic zone (e.g., 74% in Laroche et al., 2020) and zooplankton nets in surface layers (e.g., 69.5% in Schroeder et al., 2021). This suggests that, although the sampling method is appropriate, there is a discrepancy in the assigned percentage due to the molecular markers used and the reference database employed to assign metazoans. Our study showed that less than half of the deep zooplankton diversity could be assigned at the species level. In coastal zones, investigations conducted in the epipelagic layer have assigned 45–80% of the OTUs obtained at the species level (Sommer et al., 2017; Schroeder et al., 2020; Singh et al., 2021). Therefore, there are no complete reference databases for coastal zones, and this situation becomes critical when examining the deep ocean (Chaudhary et al., 2016). Furthermore, identification at the species level assumes a percentage of similarity within species (e.g., 97%) that can vary between lineages, resulting in different degrees of taxonomic classification, particularly in zones where most the species are unknown (Lindeque et al., 2013; Brown et al., 2015; Hirai and Tsuda, 2015). Moreover, the absence of reference databases and different rates of evolution in deep-sea communities may reduce the taxonomic resolution in diversity studies (Pinheiro et al., 2019). This problem can be addressed in the future by generating a reference database of the main deep-sea species, in collaboration with expert taxonomists. This would allow for a faster understanding of the diversity of this environment (Wheeler, 2018; Pinheiro et al., 2019).

Recent diversity analyses have employed a combination of 18S and COI genes as biomarkers to characterize the zooplankton community structure (Cowart et al., 2015; Djurhuus et al., 2018; Stefanni et al., 2018; Zhang et al., 2018; Carroll et al., 2019). However, the results of these studies indicated a discrepancy in the level of taxonomic classification among genes, yielding up to three times higher diversity in the COI gene, coinciding with the findings of other studies (Tang et al., 2012; Clarke et al., 2017). Furthermore, the 18S V9 hypervariable region has a low resolution capability at the species level, providing information on higher taxonomic levels, such as genera or families (De Vargas et al., 2015; Djurhuus et al., 2018; Bucklin et al., 2019; Sawaya et al., 2019; Blanco-Bercial, 2020). In contrast, the COI gene has a greater taxonomic resolution for species identification, but a lower amplification success rate in taxonomically diverse groups (Hebert et al., 2003; Leray et al., 2013; Deagle et al., 2014). This indicates that there is currently no single marker capable of fully resolving the diversity in communities. However, the development of more comprehensive reference databases, standardization, and application of multiple markers can contribute to a more powerful approach, thus allowing a more general view of the deep ecosystem to be conserved in the face of environmental changes (Cristescu, 2014; Deagle et al., 2014; Bucklin et al., 2016; Clarke et al., 2017).

Several studies investigating the epipelagic zone have indicated that the diversity and distribution of zooplankton can be driven by the vertical gradients of temperature and oxygen (Beaugrand et al., 2002; Richardson and Schoeman, 2004; Peterson et al., 2006; Aronés et al., 2009; Gewin, 2010; Tittensor et al., 2010; Seibel, 2011). However, in the deep sea, hydrostatic pressure and sources of organic carbon and nitrogen may be the main factors in structuring communities in the deep ocean (Danovaro et al., 2003; Jamieson et al., 2010; Tamburini et al., 2013; Ichino et al., 2015). Our data showed that, despite the significant correlations between diversity and biogeochemical variables (POC, PON, and C/N), temperature and oxygen were the main correlated variables. In the same context, it is worth noting that the highest species diversity was found in the 1000 and 2000 m layers, far below the influence of the OMZ, which is located between 200 and 500 m in this area (Ulloa et al., 2012). The OMZ is known to act as an ecological barrier to species dispersion and migration from the photic zone to the layers below (Donoso and Escribano, 2014; Wishner et al., 2020). It may then be suggested that such a barrier promotes a higher diversity than that reported for the upper layers in this region. It is also important to note that the two stations were different when comparing the peaks of diversity. In this regard, the hydrographic characteristics of the Southeast Pacific indicate that temperature and oxygen are not linearly linked to depth; rather, they differ between sampling stations owing to the spatial variability in water mass properties, as reflected by the dissolved oxygen concentration between 2500 and 3500 m. This allowed for disentanglement of the influence of depth from the environmental variables in this study. In addition, optimal ranges of environmental conditions for some species cause community distribution to vary across the water column (Laakmann et al., 2009; Kosobokova and Hopcroft, 2010), and they are not always tightly linked to an increase in depth. The case of this study robustly demonstrates this conclusion using both the entire zooplankton community and specific groups.

When interpreting the mechanisms linking diversity and environmental drivers, it must be recognized that temperature is an environmental parameter that influences the rates of development, respiration, hatching, mortality, and zooplankton distribution (Nogueira et al., 2018; Jo et al., 2019; Biard and Ohman, 2020). Oxygen concentration is a key factor influencing the growth, metabolic processes, and behavioral processes of various pelagic organisms (Ekau et al., 2010). The findings of this study suggest that these variables might modulate the community, at least in the upper 2000 m. In the deep ocean of the SEP, variable water masses, including the ESSW, AAIW, and PDW, generate a gradient across the water column (50–2000 m), showing an intense OMZ in the surface and subsurface layers and a stronger temperature gradient than that observed in the deeper layers (Morales et al., 1996; Escribano et al., 2004; Ulloa et al., 2012). Strikingly, the weak gradients in deep waters (>2000 m) might still contribute to a higher rate of allopatric specialization by generating suboptimal conditions for dispersal, thus promoting high diversity along the Chilean coast, as suggested for other regions (Wilson and Hessler, 1987; Grassle, 1989; Ekau et al., 2010). These factors acting across the water column can have different effects on the benthic system, where the characteristics of the ocean bottom (e.g., muddy, rocky, sandy) and processes, such as burial and sediment transport, play critical roles in influencing life diversification and species establishment. Nevertheless, similar vertical patterns of diversity between benthic and pelagic ecosystems suggest the existence of a common ecological process that controls species diversification in the ocean. The vertical diversity pattern observed in our study was derived from data collected from only two sampling stations, which could be considered a relatively small sample size from a statistical perspective. However, it is important to note that these sampling stations represented patterns observed in large water masses, with a substantial volume of water (> 20,000 m³) sampled by MOCNESS in each stratum. Despite these limitations, the findings of this study are unique because comparative studies conducted in other regions at similar depths are lacking.

The significant correlations between zooplankton diversity and environmental variables should be interpreted with caution, this because some sources of bias may have affected the study results. For instance, the vertical partitioning of the water column every 1000 m was based on an arbitrary design and did not necessarily reflect the actual ecological structure of zooplankton distribution across the vertical axis. However, some zooplankton species perform daily vertical migrations mostly across the upper water column, which can cause daytime-nighttime species exchange and mixing between the upper hundred meters where the environmental gradient is steepest (Sommer et al., 2017; Tutasi and Escribano, 2020). These day-night effects due to vertical migration would be, however, not significant for the comparison between layers with 1000 m spacing. An additional source of bias may have arisen from the characteristics of the sampling gear, such as the mesh size of the nets, towing speed, and opening size, which can affect the quality and quantity of samples (Skjoldal et al., 2013). For instance, larger, or smaller animals could have been under-sampled, and even some dead animals or parts of them could have been caught in the nets. However, it is reasonable to presume that the study samples adequately represented the living mesozooplankton community (0.2-2.0 mm) at each strata. Therefore, in future research pertaining to the midnight zones, it will be of utmost importance to consider these factors and compare the observed vertical patterns with those of various regions worldwide. Therefore, in future research pertaining to the midnight zones, it will be of utmost importance to consider these factors and compare the observed vertical patterns with those of various regions across the globe.