INTRODUCTION
Zooplankton form some of the most abundant and diverse biological communities in aquatic environments and constitute an integral component of both marine and freshwater ecosystems (Bucklin, Lindeque, Rodriguez-Ezpeleta, Albaina, & Lehtiniemi, 2016; Djurhuus et al., 2018; Xiong et al., 2019). Estimating the richness of zooplankton species remains a priority in several fields of research, from ecology to conservation biology (Carugati, Corinaldesi, Dell’Anno, & Danovaro, 2015; Gaston, 2009). Such studies are particularly important in marine ecosystems where ecological processes are often maintained by both complex biotic interactions and environmental drivers (Palumbi et al., 2009; Steinberg & Landry, 2017). For example, zooplankton communities play a central role in regulating biogeochemical cycles by transferring carbon from primary producers to higher trophic levels (Bucklin et al., 2019; Everaert, Deschutter, De Troch, Janssen, & De Schamphelaere, 2018; Lindeque, Parry, Harmer, Somerfield, & Atkinson, 2013; Steinberg & Landry, 2017). Moreover, due to their small size, limited capacity for self-dispersal, and high sensitivity to environmental change, zooplankton are also considered to be useful bio-indicators that may be used to evaluate the health of marine ecosystems (Buttay, Miranda, Casas, González-Quirós, & Nogueira, 2015; Johnston, Mayer-Pinto, & Crowe, 2015; Parmar, Rawtani, & Agrawal, 2016; Yang & Zhang, 2019). Long-term studies have reported that both environmental change and ecological processes shape the spatial and temporal abundance of zooplankton and the taxonomic composition of their communities (Buttay et al., 2015). As such, environmental change may lead to the loss of zooplankton biodiversity (Gazonato Neto, Silva, Saggio, & Rocha, 2014; Parmar et al., 2016), which in turn may affect ecosystem services and result in economic consequences (Beaugrand, Edwards, & Legendre, 2010; Bucklin et al., 2016; Everaert et al., 2018; Johnston et al., 2015).
The Gulf of Mexico (GoM) is an example of a system that is subject to changing environmental conditions. For example, the large scale near-surface circulation in the GoM is largely dominated by the energetic Loop Current (LC; Damien et al., 2018). The northward penetration of the LC is often accompanied by the formation and release of anticyclonic LC eddies (LCEs), which travel towards the western boundary of the GoM (Damien et al., 2018; Sheinbaum, Athié, Candela, Ochoa, & Romero-Arteaga, 2016). Together with other smaller cyclonic and anticyclonic eddies that are formed within the gulf, LCEs are considered to constitute the principal source of mesoscale variability within this ecosystem (Damien et al., 2018; Hamilton, 2007; Jouanno et al., 2016). Moreover, LCEs and smaller eddies have been found to constrain the distributions of nutrients, zooplankton, and other planktonic species in the gulf (Biggs & Ressler, 2001; Linacre et al., 2015). Based on the distribution of chlorophyll (Damien et al., 2018; Muller-Karger et al., 2015; Salmerón-García, Zavala-Hidalgo, Mateos-Jasso, & Romero-Centeno, 2011), the GoM may be conceptually divided into two main areas: (1) a central oligotrophic area and (2) the eutrophic inshore waters that semi-surround it (Damien et al., 2018; Uribe-Martínez, Aguirre-Gómez, Zavala-Hidalgo, Ressl, & Cuevas, 2019). Nevertheless, few studies have evaluated plankton distributions in the eutrophic inshore waters of the GoM (Biggs & Ressler, 2001; Damien et al., 2018) nor the distribution patterns of zooplankton and other planktonic species over the entire gulf. This constitutes the principle limitation in determining if the zooplankton distribution within the GoM also follows the pattern that has been established for chlorophyll.
The main challenges in evaluating the spatial and temporal patterns of zooplankton communities lie in the taxonomic complexity of their assemblages, which include a considerable number of morphologically diverse and cryptic species and a lack of diagnostic characteristics for immature and larval developmental stages (Bucklin et al., 2016). In recent years, molecular methods have generated new and powerful approaches for assessing zooplankton biodiversity by overcoming the main limitations associated with traditional taxonomic surveys, such as the amount of time needed to identify specimens and the need for advanced taxonomic expertise (Creer et al., 2016; Cristescu, 2014; Zhang, Chain, Abbott, & Cristescu, 2018). In particular, the metabarcoding approach (Creer et al., 2016) is considered to be one of the most comprehensive means to holistically evaluate zooplankton assemblages, as it combines DNA barcoding and high-throughput sequencing to evaluate the taxonomic composition of a sample with any source of environmental DNA (eDNA; Stefanni et al., 2018; Zhang et al., 2018). However, metabarcoding results strongly depend on having adequate taxonomic coverage and a systematic resolution of the chosen molecular marker (Bucklin et al., 2016; Larke, Beard, Swadling, & Deagle, 2017; Zhang et al., 2018).
Nuclear-encoded ribosomal RNA fragments, especially hypervariable regions of the 18S rRNA gene, were prime targets in early studies because they are able to provide conserved primer binding sites with broad taxonomic coverage across the eukaryotic domain (Larke et al., 2017; Lindeque et al., 2013). Nevertheless, the 18S rRNA gene often lacks the taxonomic resolution of protein-coding genes, such as mitochondrial cytochrome oxidase c subunit I (COI; Larke et al., 2017; Machida, Leray, Ho, & Knowlton, 2017; Zhang et al., 2018). Indeed, COI undergoes faster rates of evolution than that of the 18S rRNA gene; hence, its great genetic variability can facilitate the systematic classification of even closely related taxa (Zhang et al., 2018). Nevertheless, the wobble effect (i.e., meaningless nucleotide change in the position of the third codon) can increase the occurrence of primer mismatches among taxa and consequentially induce taxonomic bias during PCR amplification (Larke et al., 2017; Piñol, Mir, Gomez-Polo, & Agustí, 2015). These gene-related limitations have led many researchers to propose the need to implement a multi-locus approach in metabarcoding surveys, as the synergistic information of different loci may represent a trade-off between increasing the resolution needed to identify metazoan taxa (Carroll et al., 2019) and reducing the occurrence of false negative and/or false positive outcomes (Bucklin et al., 2016; Zhang et al., 2018). Accordingly, it has become increasingly accepted that the combined use of nuclear and mitochondrial genes may facilitate the assessment of zooplankton communities due to the resulting improved sensitivity for detecting cryptic as well as intra-species genetic diversity (Carugati et al., 2015).
Despite the ecological importance of zooplankton and an improved understanding of the dynamics of the GoM, few efforts have been made to describe the composition and distribution of the entire zooplankton community that is dispersed throughout the deep-water region of the gulf in the Exclusive Economic Zone (EEZ) of Mexico. Moreover, the extent to which the physical environment may shape the structure of those communities remains uncertain. In this study, we hypothesized that the structure of the zooplankton community would reflect regional and seasonal environmental features. Accordingly, the main goals of this study were to (i ) assess the spatial and temporal variability of the zooplankton community from the deep water region of the GoM within the Mexican EEZ, (ii ) explore the possibility that structural changes in the zooplankton community might be related to environmental factors, and (iii ) describe the regional and temporal diversity of the community.
In order to address these hypotheses, we characterized the structure and variability of the zooplankton community using a metabarcoding approach based on the genetic information of both the hypervariable V9 and Leray_Folmer 1 regions of the 18S rRNA and COI genes, respectively. We chose this combination of genes because it currently provides the best compromise between the detection and resolution of zooplankton taxa. Specifically, the hypervariable V9 region may be used to amplify across all known and unknown metazoan phyla due to its highly conserved flanking regions (Amaral-Zettler, McCliment, Ducklow, & Huse, 2009). In contrast, we posited that COI could provide greater taxonomic resolution (or α-biodiversity) due to its high variability (Heimeier, Lavery, & Sewell, 2010). For both molecular markers, detected amplicon sequence variants (ASVs) were used to estimate zooplankton alpha and beta diversity across the various sampled areas of the GoM. First, we analyzed the genetic information from three cruises, and we pooled all data from each locus to evaluate the potential zooplankton patterns in a comprehensive analysis. Our efforts were directed towards investigating the effects of temperature, salinity, oxygen, depth, latitude, longitude, water density, and florescence (as a proxy for chlorophyll abundance) in shaping the structure of the zooplankton community in order to better understand the physical-biological interactions controlling spatial, seasonal, and inter-annual changes of zooplankton assemblages.