Metabarcoding analysis of marine zooplankton confirms the ecological role of a sheltered bight along an exposed continental shelf

Zooplankton plays an essential role in marine ecosystems as the link between primary producers (phytoplankton) and higher trophic levels in food webs, and as a dynamic pool of recruits for invertebrates and fish. Zooplankton communities are diverse with a patchy distribution at different spatial scales, influenced by oceanographic processes. The continental shelf of eastern South Africa is narrow and exposed to the western‐boundary Agulhas Current, with some shelter against strong directional flow provided by the broader KwaZulu‐Natal Bight, a coastal offset adjacent to an estuary. We compared zooplankton species richness, diversity and relative abundance of key taxa among sheltered and exposed shelf areas using metabarcoding and community analysis, to explore the ecological role of the bight in a highly dynamic ocean region. Metabarcoding recovered higher richness and diversity at a finer resolution than could previously be achieved with traditional microscopy. Of 271 operational taxonomic units (OTUs) recovered through metabarcoding, 63% could be matched with >95% sequence similarity to reference barcodes. OTUs were dominated by malacostracan crustaceans (161 spp.), ray‐finned fishes (45 spp.) and copepods (28 spp.). Species richness, diversity and the relative abundance of key taxa differed between sheltered and exposed shelf areas. Lower species richness in the bight was partly attributed to structurally homogeneous benthic habitats, and an associated reduction of meroplanktonic species originating from local benthic–pelagic exchange. High relative abundance of a ray‐finned fish in the bight, as observed based on fish eggs and read counts, confirmed that the bight is an important fish spawning area. Overall, zooplankton metabarcoding outputs were congruent with findings of previous ecological research using more traditional methods of observation.


| INTRODUC TI ON
Zooplankton plays an essential role in marine ecosystems as the link between primary producers (phytoplankton) and higher trophic levels in food webs, and by forming a dynamic pool of recruits for multiple invertebrate and fish species (Brierley, 2017;Peijnenburg & Goetze, 2013;Richardson, 2008).Zooplankton have high taxonomic diversity and abundance and are distributed unevenly in oceans, with dispersal and aggregation driven by a combination of physical processes and biological interactions at different spatial scales (Folt & Burns, 1999;McManus & Woodson, 2012;Pinel-Alloul, 1995).
The coastal oceanography of southeast Africa is influenced by the Agulhas Current, one of the strongest western-boundary currents in the world (Lutjeharms, 2006).The Agulhas Current closely follows the shelf-edge polewards, bringing tropical waters from the Western Indian Ocean southwards along the coast into subtropical and temperate regions off eastern South Africa (Figure 1).The continental shelf of eastern South Africa is narrow and slopes down steeply after reaching 100 m depth, some 3-11 km from the coast, except at the broader KwaZulu-Natal (KZN) Bight, a coastal offset 160 km long with a maximum shelf-width of 45 km (Schumann, 1988).The KZN Bight is the only sizeable shelf area along a coastline spanning about 1000 km.Along its seaward edge, the KZN Bight diverts the Agulhas Current offshore and diminishes its velocity gradient, causing shelfedge upwelling cells, cyclonic lee-trapped eddies and counter-currents (Roberts et al., 2016).Some of the Agulhas Current boundary features are propagated southwards (downstream) of the bight as transient meanders or eddies along the coast, collectively known as the Natal Pulse (Lutjeharms, 2006;Lutjeharms & Roberts, 1988;Roberts et al., 2016).
Shelf habitats and biota differ markedly between the sheltered KZN Bight and the exposed narrow shelf to its north and south.
Whereas the narrow shelf waters are mostly oligotrophic, the inner bight receives seasonal sediment, organic matter and nutrient input from riverine discharge (De Lecea & Cooper, 2016;Scharler et al., 2016).Zooplankton biomass over the bight increases during winter and close inshore (Pretorius et al., 2016) when fish and invertebrate larvae become seasonally abundant in the zooplankton (Collocott, 2016;Hutchings et al., 2002).Particle dispersal models (Singh et al., 2018(Singh et al., , 2019) ) have demonstrated that passively drifting zooplankton will either be retained in the bight by submesoscale processes, disperse downstream along the coast, or become entrained in the Agulhas Current and lost offshore.
DNA barcodes have been developed mainly for fishes (Steinke et al., 2016) and decapod crustaceans (Govender et al., 2019(Govender et al., , 2022) ) in the region, several of which are important to fisheries.Despite their high diversity and abundance in the zooplankton, barcode records remain sparse for most holoplanktonic taxa (Singh et al., 2021).
We used metabarcoding of zooplankton tow net samples collected in surface waters of eastern South Africa combined with a community analysis approach to test the hypothesis that distinct zooplankton communities occur in sheltered (KZN Bight) and exposed narrow shelf areas.The hypothesis was based on the findings of previous oceanographic and ecological studies of the region, which described semipermanent eddies and counter-currents within the bight (Roberts et al., 2016;Schumann, 1988), that it is a spawning area for fishes (Hutchings et al., 2002) and that biological productivity in the bight is enhanced by nutrients originating from shelf-edge upwellings and land-derived inputs (Fennessy et al., 2016).The indicators used to test the hypothesis were zooplankton species richness and diversity, and relative abundance of key taxa at exposed and sheltered areas.
The congruence of metabarcoding outputs with findings from previous ecological research relating to retention of zooplankton over the bight and its role as a spawning area was explored.

| Study area
Four cross-shelf transects were sampled for zooplankton along the east coast of South Africa, at Isimangaliso, Thukela, Durban and Aliwal Shoal (hereafter called Aliwal; Figure 1).The transects represented narrow shelf and sheltered bight locations, with different levels of exposure to the inshore edge of the Agulhas Current (Table 1).Pelagic habitats ranged from warmer oligotrophic waters at Isimangaliso and along the offshore edges of transects to bight waters enriched with land-derived minerals and nutrients deposited by rivers in flood.Benthic habitats adjacent to transects (a potential source of meroplankton) differed between exposed shelf and sheltered bight transects (Table 1).

| Sampling strategy
Each of the four cross-shelf transects was sampled for zooplankton at four discrete sampling stations, located at the 20-, 50-, 100and 200-m depth isobaths (Figure 1).Sampling was conducted at night, when zooplankton typically rise to the surface to feed (Lampert, 1989).Sampling gear consisted of a plankton ring net (500μm mesh; 0.8-m ring diameter) which was towed <5 m below the sea surface, for 5 min at a time, at a boat speed of 2-3 knots.Three replicate tows were made at each sampling station, so that the variability inherent in the method and sampling gear could be quantified at the scale of sampling stations.Samples were washed from the cod-end into a jar with 95% ethanol and stored at −20°C.
Adverse weather at the Isimangaliso transect restricted sampling to a single tow per sampling station and none at the 200-m sampling station.All field sampling took place between August and November 2018.

| Extraction and quantification of genomic DNA
Individual tow-net samples (n = 37) were homogenized for 45 s using a consumer blender (Defy PB7354X, 350 W and 22,000 rpm).The blender blades and container were sterilized with 4% industrial bleach and washed with 95% ethanol between applications.Three subsamples per homogenate were taken to increase sequencing depth and improve diversity estimates (Lanzén et al., 2017).Each subsample (10 ml of zooplankton) was centrifuged at 1200 x g for 1 min, repeated to remove excess ethanol; thereafter 40 mg of tissue was transferred to a sterile tube.DNA was extracted using the Qiagen DNeasy Blood and Tissue Kit by adding 180 μl ATL buffer and 40 μl proteinase K to the tissue for overnight lysis at 56°C following the manufacturer's standard protocol.
The extracted DNA from the three subsamples per homogenate were

| PCR amplification, library preparation and high-throughput DNA sequencing
Polymerase chain reactions (PCRs) were performed in triplicate to minimize bias and amplification errors (Dopheide et al., 2019).

Illumina sequencing was performed at the KZN Research and
Innovation Platform (KRISP).Libraries were purified using 1.8× AmpureXP purification beads (Beckman Coulter).Index PCR was performed using the Nextera XT Index Kit (Illumina).Libraries were further purified using 0.6× AmpureXP purification beads (Beckman TA B L E 1 Four cross-shelf transects sampled for zooplankton in surface waters showing the general bottom topography, exposure to the inshore edge of the Agulhas Current (AC), and pelagic and benthic habitats

| Assignment of operational taxonomic units (OTUs)
The dada2 algorithm (Callahan et al., 2016)  Rarefaction curves of sequences obtained from qiime2 were generated using vegan version 2.5.6 (Oksanen et al., 2018) in R to determine whether samples were sequenced to a sufficient depth.
Sequences generated using qiime2 were searched against BOLD systems version 4 (https://www.boldsystems.org)and NCBI GenBank (https://www.ncbi.nlm.nih.gov/genbank/).Matches above 95% sequence similarity were assigned to OTUs.We selected 95% to maximize genetic diversity recovered while controlling for the effects of sequence errors (Porter & Hajibabaei, 2020).Below the 95% threshold a multiple sequence alignment was performed with the default settings in mafft version 7.470 (Katoh et al., 2019)
Species richness and relative abundance were compared between individual sampling stations on each transect (alpha-diversity) and between the four transects (beta-diversity) using phyloseq and plotted with ggplot2 (Wickham, 2009).We used ACE and Chao1 for species richness and Shannon, Simpson, and Fisher's alpha to determine alpha-diversity.Nonmetric multidimensional scaling (NMDS) analyses based on the Bray-Curtis dissimilarity index was performed in phyloseq to determine beta-diversity.The Bray-Curtis  Permutational analyses of variance (PERMANOVAs) using the "adonis" function in vegan was used to assess whether benthic habitat affected zooplankton community composition.Clustering of sampling stations was investigated by using the unweighted pair group average (UPGMA) in phyloseq.

| High-throughput sequencing results
Sequencing was efficient with minimal filtering needed when merging the paired-end reads for all 37 zooplankton libraries.A total of 2.6 million read counts were consolidated into 102,875 merged reads and 51,439 paired reads (Table 2).Consolidation of 1627 sequences resulted in 271 OTUs of which 170 (63%) could be matched with >95% sequence similarity to sequences on BOLD or GenBank.
Rarefaction curves based on sequences were asymptotic for all libraries (Figure 2), indicating that the data approached saturation and that sequencing depth was sufficient.

| Relative abundance of key taxa
Actinopterygii accounted for most reads overall (41% of all reads) indicating that they had the highest relative abundance in all samples combined (Figure 4a) and dominated at Thukela where they made up 70% of reads (Figure 4b).Malacostraca were dominant at Isimangaliso (53% of reads) followed by Actinopterygii (25% of reads).At Durban, Actinopterygii (57% of reads) and Malacostraca

| Community analysis
Both the NMDS ordination plot (Figure 6a) and UPGMA cluster analysis (Figure 6b) indicated an overlap between the Thukela and Durban zooplankton communities but found those at Isimangaliso and Aliwal to be discrete.The locational PERMANOVA test indicated that the community composition differed significantly between the four transects (p < .001,R 2 = .239,F = 7.408), but the low R 2 value suggested that the model could only account for a small part of the observed variation.PERMANOVAs on physical variables collected during sampling (sea-surface temperature, salinity, pH and dissolved oxygen levels) likewise failed to account for substantial variation in zooplankton distribution profiles (data not shown).Benthic habitats differed between the four transects (Table 1)-and although superimposed on the NMDS plot (Figure 6a), they were not modelled in conjunction with locations because of collinearity.

| DISCUSS ION
Metabarcoding has been used in recent ecological research to determine the effects of latitude, season, water temperature, biogeographical boundaries and ocean currents on zooplankton community structure and distribution patterns (Blanco-Bercial, 2020; Casas et al., 2017;Ershova et al., 2019;Hirai et al., 2020;Macher et al., 2020 ;Pitz et al., 2020).The present study is the first to apply it to ecological research in the Western Indian Ocean, a tropical/subtropical province of the Western Indo-Pacific region with exceptionally high marine biodiversity and many endemic species (Spalding et al., 2007).We used metabarcoding to test the hypothesis that zooplankton communities in the sheltered KZN Bight differ from those over a steep narrow shelf exposed to the Agulhas Current.
The hypothesis was based on findings from a range of previous studies, showing the occurrence of submesoscale eddies and countercurrents over the bight (Roberts et al., 2016;Schumann, 1988); temporal retention of passively drifting particles over the bight, using satellite-tracked drifters and particle dispersal models (Guastella & Roberts, 2016;Singh et al., 2018); seasonally greater zooplankton biomass over the bight than other areas of the greater Agulhas Current system (Pretorius et al., 2016); and use of the bight by fish and invertebrates for spawning (Hutchings et al., 2002).The congruence of metabarcoding outputs with findings of previous ecological research from the region (summarized by Fennessy et al., 2016) was evaluated.
We estimated species richness and relative abundance, with the latter a key component in diversity estimates and ecosystem health assessments (Nichols et al., 2018).Lower confidence in relative abundance estimates (measured as read counts in metabarcoding) has been expressed by Elbrecht and Leese (2015), Piñol et al. (2015) and Jusino et al. (2019), although these estimates have been successfully used for nematode communities (Schenk et al., 2019), arthropod taxa in songbird diets (Verkuil et al., 2020) and zooplankton research (Macher et al., 2020;Pitz et al., 2020).In an experimental approach, Govender et al. (2022) found a strong relationship between known zooplankton species ratios (including abundant and rare species) and read counts derived from metabarcoding using COI minibarcode primer sets.Based on the above, and using the same minibarcode primers, we considered proportional read counts as good indicators of relative abundance of key taxa in the present study.
Metabarcoding recovered a much higher species richness at finer resolution (271 OTUs of which 170 matched reference sequences at >95%) than could previously be achieved with traditional microscopy (Singh et al., 2021;Turon et al., 2020).The derived community structure was dominated by Malacostraca, a faunal class specifically targeted with COI minibarcode primers designed for marine lobsters (Govender et al., 2019), prawns (Dendrobranchiata), shrimps (Caridea) and crabs (Brachyura) (Govender et al., 2022).Minibarcode primers amplify smaller gene fragments and can increase species detection rates when DNA is degraded, as expected from bulk tow net-and environmental DNA samples.The use of universal barcode primers for fish (Leray et al., 2013;Ward et al., 2005) combined with good representation of fish species on online barcode reference databases (Steinke et al., 2016) resulted in high richness estimates of ray-finned fish (Actinopterygii).The species richness estimates for Copepoda, Gastropoda, Hydrozoa, Sagittoidea, Ostracoda, Thaliacea and Branchiopoda were lower at all sampling stations and potentially underestimated, because these taxa were not targeted by the taxon-specific primer cocktails used.Use of multigene markers (Questel et al., 2021), taxon-specific primers for all groups (Govender et al., 2022), and addition of regional and endemic species records to Western Indian Ocean barcode reference databases will undoubtedly increase the numbers of species that can be identified in the region (Singh et al., 2021).
Distinct zooplankton communities were found at the four transects, with overlap between the sheltered bight communities at Thukela and Durban.The two exposed shelf communities at Isimangaliso and Aliwal (north and south of the bight, respectively) differed from the two sheltered communities and from each other.
Species richness and diversity were comparatively lower in sheltered than exposed communities.Abundance was dominated by Actinopterygii at both the sheltered bight communities, but they were relatively less abundant at exposed communities, which were dominated by Malacostraca.Copepod abundance increased markedly at Aliwal, where Actinopterygii abundance was lower than elsewhere.
Higher zooplankton species richness at exposed than sheltered transects could partly be explained by benthic habitat heterogeneity and benthic-pelagic exchange of meroplankton (Boero et al., 1996;Ershova et al., 2019;Lindley et al., 1995;Marcus & Boero, 1998).The underlying assumptions were that species richness in benthic and pelagic environments would be correlated; and that heterogenous and structurally complex benthic habitats (such as the reefs at Isimangaliso and Aliwal) would have higher species richness than structurally homogenous habitats with fewer niches (such as the muddy substrates with few rocky outcrops at Thukela; Zeppilli et al., 2016).The high richness and diversity of zooplankton at Aliwal was attributed to a heterogenous benthic environment, dominated by a mix of high-profile coral reefs (Olbers et al., 2009;Schleyer, 2000), sandy seafloor and even muddy habitats at the mouth of the nearby Umkomaas River (see Table 1).The richness and diversity of zooplankton were also high at Isimangaliso, where the benthic environment was also heterogenous, including reefs with hard and soft corals.Even so, zooplankton richness and diversity at Isimangaliso were potentially underestimated compared to the other three transects, because fewer samples were available.
Greater similarity among zooplankton communities at the Thukela and Durban transects within the KZN Bight, compared to exposed narrow-shelf Isimangaliso and Aliwal transects, suggests a higher level of connectivity or exchange of organisms within the bight.Circulation of shelf waters in the southern bight is characterized by the semipermanent cyclonic Durban Eddy which occasionally "swirls" northwards to the Thukela River mouth in the central bight (Guastella & Roberts, 2016;Roberts et al., 2016).The resulting semiclosed circulation system over the bight is expected to facilitate northward advection and retention of zooplankton in inshore counter-currents, thus explaining similarities in zooplankton composition between the two transects, relative to the exposed transects with mainly directional currents.Acoustic Doppler Current Profile (ACDP) data and satellite-tracked drifters confirmed that exchange of drifting particles would occur at a higher intensity between sheltered transects in the bight than between sheltered (Durban) and exposed (Aliwal) transects, the latter via a southern extension of the Durban Eddy (Guastella & Roberts, 2016).
Highly abundant fish eggs observed in tow-net samples collected at the Thukela and Durban transects, combined with metabarcoding of the samples, could confirm the findings of previous studies (Collocott, 2016;Hutchings et al., 2002), identifying the KZN Bight as an important spawning area for fish.Read counts in the present study indicated a high abundance of Actinopterygii at the two sheltered transects, of which most were Scomber japonicus, a common small pelagic fish from the region.S. japonicus is a broadcast spawner that produces large numbers of small eggs and larvae that form dense patches in the zooplankton during winter and spring (August to November; Beckley & Leis, 2000;Connell, 2001;Hutchings et al., 2002).Tow-net samples for the present study were collected between August and November 2018, thus overlapping the known spawning season of S. japonicus.The metabarcoding results were therefore congruent with previous information identifying the KZN Bight as an important spawning area for fish (Beckley, 1993;Hutchings et al., 2002;Lamberth et al., 2009).
In conclusion, metabarcoding of zooplankton communities recovered much higher species richness and diversity (measured as OTUs) than could previously be achieved with microscopy.Zooplankton species richness and relative abundance of key taxa (measured as read counts) differed between sheltered bight and exposed narrowshelf transects.Existing theory of zooplankton retention within the bight by eddies and counter-currents, and increased connectivity was supported by community analysis.Lower species richness in the bight was attributed to structurally homogenous benthic habitats, and an associated reduction in local meroplanktonic species originating from benthic-pelagic exchange.Metabarcoding identified abundant fish eggs collected in the bight as S. japonicus, supporting existing theory that the bight is an important fish spawning area.
Overall, zooplankton metabarcoding outputs were congruent with the findings of previous ecological research using more traditional methods of observation.

F
Location of cross-shelf transects at Isimangaliso, Thukela, Durban and Aliwal, with sampling stations at 20-, 50-, 100-and 200-m depth isobaths.The proximity of the Agulhas Current to the coast and eddies/counter-currents over the sheltered KwaZulu-Natal Bight are shown then pooled and stored at −20°C.Laboratory blanks were used consistently during DNA extractions to test for contamination.
and optimized manually to ensure homology.Genetic distance and neighbor-joining trees (NJ) were used to derive clusters which were assigned to OTUs at the class level.Clusters assigned to the same species were merged by summing up reads using phyloseq (McMurdie & Holmes, 2013) in R version 4.0.2(R Development Core Team) to prevent analyses of intraspecific variation.

TA B L E 2
High-throughput summary statistics across all samples collected at cross-shelf transects at Isimangaliso, Thukela, Durban and Aliwal F I G U R E 2 Rarefaction curves of samples collected at individual sampling sites along the four cross-shelf transects at Isimangaliso, Thukela, Durban and Aliwal was based on transformed data using the TSS normalization protocol in phyloseq.

(
38%) were abundant, and at Aliwal, Copepoda were abundant (47% of reads).A comparison of relative abundance of key taxa at F I G U R E 3 (a) A Venn diagram with species composition plots of shared species clusters and (b) beta-diversity indices comparing the difference between Isimangaliso, Thukela, Durban and Aliwal transects Species composition bar graphs for (a) all data combined (overall), (b) by cross-shelf transect (Isimangaliso, Thukela, Durban and Aliwal) and (c) by sampling station stations along each transect (Figure 4c) indicated high variability across the shelf at Thukela, where Actinopterygii were highly abundant at the 20-m (86% of reads) and 100-m (90% of reads) stations, but Malacostraca were more abundant at the 50-m station (57% of reads).At Durban, Actinopterygii dominated at the 50-m (80% of reads) and 200-m stations (76% of reads), but Malacostraca were more abundant at the 100-m station (70% of reads).These results indicate spatial patchiness at a cross-shelf scale.F I G U R E 5 Boxplots of observed species richness, Shannon's diversity index, Simpson's diversity index and Fisher's alpha across the four transects Beta-diversity analysis to estimate the dissimilarity and similarity between Isimangaliso, Thukela, Durban and Aliwal and their respective sampling stations.(a) NMDS ordination plot based on Bray-Curtis dissimilarity, stress value = 0.09939, with the underlying benthic habitats most common at each transect superimposed, and (b) UPGMA cluster analysis based on Bray-Curtis dissimilarity