Taxon-specific mini-barcode primers were designed following Govender et al. (2019). For in silico analysis to determine the most variable COI portion, we downloaded sequences of 954 copepod species (1,910 sequences; 110 families; 345 genera), 65 euphausiid species (130 sequences; 2 families; 11 genera), 25 chaetognath species (focusing on Sagittoidea; 50 sequences; 5 families; 13 genera), and 404 hydrozoan species (808 sequences; 83 families; 200 genera) from GenBank and BOLD (www.ncbi.nlm.nih.gov/genbank; http://www.boldsystems.org/; accessed: 1 December 2023; Suppl. material 1). Two individuals per species from different geographical areas were included, where available, to capture within-species variation.

Sequences were aligned separately for each of the four datasets using Clustal X2.1 (Larkin et al. 2007) and manually optimized to ensure homology using BioEdit 7.2.5 (Hall 1999). All possible mini-barcode fragments per dataset were estimated using sliding window analysis (SWAN; Proutski and Holmes 1998) in the Species Identity and Evolution (SPIDER; Brown et al. 2012) package in R (http://www.r-project.org). The “slideAnalyses” function was used to generate windows varying in size from 100 to 300 bp, which were moved along alignments at 10 bp intervals. Two mini-barcode fragments per window length were chosen for further analysis based on: (1) high mean Kimura 2-parameter (K2P) distance; (2) few zero pairwise non-conspecific distances (minimum interspecific distance for each individual); and (3) a high proportion of clades shared between neighbor-joining (NJ) trees derived from the mini-barcode region compared to reference trees constructed from full-length DNA sequence alignments (Govender et al. 2019).

Based on the SWAN analysis, a total of 42 potential mini-barcode fragments were selected for each dataset. Maximum likelihood (ML) analysis was conducted on potential mini-barcode fragments (comparison trees) and full-length sequence alignments (reference trees) for each dataset using Garli 0.951 (Zwickl 2006). The K2P model (Kimura 1980) of sequence evolution was used for all ML analyses, as it is the most often applied by the DNA barcode community and BOLD. Ktreedist (Soria-Carrasco et al. 2007) was used to compare the mini-barcode ML trees to the reference trees. Ktreedist calculates the minimum branch length distance (K-score) from one phylogenetic tree to another. The K-scores (topology and branch length differences) and Robinson–Foulds symmetric differences (topological differences) were also calculated. For both methods, low values imply high similarity between mini-barcode and reference trees. Mini-barcode regions for each dataset were selected based on the lowest K-scores and Robinson–Foulds symmetric differences.

DNA barcode gap analyses were conducted to confirm that the selected regions would be able to distinguish between species. The K2P nucleotide substitution model in MEGA 6.0 (Tamura et al. 2013) was used to calculate intra- and interspecific genetic distances. The maximum intraspecific distance value was subtracted from the minimum interspecific distance value to determine the barcoding gap (Meier et al. 2006). The Jeffries–Matusita distance (J–M) statistic was used to test whether the intra- and interspecific genetic distance classes were distinguishable by considering the distance between their means and the distribution of values around the means (Dabboor et al. 2014). The J–M distance is asymptotic to 1.414, with higher values indicating that intra- and interspecific genetic distances are statistically separable (Trigg and Flasse 2001). After confirming the presence of a barcode gap, primers were designed manually for regions flanking the selected mini-barcodes for each of the four groups. When high variability of flanking regions prevented the design of a single primer pair for amplification of all high-level target taxa, multiple primers were designed to form a primer cocktail, specifically for chaetognaths, where two reverse primers were designed. To avoid inadvertent amplification of nuclear mitochondrial insertions (numts), the specificity of each primer was assessed using NCBI’s Primer-BLAST tool, confirming that the primers exclusively targeted mitochondrial genes of the intended taxa (https://www.ncbi.nlm.nih.gov/tools/primer-blast/; Suppl. material 1: table S1).

To evaluate the effectiveness of the mini-barcode primers for copepods (271 bp), euphausiids (310 bp), chaetognaths (283 bp), and hydrozoans (266 bp) (Table 1), amplification success was tested against the universal COI primer set of 658 bp (Folmer et al. 1994) and the miniCOI barcode (Leray et al. 2013) using Sanger sequencing. DNA from 4–5 individual voucher specimens for each of the four groups was extracted using the Zymo Quick-DNA Universal Kit (Zymo Research). An overnight incubation step at 55 °C was added to the standard protocol to ensure complete digestion of tissue in lysis buffer and proteinase K.

PCR reactions (25 μl) contained 20 ng/μl genomic DNA, 12.5 μl OneTaq Quick-Load Master Mix (1X, BioLabs, New England), 0.50 μl forward and reverse primer (10 nmol/L), 7.5 μl sterile nuclease-free water, 1 μl MgCl2 (25 μmol/L), and 1 μl bovine serum albumin (BSA; 1 mg/ml). For primer cocktails, the 0.50 μl primer volume was divided by the number of primers for each forward and reverse reaction; for example, for two reverse primers, 0.25 μl was added for each. Thermal cycling for COI and the newly designed mini-barcode primers comprised an initial denaturation at 94 °C for 2 minutes, 35 cycles of denaturation at 94 °C for 30 seconds, annealing at different primer-specific temperatures for 30 seconds, and extension at 68 °C for 1 minute. The final extension was at 68 °C for 5 minutes. Annealing temperatures were 50 °C for the COI primer, 52 °C for copepods, 56 °C for both euphausiids and hydrozoans, and 53 °C for chaetognaths. A “touchdown” PCR was used for the miniCOI barcode (Leray et al. 2013) to minimize the probability of non-specific amplification. The conditions comprised an initial denaturation at 94 °C for 2 minutes, followed by 16 cycles of denaturation for 10 seconds at 94 °C, annealing for 30 seconds at 62 °C (−1 °C per cycle), extension for 60 seconds at 72 °C, and then an additional 25 cycles at an annealing temperature of 46 °C. The final extension was carried out at 72 °C for 5 minutes. All PCRs included a no-template negative control.

PCR products were visualized on a 2% (w/v) TBE agarose gel containing SafeView™ Classic (Applied Biological Materials Inc., Cat. No. G108). Amplicon size was determined using a 100 bp molecular weight marker (BioLabs, New England). PCR bead clean-up and Sanger sequencing were performed at the Central Analytical Facilities (CAF) at the University of Stellenbosch (South Africa). Sequences were edited using Geneious Prime 2025.0.3, and the percentage of nucleotides with Phred quality scores of at least 30 was calculated for each sequence. Phred scores indicate the probability that a nucleotide has been correctly identified, serving as a standard measure of sequence quality. A Phred score of 30 corresponds to a 1 in 1000 chance of an incorrect base call and is commonly used as a threshold for high-quality sequence data. The nucleotide BLAST tool (BLASTn) on GenBank and BOLD was used for species identification, with a 97% sequence identity threshold. Sequences generated using the newly designed primers were translated into amino acid sequences using the ExPASy Translate tool (https://web.expasy.org/translate/) and screened for stop codons to verify that nuclear mitochondrial insertions (numts) were not inadvertently amplified. The resulting amino acid sequences were subsequently queried against GenBank using the protein BLAST tool (BLASTp) to confirm that only mitochondrial genes were amplified.

To evaluate the versatility of the newly designed primers (Table 1) as well as those chosen from literature for Malacostraca (lobsters, crabs, shrimps, prawns, and isopods; Govender et al. 2019, 2022), Actinopterygii (ray-finned fishes; Folmer et al. 1994; Ward et al. 2005; Leray et al. 2013), and the universal primer set (Folmer et al. 1994; Leray et al. 2013), primer testing was carried out on bulk zooplankton samples that were collected with two vertical (200 μm mesh) and two oblique (500 μm mesh) Bongo net tows off eastern South Africa in July 2022 (36.63°S; 23.75°E). For details of previously published primers, refer to Govender et al. (2022).

Whole zooplankton samples were preserved in 97% ethanol during sampling, with the ethanol replaced after 24 hours for long-term storage. In the laboratory, individual samples were homogenized in the 97% ethanol solution for 1 minute using a consumer blender (Milex; 1500 W at 22,000 rpm). Between samples, the blender was washed to remove residual material and rinsed with a 10% bleach solution and 70% ethanol to degrade any remaining DNA. Three subsamples were taken from each homogenate to improve diversity estimates (Lanzén et al. 2017). Each subsample (10 ml of zooplankton homogenate) was centrifuged at 3000 rpm for 1 minute, and excess ethanol was removed; thereafter, 40 mg of tissue was transferred to a sterile tube. The remaining DNA extraction process was carried out as described in Govender et al. (2023).

PCRs were performed in triplicate to reduce the effects of stochasticity, improve accuracy, and minimize bias and amplification errors (Dopheide et al. 2019). PCR reactions (15 μl) contained 20 ng/μl genomic DNA, 7.5 μl Q5 high-fidelity Taq, 0.6 μl forward and reverse primers (10 μM; Table 1), 3.3 μl sterile nuclease-free water, 0.6 μl MgCl2 (25 μM), and 1.2 μl BSA (1 mg/ml). Where primer cocktails were used, the 0.6 μl primer volume was divided by the number of primers for each forward and reverse reaction.

All primers used the same thermal cycling program: initial denaturation at 98 °C for 30 seconds, denaturation at 98 °C for 10 seconds, annealing at different primer-specific temperatures for 30 seconds, and extension at 72 °C for 30 seconds. The final extension step was carried out at 72 °C for 4 minutes. Annealing temperatures were 46 °C for malacostracans, fishes, and universal primers. The annealing temperatures for the newly designed primers were as per the primer testing section above. All PCRs included a no-template negative control. PCR products were visualized on a 2% (w/v) TBE agarose gel containing SafeView™ Classic. Amplicon size was determined using a 100 bp molecular weight marker. The triplicate PCR products for each of the seven primer sets were pooled and quantified using a Qubit 2.0 Fluorometer (Life Technologies, California, USA), and each of the seven different amplicon pools was further consolidated into a single sample per tow-net haul to create four libraries (one per individual tow) with equimolar concentrations (5 ng/μl). Each library was cleaned using 1.8× Ampure XP purification beads (Beckman Coulter, High Wycombe, UK). Index PCR was performed using the Nextera XT Index Kit (Illumina, San Diego, USA). Thereafter, libraries were cleaned using 0.6× Ampure XP purification beads (Beckman Coulter) and quantified using the Qubit dsDNA High Sensitivity assay kit on a Qubit 4.0 instrument (Life Technologies). The four libraries were sequenced on the Illumina MiSeq platform using the MiSeq Nano Reagent Kit v.2 (500 cycles), following the protocols described by Govender et al. (2023) at the KwaZulu-Natal Research and Innovation Platform (KRISP), South Africa.

The DADA2 algorithm (Callahan et al. 2016) implemented in QIIME2 v.2019.10 (Bolyen et al. 2019) was used to conduct quality control checks, chimera removal, filtering, trimming of primers, truncation of forward and reverse reads, and merging of the paired-end reads into amplified sequence variants (ASVs). The first 20 bases were trimmed, and the sequences were truncated at 200 bp for both forward and reverse reads to exclude low-quality regions. The ASVs were queried against BOLD and NCBI GenBank in August 2024 using a 97% sequence identity threshold. ASVs assigned to the same species were merged manually using MS Excel. Species detected by metabarcoding were cross-referenced with occurrence records obtained from online databases such as the World Register of Marine Species (WoRMS; https://www.marinespecies.org), the Ocean Biodiversity Information System (OBIS; https://obis.org), the Global Biodiversity Information Facility (GBIF; https://www.gbif.org), and online literature.

The most informative region of the COI gene was identified for each of the four groups. The smaller SWAN window sizes (120–140 bp for copepods and hydrozoans, 100–140 bp for euphausiids, and 100–130 bp for chaetognaths) had higher mean K2P distance and lower zero non-conspecific values. Larger window sizes (260–300 bp for copepods, 270–300 bp for euphausiids, 220–300 bp for chaetognaths, and 250–300 bp for hydrozoans) showed better congruence of NJ trees and generated lower K- and R–F scores when compared with reference trees. The intraspecific K2P pairwise distances ranged from 0.00 to 0.14 for copepods, 0.00 to 0.13 for euphausiids, 0.00 to 0.35 for chaetognaths, and 0.00 to 0.16 for hydrozoans. The interspecific distances ranged from 0.00 to 0.67 for copepods, 0.00 to 0.32 for euphausiids, 0.00 to 0.52 for chaetognaths, and 0.00 to 0.79 for hydrozoans. DNA barcode gap analyses based on in silico data showed minimal overlap between intra- and interspecific K2P pairwise distances in all four groups for the selected mini-barcode regions (Fig. 1). All J–M distances were > 1.414, providing statistical support for the existence of a DNA barcode gap in all four datasets.

Figure 1.

Frequency distribution of intra- and interspecific pairwise K2P genetic distances calculated using the selected mini-barcode regions for (a) Copepoda, (b) Euphausiacea, (c) Chaetognatha, and (d) Hydrozoa. The frequency data (Copepoda = n/120000, Euphausiacea = n/1200, Chaetognatha = n/100, Hydrozoa = n/20000) was normalized to obtain a range between 0 and 1.

Primers were designed within conserved regions flanking the mini-barcode regions (Table 1; Suppl. material 1: table S1). The mini-barcode regions were successfully amplified for all four target taxa using the designed primers, with each primer set tested across the respective samples. BLAST and BOLD search results confirmed that the mini-barcode sequences matched the morphologically identified specimens to at least a family level (93–100% sequence similarity). The amplification success rate of the newly designed mini-barcode primer sets was similar to those of the standard COI primer set and the universal miniCOI primer set (Fig. 2). The sequence quality of the new mini-barcode primers outperformed the standard and universal miniCOI primers by consistently recovering higher Phred quality scores of resulting sequences (Fig. 3). The final primer cocktails used for metabarcoding (Table 1) were based on the amplification and sequencing success of the voucher specimens.

Figure 2.

Amplification of mini-barcode primers against COI and miniCOI primers for voucher specimens of (a) Copepoda, (b) Euphausiacea, (c) Chaetognatha, and (d) Hydrozoa. Lane 1 shows the 100 bp molecular weight marker. Lanes 2–6 (a, b) and 2–5 (c, d) display PCR products from DNA amplified using universal COI primers. Lanes 8–12 (a, b) and 7–10 (c, d) show the same samples that were amplified with miniCOI barcode primers. Lanes 14–18 (a, b) and 12–15 (c, d) show amplification with taxon-specific mini-barcode primers. All PCRs included a no-template negative control in the test gel, although this is not reflected in the gel shown above.

Figure 3.

Box-and-whisker plot comparing Phred quality scores of COI, mini-barcodes, and universal miniCOI primers for each taxon.

In total, four zooplankton community libraries were sequenced with Illumina MiSeq. Sequencing was efficient, requiring minimal filtering for both forward and reverse reads during the merging of the paired-end reads for all four zooplankton libraries. For the four libraries, a total of 712,084 read counts were consolidated into 126,449 merged reads, of which 1,456 ASVs were available for analysis across all groups amplified. These ASVs were assigned to a species level at > 97% sequence similarity against reference sequences on BOLD or GenBank and subsequently collapsed into 220 species (Table 2; Suppl. material 1: table S2; Suppl. material 3).

Species richness was greatest for copepods (80 spp.; 36%), followed by euphausiids (25 spp.; 11%), molluscs (24 spp.; 11%), fishes (23 spp.; 11%), hydrozoans (21 spp.; 10%), decapods (14 spp.; 6%), and chaetognaths (12 spp.; 5%; Fig. 4a). Taxon-specific primers were used for all these groups except for molluscs. Groups for which taxon-specific primers were not used had lower species richness in metabarcoding outputs (<3% per taxon). These taxa included amphipods, echinoderms, ostracods, and scyphozoans. Copepods (41%) and euphausiids (35%) dominated the read counts as a proxy for relative biomass (Fig. 4b), followed by decapods (10%). Chaetognaths (4%) and hydrozoans (5%) had more read counts than any of the other groups for which taxon-specific primers were not used (<2% of read counts per group).

Figure 4.

Bar graphs showing (a) overall species richness and (b) overall read counts. Blue bars indicate taxa for which specific primers were used.

The authors have declared that no competing interests exist.

No ethical statement was reported.

No use of AI was reported.

This work was supported by the National Research Foundation (NRF).

Conceptualization: AG, JG. Data curation: AG, SH, JG. Formal analysis: JG, SH, AG. Funding acquisition: AG, JG. Investigation: JG, AG. Methodology: JG, AG. Project administration: AG, JG. Resources: JH, AG. Software: AG. Supervision: SWM, AG, JG. Validation: AG, SH. Visualization: AG. Writing - original draft: SH. Writing - review and editing: AG, SWM, JG, JH.

Saiesha Harpal https://orcid.org/0000-0002-6628-1453

Johan Groeneveld https://orcid.org/0000-0002-9831-9073

Sandi Willows-Munro https://orcid.org/0000-0003-0572-369X

Jenny Huggett https://orcid.org/0000-0001-9315-8672

Ashrenee Govender https://orcid.org/0000-0002-2860-4610

All raw sequence reads used to perform analyses have been uploaded to NCBI under accession number PRJNA1274484.