Two primer sets were developed using PrimerMiner (Elbrecht and Leese 2016) and a previously generated sequence alignment of 15 bioassessment relevant freshwater macroinvertebrate groups (Elbrecht and Leese 2017). The novel fwh1 and fwh2 primer sets amplify a short region of the cytochrome c oxidase I (COI) region of 178 and 205 bp in length respectively (Fig. 1, A). Both primer sets were based on primer sequences previously published in the literature (Table 1, Suppl. material 1), but primer degeneracy was increased to better match freshwater invertebrate taxa. For sequencing, the primers were ordered to include Illumina tails and individual inline barcodes for multiplex sequencing on the MiSeq system (Suppl. material 2, see Elbrecht and Leese 2015 for details on the “fusion primer” method). Using a 6 bp inline barcode for demultiplexing, the developed fusion primers can be used to individually tag up to 36 samples per primer set (Suppl. material 3).

Figure 1.

Developed primer sets and samples sequenced for primer validation. Two independent PCR replicates were run and sequenced for each sample. A: Binding sites of the two primer sets (fwhF1+fwhR1 and fwhF2+fwhR2) targeting a 178 and 205 bp fragment internal to the COI Folmer barcoding region (Folmer et al. 1994). The fwhR2 primer was affected by a design error, thus an improved version (fwhR2n) was developed. B: Overview of the sequenced benthic communities and mock samples to test and validate the developed primer sets. Five mock communities (four multi and one single species) from Elbrecht and Leese (2015), as well as three kick samples collected from streams in Romania (Călățele River: L2, Almaşul River: R2, Valea Racilor River: Z2), were collected and tested using the fwh1 and fwh2 primer sets (except for sample DceM that could only be amplified using the fwh1 primer set).

To explore the broader performance of the newly developed primers compared to the commonly used primers sets (Suppl. material 1), all primers were evaluated in silico for insect groups (following the taxonomy by Misof et al. 2014). Insect COI reference data was obtained in April 2016 and clustered into OTUs from NCBI and BOLD using PrimerMiner v0.3 as described previously (Elbrecht and Leese 2016, Elbrecht and Leese 2017). Sequence alignments were generated and used to evaluate the penalty scores given for primer mismatches using PrimerMiner v0.13 with the default settings (mm_position ="Position_v1", mm_type = "Type_v1"). Only orders with at least 100 OTUs were used to calculate the average penalty score for the respective primer.

The performance of the fwh1 and fwh2 primer sets was evaluated using four previously used mock communities each containing 52 different freshwater taxa (sample A, B, C and D) and one single species mock sample with 31 specimens with unique haplotypes and known biomass (Elbrecht and Leese 2015). Additionally, kick samples from three Romanian rivers (Almaşul, Călăţele and Valea Racilor, Suppl. material 4) were analyzed using both primer sets. The kick samples were collected in fall 2016, preserved in 95% ethanol and stored at -20°C for later molecular analysis. For the kick samples, no morphological identification of the macroinvertebrates was performed. Prior to DNA extraction, specimens were size sorted into small (S, body size < 2.5 x 5 mm), medium (M, up to 5 x 10 mm) and large (L, maximum size of 10 x 20 mm) specimens (Suppl. material 5, also see Elbrecht et al. 2017a).

Specimens of each size category (S, M & L) were dried overnight in sterile Petri dishes to remove the ethanol. Specimens from each category were homogenised using an IKA ULTRA-TURRAX Tube Drive control system (IKA, Staufen, Germany) with sterile 20 mL tubes and 10 steel beads (5 mm Ø) by grinding at 4000 rpm for 30 minutes. Approximately equal amounts of grinded tissue from each category were digested following a modified salt DNA extraction protocol (on average 13.41 mg of tissue, SD = 12.34 mg, Sunnucks and Hales 1996, Elbrecht et al. 2017a). Next, the lysate was pooled proportionately to the abundance of individuals in each size category to reduce the overrepresentation of large specimens (see Elbrecht et al. 2017a for details). Further, 20 μl of the extracted DNA from each respective sample was digested with 1 μL RNase A (10 mg/mL, Life Technologies, Darmstadt, Germany) and cleaned up using a MinElute Reaction Cleanup Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. DNA concentrations were quantified fluorometrically using Qubit fluorometer (HS Kit, ThermoFisher Scientific, MA, USA) and concentrations for all samples were adjusted to 25 ng/μL for PCR.

The five mock communities and three kick samples were amplified in duplicates in a one-step PCR using the developed freshwater primers (Table 1). Suppl. material 6 gives an overview of fusion primer combinations used to uniquely tag each sample. Each PCR reaction was composed of 1× Multiplex PCR Master Mix (Qiagen Multiplex PCR Plus Kit, Qiagen, Germany) 0.5 μM of each primer, 25 ng DNA, filled up with HPLC H₂O (Carl Roth GmbH, Karlsruhe, Germany) to a total volume of 50 μL. PCR reactions were run in a Biometra TAdvanced Thermocycler (Biometra, Göttingen, Germany) using the following program 95°C for 5 min, 34 cycles of 95°C for 30 s, 52°C (for the fwhF1+fwhR1 primer pair) or 58°C (for the fwhF2+fwhR2 primer pair) for 30 s, 72°C for 2 min, and 72°C for 10 min. The annealing temperatures for both primer sets were established by first running a gradient PCR on DNA from the multi species mock communities (gradient temperature 43.7 - 70.3°C, Suppl. material 7). The annealing temperatures for the respective primer pair were chosen a few degrees below the temperature of the last visible band, to ensure efficient and consistent amplification. PCR products from the one-step PCR were purified and left size selected using SPRIselect (Beckman Coulter, CA, USA) with a ratio of 0.76x and the DNA concentration was quantified with a Qubit fluorometer, High Sensitivity Kit (Thermo Fisher Scientific, MA, USA) and Fragment Analyzer Automated CE System using NGS Standard Sensitivity kit (Advanced Analytical, Heidelberg, Germany). The mean DNA concentration from both measurements were used to pool PCR products by equal molarity. This final library was additionally purified with the MinElute Reaction Cleanup Kit (Qiagen, Hilden, Germany), as a precaution due to BSA interfering with the PCR clean-up using SPRIselect (Elbrecht et al. 2017a). Sequencing was done on two runs of an Illumina MiSeq system using a 250 bp paired end read kit (v2) and 5% PhiX spike-in. Sequencing was carried out by GATC Biotech GmbH (Konstanz, Germany). Raw sequence data were processed using a modified version of the UPARSE pipeline (Edgar 2013, v9.0.2132), which is available at GitHub (JAMP version 0.17 - https://github.com/VascoElbrecht/JAMP/). The exact commands run to reproduce the analysis are available in Suppl. material 8. In short, reads were demultiplexed, paired-end merged using usearch, reverse complement sequences generated where necessary, quality filtered (maxee = 0.05) and pre-processed (primer removal, Cutadapt v1.9 (Martin 2011), discarding of reads +/- 10 bp of the expected length, dereplication with removing singletons; minsize = 2). Before applying clustering (97% similarity) all retained sequences for the A, B, C and D samples were pooled and the three Romanian samples were also pooled. Reads including singletons were remapped against the OTUs and clusters with at least 0.003% abundance in one sample retained (both replicates). OTUs were identified using sequences from previous studies as references and comparison against BOLD and NCBI reference databases with JAMP. For the single species mock samples (DceM) filtered dereplicated reads of exact 178 bp length were directly mapped against the expected haplotypes (Suppl. material 8) and all matching hits with at least 0.003% abundance retained.

Two primer sets were developed targeting short COI fragment lengths of 178 bp and 205 bp respectively (Fig. 1). In silico evaluation of the developed primer sets on insect orders was only carried out when preparing this manuscript, and it became evident that the fwhR2 primer had a design flaw (Fig. 2). At position 9 from the 3' end (Table 1), an Adenine (A) was used instead of a Thymine (T), leading to a poor estimated primer performance (mean penalty score of 144.7). This mistake was corrected afterwards in the fwh2n version of this primer (Table 1), which shows a decreased average penalty score (53.4). However, the new improved version of the primer was only tested in silico and all laboratory tests were carried out using the flawed fwhR2 version of this primer.

Figure 2.

In silico evaluation of insect groups (after Misof et al. 2014) for selected metabarcoding primer pairs. COI reference sequences for primer evaluation were obtained from BOLD and NCBI using PrimerMiner and processed into OTUs (3% similarity). For primer-template mismatches, penalty scores were calculated using PrimerMiner (lower penalty score = better expected primer performance). The individual mean penalty scores are given in bar plots for each primer and insect order. The average penalty score was calculated for each primer for orders with at least 100 OTUs for the respective primer pair. The typically used primer combinations are indicated by connected grey lines on the left, as well as black lines for the newly developed primer pairs. Templates for primer development are indicated with blue boxes, while the newly developed primers are highlighted with blue backgrounds. The fwhR2 primer had a design error and is highlighted in grey.

The other evaluated primer sets showed mixed performances depending on the degeneracy of the respective primers. A lack of degeneracy resulted in rather high penalty scores, as was the case for the LCO1490+HCO2198 and ZBJ-ArtF1c+ZBJ-ArtR2c primer sets (scores above 100, Zeale et al. 2010, Folmer et al. 1994). Primers incorporating an Inosine, e.g. ArF5+ArR5 (Gibson et al. 2014), or a high degeneracy, e.g. the BF+BR primer sets (Elbrecht and Leese 2017), showed low average penalty scores (below 40). The universal BF+BR and mlCOIintF primers showed increased penalty scores for a few groups that have more variable primer binding regions in the template DNA (Thysanoptera, Phasmatodea or Raphidioptera). Some of the primers binding at the LCO1490 binding site showed high scores due to misaligned sequences or low number of OTUs. The newly designed primers fwhF1, fwhF2 and fwhR1 had lower penalty scores than the primer sets they are based on, while the fwhR2n primer set showed a higher penalty score (53.4) than the ArR5 primer set with a score of 6.9.

Both fwh primer sets successfully amplified the four multispecies mock communities (A, B, C and D) as well as the three Romanian stream kick samples. The fwh2 primers only produced a weak amplicon band on the agarose gel for the DceM sample, which was therefore only sequenced using the fwh1 primer set. Both Illumina MiSeq runs were successful for all sequenced samples with an average number of 1.40 (fwh1) and 0.74 (fwh2) million sequences obtained for each replicate (SD = 0.26 and 0.13, Suppl. material 9). Raw sequencing data are available on NCBI SRA (SRR5295658 and SRR5295659). OTU tables including assigned taxonomy and OTU sequences are available as supporting information (Suppl. material 10).

For the four mock communities, most of the taxa were recovered by both primer sets. While fwh1 primers detected 48 taxa out of 52, the fwh2 performed poorer, recovering 46 taxa (Fig. 3). Variation in logarithmic taxa read abundance was much lower for the fwh1 amplicons (SD = 0.62) than for the fwh2 primer set (SD = 0.97) across the mock community samples (Table 2). The fwh1 primer set also showed the highest precision (deviation from expected read abundance). For DceM mock community, only the fwh1 primer set produced an amplicon, as the fwh2 primer set did not amplify Perlidae efficiently (see also Fig. 3). Because the fwh1 fragment is shorter than the previously sequenced Folmer COI fragment (Elbrecht and Leese 2015), only 15 of the original 31 haplotypes could possibly be distinguished (Fig. 4). All 15 expected haplotypes in the DceM community were recovered with the fwh1 primer set. Both PCR replicates showed the same trend in the relative sequence abundance with an expected ratio of relative haplotype abundance approximately equal to 1 (Fig. 4, B). Also, both replicates had very similar read composition, with only rare reads being unique to specific samples (Suppl. material 11).

Figure 3.

Comparison of fwh1 (A) and fwh2 (B) primer performance, both tested with the same four bulk samples with two independent PCR replicates for each sample. Each respective sample contained 52 morphologically distinct macroinvertebrate taxa ("TierMix": A, B, C & D). The 52 taxa are shown on the x-axis with the number of reads obtained for each morphotaxon indicated by black dots on the logarithmic y-axis (mean relative abundance of detected morphotaxa is indicated by red circles, replicates are plotted). Sequence abundance was normalised across the samples and the amount of tissue used in each DNA extraction. Only OTUs which had a minimum abundance of 0.003% in at least one of the four samples were included in the analysis. Number of samples for which a morphotaxon was not detected is indicated by orange and red numbers in each plot. A thick vertical line in light red indicates if a morphotaxon was not detected.

Figure 4.

Detection of haplotypes in the tested single species mock community (DceM) using the fwh1 primer set. Sequences below 0.003% relative read abundance were discarded. A: Relative abundance of detected haplotypes in both PCR replicates plotted against cumulative specimen weight (red line indicates linear regression). Because the fwh1 fragment is shorter than the previously sequenced Folmer COI fragment (Elbrecht and Leese 2015), only a maximum of 15 haplotypes can be detected with the short COI fragment. B: Ratio of relative haplotype abundance when dividing replicate A by replicate B with a red line indicating the expected value of 1.

For the three Romanian samples the ecological quality state of the rivers was assessed only on the expert judgment (visual assessment, Suppl. material 4) and not based on a standardised assessment using morphologically identified macroinvertebrate taxa from kick samples (see Suppl. material 5 for pictures of samples composition). However, by analyzing the diversity of EPT taxa (Ephemeroptera, Plecoptera, Trichoptera, typically highly pollution-sensitive taxa), it is possible to get a proxy for the ecological condition of the streams. For the study sites L2 and Z2 (good to mediocre ecological status according to expert judgment) 42.66% and 46.44% of the OTUs were identified as EPT, while at the R2 site (poor ecological state) only 7.82% EPT taxa were detected (Suppl. material 12). For fwh2 primer set we obtained very similar results; for Z2 and L2 sites the EPT is represented by 18.53% and 29.23%, while for R2 site 8.47% of the OTUs were assigned to EPT taxa. Taxonomic richness of the streams communities is in good agreement with their ecological state. Our primer pairs also amplified non-target species, with high identity matches (>= 97%) to the reference databases, such as hop aphids, moths and few freshwater fish species (e.g. gudgeon, minnows and stone loaches). The principal component analysis of the macroinvertebrate OTUs obtained from the fwh1 and fwh2 primer sets showed clear differentiation between the three Romanian samples, while consistently grouping PCR replicates of the same sites together (Fig. 5).

Figure 5.

Principal component analysis (PCA) of freshwater macroinvertebrate OTUs detected with the fwh1 (A) and fwh2 (B) primer set in the three Romanian river samples. PCR replicates of the identical samples are shown with the same colour.

Using PrimerMiner we have developed two short universal metabarcoding primer sets targeting freshwater macroinvertebrates. As previously reported, a short 150 bp barcode marker is sufficient to identify most insect taxa on species level (Meusnier et al. 2008). Also, PCR with short amplicons is expected to work better when dealing with highly degraded DNA (Dalvin et al. 2010, Mitchell 2015, Schäffer et al. 2017). Additionally, in contrast to previously developed longer universal markers like the BF2+BR2 primer set (Elbrecht and Leese 2017), fragments of ~200 bp length can be paired-end sequenced on the Illumina NextSeq system increasing throughput ten-fold compared to the MiSeq/HiSeq system, which is commonly used for amplicon sequencing (e.g. Schöfl et al. 2017). The MiniSeq system can also be used if only a few samples have to be sequenced. Additionally, completely overlapping amplicons can reduce sequencing errors when paired-end merged (Kozich et al. 2013, Eren et al. 2013).

In silico evaluation of the newly developed primer sets revealed that all of them showed low penalty scores, except for the fwh2R primer, where a design error was introduced causing mismatches across most taxa. Unfortunately, in silico evaluation was only carried out when preparing this manuscript, therefore a corrected version of the primer (fwhR2n) could not be tested in vitro with mock and kick samples. Compared to the fwh1 primer set, the fwhR2n primer still shows a high average penalty score of above 50, thus even the improved reverse primer version is not likely to perform particularly well for amplification of insects. The in silico evaluation also showed that the fwh1, BF+BR as well as the ArF5+ArR5 primer sets (Gibson et al. 2014, Elbrecht and Leese 2017) are likely to work well across most insect orders.

Both the fwh1 and the fwh2 primer sets (using the flawed fwhR2 primer) were additionally tested on several mock communities as well as complete kick samples. While both primer sets could clearly differentiate the three stream sites in a principal component analysis, the fwh2 primer set showed a higher primer bias when amplifying mock communities, each containing 52 freshwater macroinvertebrates. The primer set also failed to amplify the DceM mock community containing specimens of a single stonefly and was thus not included in sequencing. Perlidae specimens were also underrepresented in the multi species mock community with the fwh2 primer set, indicating strong primer bias for this group. Both the fwh1 and fwh2 primer sets detected 98% of the freshwater insects present in the multi species mock communities. While both universal primer sets show higher detection rates and reduced primer bias compared to the standard COI Folmer primer sets (Folmer et al. 1994, Elbrecht and Leese 2015), the recently developed BF2+BR2 primer set for a longer fragment shows less primer bias than the here proposed novel primer sets (Elbrecht and Leese 2017). Therefore, as long as DNA degradation is not a concern we recommend the use of the BF2+BR2 primer sets for DNA metabarcoding of freshwater macroinvertebrate and insect samples (except for the derived Thysanoptera group). The ArF5+ArR5 primer set also performed well in the in silico evaluation but should be further validated in vitro. Additionally, when dealing with highly degraded DNA, especially the fwh1 primer set might prove to be useful.

We also used the new fwh1 primer set to test for a linear relationship of sequence abundance to specimen biomass within species, which we previously explored for the Folmer primer sets (Folmer et al. 1994, Elbrecht and Leese 2015). In the previous study, a single mock sample containing 31 unique haplotypes of a single stonefly species was amplified and a significant linear relationship between numbers of sequences and specimen biomass was detected. While this sounds promising for estimating taxa abundance or biomass from metabarcoding data, reliable estimates are difficult to obtain due to often severe primer bias between different species (Piñol et al. 2014, Elbrecht and Leese 2015). However, the sequence / biomass relationship within species was tested here again using the shorter fwh1 primer set that has higher degeneracy and can be paired end merged with a complete overlap of forward and reverse sequencing reads, potentially reducing sequencing errors compared to the previous study where the standard Folmer primers were used (Folmer et al. 1994). As expected, with the fwh1 primer set also a significant correlation between read abundance and specimen biomass was detected. In addition, by comparing the two sequenced PCR replicates, we could analyse the false positive haplotypes generated by sequencing errors and chimeras. In particular abundant haplotypes showed hundreds of artificial haplotypes, likely derived from sequencing errors and chimera formation. These errors were consistent between PCR replicates, with only low abundant sequences being unique to the respective sample (Lange et al. 2015). Both PCR replicates showed differences in low abundant reads, which were likely to be generated by sequencing errors or chimeric sequences. However, the majority of the false positive haplotypes were shared between both replicates indicating systematic origins, for example, chimera formation or sequencing errors on high abundant haplotypes. Amplification of abundant taxa is typically very consistent, with stochastic effects mostly affecting low abundant taxa / sequences as also demonstrated in other studies (Leray and Knowlton 2017), highlighting that the use of PCR replicates might not substantially increase the reliability of DNA metabarcoding results (Smith and Peay 2014). Some of this false positives, however, might also be present due to mitochondrial heteroplasmy or the presence of 'numts', i.e. nuclear sequences of mitochondrial origin (Bensasson 2001).