Research Article
Print
Research Article
Direct PCR meets high-throughput sequencing – metabarcoding of chironomid communities without DNA extraction
expand article infoNina Röder, Klaus Schwenk§
‡ RPTU Kaiserslautern-Landau, iES - Institute for Environmental Sciences, Landau, Germany
§ LOEWE Centre for Translational Biodiversity Genomics, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
Open Access

Abstract

Aquatic emergent insect communities form an important link between aquatic and terrestrial ecosystems, yet studying them is costly and time-consuming as they are usually diverse and superabundant. Metabarcoding is a valuable tool to investigate arthropod community compositions, however high-throughput applications, such as for biomonitoring, require cost-effective and user-friendly procedures. To investigate if the time-consuming and labour-intensive DNA extraction step can be omitted in metabarcoding, we studied the difference in detection rates and individual read abundance using standard DNA extraction versus direct PCR protocols. Metabarcoding with and without DNA extraction was performed with artificially created communities of known composition as well as on natural communities both of the dipteran family Chironomidae to compare detection rates, individual read abundances and presence-absence community composition. We found that the novel approach of direct PCR metabarcoding presented here did not alter detection rates and had a minor effect on individual read abundances in artificially created communities. Furthermore, presence-absence community compositions of natural chironomid communities were highly comparable using both approaches. In conclusion, we showed that direct PCR protocols can be applied in chironomid metabarcoding approaches, with possible application for a wider range of arthropod taxa, enabling us to study communities more efficiently in the future.

Key words

cytochrome c oxidase I (COI or COX1), DNA isolation, metabarcoding, next-generation sequencing (NGS), presence-absence community composition, read abundance, size-sorting

Introduction

A key global challenge in the 21st century is the attempt to reverse the ongoing global biodiversity decline and to mitigate its consequences. DNA metabarcoding for large-scale monitoring of species rich and abundant groups, such as insects, is labour- and cost-intensive, but has become more and more achievable since its development in the early 2000s (Gostel and Kress 2022; Guo et al. 2022). One typical processing step during the preparation of insect metabarcoding samples is the extraction of DNA from homogenized bulk samples (Yu et al. 2012; Creer et al. 2016). Although it depends on the choice of protocols or commercial kits and the available capacities to automate processes, DNA extractions are often costly and time-consuming.

Typically, a laborious DNA extraction step prior to PCR is performed in order to purify and concentrate DNA from tissue samples. However, tissue from different invertebrate species has been successfully added to PCR reactions without prior DNA extraction in so-called direct PCR (dPCR) approaches (e.g., flies and starfish: Wong et al. 2014; mosquitoes: Werblow et al. 2016; spider mites: Sakamoto and Gotoh 2017; ants: Wang et al. 2018). Through the development of high-performance enzymes and buffers, PCR has become more robust against inhibitors, reducing the need to eliminate them via elaborate DNA extraction procedures. In addition, current high-throughput sequencing (HTS) technologies, such as sequencing by synthesis or nanopore sequencing, are very efficient in detecting even trace amounts of DNA. It has further been shown that fragments amplified by dPCR can be used in HTS approaches (e.g., in generating barcodes for individual chironomids: Baloğlu et al. 2018 and individual ants: Wang et al. 2018), although, to the best of our knowledge, this has neither been conducted for animal communities nor systematically tested and compared to conventional approaches before.

To investigate the suitability of direct PCR protocols in insect community metabarcoding we compared purified DNA with direct application of homogenized tissue as a substrate for PCR and subsequent metabarcoding. We chose communities of the family Chironomidae, since they are a superabundant and very species rich group of merolimnic insects. Although chironomids are important components of biomonitoring programmes worldwide, their morphological identification to species level is very difficult. Therefore, they are suitable target organisms for developing efficient metabarcoding techniques. We subjected artificial (known species composition) and natural communities (highly variable species composition) to standard vs. dPCR metabarcoding using previously established cytochrome c oxidase I (COI) markers (Elbrecht and Leese 2017). We compared the detection rates and individual read abundances of artificial communities generated using the standard and the dPCR metabarcoding approach. Based on the findings of Elbrecht and Leese (2015), we expected prior size-sorting to have an effect on both detection rates and read abundance, as our chironomid specimens varied in their individual mass up to a factor of more than 20. Therefore, we created artificial communities both with and without prior size-sorting. Furthermore, we compared the community composition of highly diverse natural chironomid communities (exposed to environmental stress) comparing both methods. We finally investigated the influence of using different amounts of purified or amplified DNA on read abundance, to assess whether (dPCR) metabarcoding can also be used for quantitative estimates.

Materials and methods

Chironomid origin

Chironomidae were retrieved from artificial ponds of the Eußerthal Ecosystem Research Station (EERES; 49°15′14″N, 7°57′42″E) near Landau, Germany, in 2019 and 2020. Adult specimens were collected from passive emergence traps once to twice a week during spring and summer. In 2020, the artificial ponds were simultaneously used to study the effect of the mosquito control agent Bacillus thuringiensis israelensis (Bti) on merolimnic insect communities (Kolbenschlag et al. 2023).

Size-sorting and tissue preparation

Chironomid samples were sorted into four different size categories with known average weight per specimen (cf. Sabo et al. 2002; Suppl. material 1), to account for differences in individual mass that occur within and among chironomid species. They were subsequently stored in 70% ethanol. For tissue preparation, samples were dried at 60 °C for 16 to 24 hours and then finely ground by a TissueLyser II (Qiagen, Hilden, Germany) at 30 Hz for 3 × 1 min using two stainless steel beads (2–3 mm diameter). PCR-grade water (1 ml) was added to each sample and thoroughly vortexed. The tissue-water mixes were frozen at -20 °C until further analysis.

Sanger-sequencing and identification of individual specimens

Size-sorting of natural chironomid samples led to cases with only one specimen per size group and sample (due to low sample size). All these single specimens were individually Sanger sequenced. We used a direct PCR approach following Wong et al. (2014; see section ‘Metabarcoding’) with either the primer combination LCO1490 & HCO2198 (5’-GGTCAACAAATCATAAAGATATTGG-3’ & 5’-TAAACTTCAGGGTGACCAAAAAATCA-3’; Folmer et al. 1994) or, additionally, C1-J-1751 & C1 N2353 (5’-GGAGCTCCTGACATAGCATTCCC-3’; Simon et al. 1994 & 5’-GCTCGTGTATCAACGTCTATWCC-3’; Lewis et al. 2005) to target overlapping COI sequences of 658 bp and 557 bp length, respectively. The obtained Sanger sequences were aligned and taxonomically assigned using the BOLD Identification System (IDS) for COI (Ratnasingham and Hebert 2007; last accessed on 09/02/2023) by querying against the Public Record Barcode Database for animal identification (Suppl. material 2).

DNA extraction

We used 550 µl of the tissue-water mix to purify DNA from each sample with two technical replicates following an adapted high salt DNA extraction protocol after Aljanabi (1997). As the adapted protocol is designed for DNA extraction from dried tissue, smaller amounts of double concentrated buffer solutions were used in the beginning to account for the excessive water in the tissue-water mixes. Briefly, 450 µl of 2× SEB and 100 µl of 20% SDS buffer were added to 550 µl of tissue-water mix. After splitting into two technical replicates, 5 µl of Proteinase K (10 mg/ml) were added to each sample followed by 1 h of incubation at 60 °C. 350 µl of 5 M NaCl were added and – after centrifuging for 30 min – 600 µl of the supernatant were transferred to a fresh tube. An equal volume of isopropanol was added and the samples stored at -20 °C overnight. Samples were then centrifuged for 20 min at 4 °C. The resulting pellet was washed with 70% ethanol, dried and finally resuspended in 25 μl 1× TE buffer. Extraction blanks were included to ensure data reliability.

Artificial community composition

To assess the detection rates and specimen-specific read abundances under different metabarcoding approaches (Experiments 1–3), we created artificial communities with known composition of different chironomid species. Tissue-water mixes of individual chironomid specimens were selected based on the specimens’ COI sequences. We aimed for taxonomically diverse artificial communities, while still being able to distinguish between specimens by their COI sequences targeted in the metabarcoding approach. We selected tissue-water mixes of 16 specimens from 2019 with at least 1.5% dissimilarity in the targeted region to create one artificial community (“Mock A”, Suppl. material 2). In a second approach, we used tissue-water mixes of 13 specimens from the 2019 sampling campaign with at least 13% dissimilarity in the target region (“Mock B”, Suppl. material 2) to investigate the effect of varying amounts of DNA input in detail. For a schematic overview of the preparation of mock communities see Fig. 1.

Figure 1.

Schematic overview of the workflows during preparation of artificial chironomid communities. Two different sets of chironomid specimens (Mock A and Mock B) were analysed using different approaches of standard and direct PCR metabarcoding, to study the effect of direct PCR or DNA extraction protocols on detection rates and individual read abundances. Experimental setups included analysing the effects of variable mass per individual and different mock compositions (Experiment 1), variable or similar amounts of different input materials (Experiment 2), and variable concentrations of mock material to assess sensitivity (Experiment 3). Mocks were created by pipetting tissue-water mixes (Mocks A-1, A-2 and B-1), purified DNA extracts (Mock B-2) or PCR products (Mock B-3) of individual specimens. A subsample of Mocks A-1 (a, b) and A-2 (c, d) was used for sequential dilution (Mock A-3). Pie charts indicate if artificial communities were created under a size-sorting scenario (aiming for even masses per individual) or under a non-size-sorting scenario (where individual masses naturally vary). Purified DNA extracts were subjected to PCR, while non-purified tissue-water mixes were used in dPCR approaches. The labels ‘Figure 2’ to ‘Figure 5’ within the figure refer to the corresponding graphical representations of the results.

Experiment 1: Effect of size-sorting and tissue preparation on read abundance of specimens in artificial communities

Mock A was used to assess the effect of size-sorting (prior to sample preparations) and two different metabarcoding approaches i.e., with and without DNA extraction, on read-abundance per species. Each of the 5 replicates of Mock A were created by pipetting either the same amount of tissue-water mix (34 µl per specimen, concentrations varying between 0.0001 and 0.0014 mg tissue per 1 µl) simulating a “without size-sorting” scenario (Mock A-1, n = 5) or an adapted amount of tissue-water mix (4–116 µl) with approx. 0.006 mg tissue per each of the 16 specimen (Mock A-2: “with size-sorting” scenario, n = 5). For DNA extraction of the resulting tissue-water mixes, we followed the approach described in section ‘DNA extraction’ (but using only 495 µl tissue-water mix, the rest was needed for direct PCR applications). We further ran a size-sorting plus dPCR approach with Mock B (Mock B-1; n = 3, see Experiment 2 for details) to test if results are consistent when different specimens are used.

Experiment 2: Effect of DNA input variation on read abundance of specimens in artificial communities

With Mock B we assessed the effect of DNA input variation, i.e., varying or equal amount per specimen of either total DNA (Mock B-2) or target fragments (Mock B-3), on read abundance per specimen. DNA was extracted and purified from the 13 individual specimens and quantified using a Qubit fluorometer (Qubit dsDNA HS Assay Kit, Thermo Fisher Scientific). Both equal (5 ng) and varying (0.7–5 ng) amounts of DNA per specimen were then pooled into artificial communities (Mock B-2; n = 4), the latter using each 1.2 µl of the purified DNA extract per specimen. In addition, direct PCR was performed with each of the specimens individually and target fragments were quantified using a Tapestation 4200 (D1000 DNA Screentape Analysis Kit; Agilent Technologies, Santa Clara, CA, USA). Consequently, we pooled either the same (14, 31, 50 or 70 ng) or varying (14–80 ng) amounts of target fragments per specimen into artificial communities (Mock B-3; n = 4) and subjected them to the second PCR for tagging and adding of Illumina adapters.

Experiment 3: Effect of DNA input dilutions on detection rates of specimens in artificial communities

To assess sensitivity of the different methods, we determined specimen detection rates with different dilutions of the mock communities created for Experiment 1, but using only three of the five replicates per artificial mock (Mock A-3; n = 3). Subsamples of artificial communities were sequentially diluted 1:2 and PCR success was checked on a 1.5% TBE agarose gel. Due to limited sample capacity, five to six dilutions were chosen based on the quality of resulting bands to cover a representative range of dilutions that still yielded in PCR success. Further, only one technical replicate was used. We used sequential dilutions of up to 1:32 or 1:64 of the original mock for those resulting from direct PCR with or without prior size-sorting, respectively. For mocks based on purified DNA we chose dilutions between 1:8 and 1:512 (with prior size-sorting) and 1:16 up to 1: 1024 (without prior size-sorting) dilutions of the original mock.

Experiment 4: Natural communities under environmental stress

To assess the applicability of the dPCR approach compared to common metabarcoding protocols on natural chironomid communities, we used a sub­sample (eight out of 12 artificial ponds, five out of 26 sampling dates, N = 40) of an ongoing ecotoxicological study. Half of the eight artificial ponds had been treated with the mosquito control agent Bti (for details see Kolbenschlag et al. 2023), the five sampling dates corresponded to a time period of 17 days in June 2020. The adult specimens sampled were size-sorted and their tissue prepared (see section ‘Size-sorting and tissue preparation’). Tissue-water mixes of different size-groups were pipetted to achieve even masses per specimen in each sample. All samples were then analysed using both metabarcoding approaches, i.e., with and without DNA extraction.

Metabarcoding

We followed a two-step PCR metabarcoding approach using the primers BF2 & BR2 (5’-GCHCCHGAYATRGCHTTYCC-3’ & 5’-TCDGGRTGNCCRAARAAYCA-3’; Elbrecht and Leese 2017) to amplify a 421 bp COI fragment. PCRs included negative controls and two technical replicates per sample. When purified DNA was applied, each of the two purified DNA extracts per sample served as source for one of the technical replicates. For the initial PCR, we used per reaction either 1 µl of purified DNA extract or 5 µl of tissue-water mix, added 2 µl 10× buffer (TaKaRa, Shuzo, Japan), 1.2 µl dNTPs (2.5 mM), 1 µl of each primer (10 µM) and 0.15 μl of Ex Taq DNA Polymerase (TaKaRa, Shuzo, Japan), and filled up with PCR-grade water to achieve a final reaction volume of 20 µl (Wong et al. 2014). For the second PCR, the first-PCR product was 1:20 diluted (to reduce the amount of primers) and 1 µl was used as template. Primers were replaced by the corresponding fusion primers, including Illumina adapters for sequencing (P5 or P7) and inline barcodes of different length for individual tagging of samples (Elbrecht and Steinke 2018). PCRs had following cycling conditions: 94 °C for 5 min, 42 cycles (12 instead of 42 cycles in the second PCR) of 94 °C for 30 s, 51 °C for 30 s and 72 °C for 1 min, ending with 72 °C for 10 min (adapted from Wong et al. 2014). PCR success was checked on a 1.5% TBE agarose gel. The target fragment concentration was quantified using a Tape­station 4200 (D1000 DNA Screentape Analysis Kit; Agilent Technologies, Santa Clara, CA, USA). For each sample, the PCR product was pooled into a library proportionally to the number of specimens to ensure equal sequencing depth per specimen over all samples. Final libraries were purified and size-selected with magnetic beads (ratio: 0.65×, SPRIselect, Beckmann Coulter, Bread, CA, USA), retaining mainly fragments >300 bp. Libraries were sent to an external laboratory (CeGaT, Tübingen, Germany) for 2 × 250 bp (v2) paired‐end sequencing on a MiSeq Illumina system. Samples were included in four different libraries: Sequencing of the library containing Mock A resulted in a total of 13.763.845 reads, sequencing of the library with Mock B resulted in 12.047.150 reads in total. Sequencing of the two libraries containing the natural community samples resulted in 5.928.530 and 9.336.106 reads, respectively. Each of the latter two libraries contained all samples from two treatment and two control ponds.

Bioinformatic analysis

Bioinformatic analysis was performed following the JAMP approach (https://github.com/VascoElbrecht/JAMP, package version 0.77). In short, raw data were demultiplexed and adapters were removed. Consequently, paired end-sequences were merged using usearch (https://drive5.com/usearch, version 11.0.667), then primers were removed using cutadapt (https://github.com/marcelm/cutadapt, version 3.5). After quality filtering, where sequences beyond the target length of 421 ± 10 bp (cutadapt) and more than 1 expected error (usearch) were discarded, operational taxonomic units (OTUs) were clustered with a 3% radius using usearch and vsearch (https://github.com/torognes/vsearch, version 2.21.1).

In the case of Mock A, specimens were too similar in the target COI region to be identified by the OTU clustering approach, i.e., two specimens could be clustered into one OTU. Therefore, we used the dada2 pipeline (Callahan et al. 2016) for these samples to analyse them after demultiplexing and adapter removal. Sequences were filtered and trimmed, then error rates were learned for sample inference. Finally, paired reads were merged and chimeras removed, resulting in the final amplicon sequence variant (ASV) table. The resulting 1184 ASVs were compared to nt database (last updated on 20.06.2022) of NCBI using BLAST+ 2.12.0 and 881 putatively non-dipteran ASVs were identified (i.e., 123 non-dipteran insects, 14 non-insect arthropods, 75 vertebrate, 3 non-arthropod invertebrates, 553 non-matching and 113 non-animals, mainly fungi, bacteria and plants) that accounted for 2.8% of the total reads in the samples. A small database consisting of only the Sanger sequences of the 16 specimens used for Mock A was generated in order to identify the ASVs with the assignTaxonomy function (bootstrap value 80). In this way, ASVs that resulted from minor variations in the sequences due to sequencing errors could be assigned to the corresponding specimen (OTU clustering). None of the putative non-dipteran ASVs were matched to one of the Sanger sequences. Of the putative dipteran ASVs, 254 ASVs were matched to one of the Sanger sequences, accounting for 99.96% of the reads, 49 ASVs were not matched and discarded. We compared the assigned ASVs with the respective Sanger sequence and discarded additional 21 ASVs (representing 0.03% of reads), which differed more than 1% (genetic distance) from the respective Sanger sequence. In this way, 1–27 ASVs were clustered into one of the 16 OTUs, each representing a specimen in Mock A. Finally, the mean numbers of reads (based on two technical replicates) for the 16 OTUs were calculated.

Data preparation

Mock B

Each of the 13 most abundant OTUs over all Mock B samples matched with one of the Sanger sequences of the 13 individuals (100% identity). In each sample, 96.7 to 99.8% of the reads were assigned to the 13 OTUs, the rest of the reads were mainly associated with spurious non-chironomid DNA. Except for one sample, where one of the replicates showed a contamination by another OTU (9.2% of reads) and therefore only 89.4% of the reads were assigned to the 13 OTUs. However, the number of reads per OTU was in the same range as in the other technical replicate, thus the replicate was not discarded. We calculated the mean number of reads (based on two technical replicates) for the 13 relevant OTUs; all other reads/OTUs were deleted.

Natural communities

The 160 OTUs resulting from bioinformatic processing were compared to nt database (last updated on 20.06.2022) of NCBI using BLAST+ 2.12.0 and 41 putatively non-dipteran OTUs were discarded (i.e., 7 non-dipteran insects, 7 non-insect arthropods, 1 vertebrate, 1 parazoa, 13 non-matching and 12 non-animals, mainly fungi, bacteria and plants). Negative controls showed no sign of cross-contamination among samples as the sum of reads per negative control was always low (< 112 reads, around 10 reads on average; all non-control samples contained > 1500 reads, around 22.700 on average). However, we detected 209 low-read false positives in the negative control samples (95.7% of them with less than 10 reads) presumably derived from tag switching. Additionally, read abundances per OTU were generally in the same order of magnitude in each of the two technical replicates per sample (i.e., less than 1 order of magnitude apart in 97.8% of the comparisons). In 95.9% of cases where only one of the two technical replicates contained reads for an OTU, the read numbers were spurious (i.e., < 10). Therefore, for artificial communities, read abundances were calculated as the average of both technical replicates only when both technical replicates contained reads. To further prevent false-positives, that can be introduced in any step of the metabarcoding procedure (for example spurious contamination or tag-switching, cf. Drake et al. 2021), we defined minimum sequence copy thresholds utilizing two technical replicates per sample and a relatively high number of negative controls. We set this threshold by identifying the percentage of average read abundance per OTU which allowed us to keep as many putative target reads (i.e., those with reads in both technical replicates) as possible while removing reads in the controls (false-positives) as far as possible. The threshold was set at 1% of the average read abundance per OTU in samples with a read abundance larger than zero in both technical replicates. Read abundances below this threshold, i.e., that were 100 times smaller than the average, were set to zero. In this way, 91.4% of reads in all negative controls were removed, while 97% of predominantly target reads (i.e., those with reads in both technical replicates) could be maintained. We assumed that OTUs with a read abundance larger than this specific threshold in both technical replicates indicated true presence of the species in a sample. Community compositions of natural communities were finally assessed based on presence-absence data and OTU sequences were taxonomically assigned to species level using the BOLD Identification System (IDS) for COI (Ratnasingham and Hebert 2007; last accessed on 14/02/2023) by querying against All Barcode Records on BOLD for animal identification. Five OTUs (out of 119) were not assigned to Chironomidae (3× Ceratopogonidae, 2× Chaoboridae) and therefore ignored in downstream analyses.

Statistical analysis

All statistical analyses were conducted in R version 4.2.0 (R Core Team 2019) using R Studio version 2022.07.1+554 (RStudio Team 2019). Plots were created using the packages dplyr version 1.1.0 (Wickham et al. 2021), ggplot2 version 3.4.1 (Wickham 2016), ggpmisc version 0.5.2 (Aphalo 2016) and ggpubr version 0.6.0 (Kassambara 2020). We used Hedges’ g of the package effsize version 0.8.1 (Torchiano 2016) as a measure of the effect size for comparison between the read abundance variation in different taxa and the variation within taxa (derived from the two different methods) in Experiment 1. For analysing the difference in read abundance of specimens in artificial communities when using equimolar amounts of purified total DNA vs. target DNA fragments (Experiment 2), assumptions for parametric hypothesis testing were checked using the functions shapiro.test and leveneTest of the package car version 3.1-1 (Fox and Weisberg 2019). As assumptions were not met, the read abundances of the individual OTUs were statistically compared using Kruskal-Wallis rank sum test (kruskal.test), while the deviation from mean read abundances in both groups (purified total DNA vs. target DNA fragments) was statistically analysed using Wilcoxon rank sum test (wilcox.test). Further, we calculated the Pearson’s correlation for linear regression analysis for analysing the effect on read abundance when using different amounts of purified total DNA and target DNA fragments (Experiment 2). For assessing the influence of the metabarcoding methods used for natural communities (Experiment 4), we conducted a Permutational Multivariate Analysis of Variance Using Distance Matrices (PERMANOVA) on the basis of Jaccard distance with 999 permutations using the function adonis2 and a correspondence analysis (CA) using the vegan package version 2.6-4 (Oksanen et al. 2019). To visualize the distribution of data, we utilized the ordiellipse function from the vegan package, which allowed us to draw ellipses based on the standard deviation of points.

Results

Bioinformatic analysis resulted on average in 22.299 (+/-4.068 SD) reads for Mock A samples (except one sample with 95.954 reads) and on average in 24.982 (+/-3.803 SD) reads for Mock B samples (except two samples with 54.942 and 38.809 reads, respectively). The natural community samples contained on average 163 (+/- 48 SD) reads per individual (except two samples with 6 and 31 reads per individual, resp., and one sample with 467 reads per individual).

Experiment 1: Effect of size-sorting and tissue preparation on read abundance of specimens in artificial communities

In total, 13 out of 16 specimens of Mock A and all specimens of Mock B were reliably detected in each replicate of every method. Three specimens of Mock A, corresponding to OTU_J, OTU_K and OTU_P, were detected in 97.5, 95 and 95% of the samples, respectively. For these three OTUs we found no reads in one of the two technical replicates in one or two samples. Variation in read abundance was higher between taxa (OTUs) within mock communities than within taxa comparing both metabarcoding approaches (Fig. 2B, D, F). The effect size Hedges’ g indicates large differences between the two groups (g > 0.8) in each case. Pre-adjusting the mass per specimen did not level off this variation (compare “within mocks” read abundance variation in Fig. 2B, D). The variation in individual read abundances exhibited substantial heterogeneity among the different taxa, which was observed in both analyses based on purified DNA and those based on dPCR (Fig. 2A, C, E).

Figure 2.

Effect of size-sorting and tissue preparation on individual read abundance. A, C, E show read abundances (logarithmic scale) per operational taxonomic unit (OTU) resulting from two different metabarcoding approaches (orange – purified DNA extracts, blue – direct PCR) with varying (2A) or equal (2C, E) mass per specimen from two different mock communities (Mock A, n = 5; Mock B, n = 3). Dots represent read abundances; means and ranges are indicated by vertical lines. B, D, F show boxplots illustrating variation in read abundance within mocks and within OTUs. Variation in read abundance is calculated as standard deviation of the compared values. The lower and upper hinges of boxplots correspond to the first and third quartiles. The upper whisker extends from the hinge to the largest value no further than 1.5 times the inter-quartile range from the hinge. The lower whisker extends from the hinge to the smallest value at most 1.5 times the inter-quartile range from the hinge. Data beyond the end of the whiskers, i.e., outliers, are plotted individually.

Experiment 2: Effect of DNA input variation on read abundance of specimens in artificial communities

When using the same amount of purified total DNA per specimen, read abundance varied significantly between the different OTUs (Kruskal-Wallis rank sum test, n = 4, p < 0.001; Fig. 3, orange). Equalizing the amount of DNA fragments after the first PCR led to more even read abundances, with no significant differences between specimens (Kruskal-Wallis rank sum test, n = 4, p = 0.42; Fig. 3, grey). Deviations from mean read abundance were significantly higher when using the same amount of purified DNA per specimen as compared to when equalizing the amount of DNA fragments after the first PCR (Wilcoxon rank sum test, n = 52, p < 0.001). While there was no linear relationship between the amount of purified total DNA used in the PCR and read abundance in any of the four mock communities (Pearson’s correlation, R2 = < 0.01, DF = 11, p = 0.83–0.93; Fig. 4, orange), read abundance showed a positive linear correlation with the amount of target fragment used in the second PCR in all four replicates (Pearson’s correlation, R2 = 0.74–0.94, DF = 11, p < 0.001; Fig. 4, grey).

Figure 3.

Individual read abundances resulting from equal amounts of purified total DNA or target DNA fragments. Read abundance of specimens in artificial communities (n = 4) after pooling equimolar amounts of purified total DNA (orange) or equal amounts of target DNA fragments (grey) per specimen. Dots represent read abundances, means and ranges are indicated by vertical lines. Overall read abundance per method is illustrated by boxplots on the right. The lower and upper hinges of boxplots correspond to the first and third quartiles. The upper whisker extends from the hinge to the largest value no further than 1.5 times the inter-quartile range from the hinge. The lower whisker extends from the hinge to the smallest value at most 1.5 times the inter-quartile range from the hinge. Data beyond the end of the whiskers, i.e., outliers, are plotted individually.

Figure 4.

Individual read abundances resulting from varying amounts of purified total DNA or target DNA fragments. Read abundance of specimens in artificial communities when using different amounts of purified DNA (orange colours represent four mock communities) or different amounts of target DNA fragments (grey, four mock communities) per specimen for community pools. DNA input units correspond to 0.05 ng of purified total DNA and 1 ng of target DNA fragments. Shown are regression lines and multiple R2 of the fitted linear regression models.

Experiment 3: Effect of DNA input dilutions on detection rates of specimens in artificial communities

When diluting the purified DNA extracts or tissue-water mixes prior to PCR, the reliability of OTU detection dropped in all approaches (Fig. 5). However, the decrease in OTU detection using purified DNA extracts is much lower than the decrease in detection rate based on tissue-water mixes, reflecting the faster deterioration of band intensity for dPCR samples visible during gel electrophoresis. While the total read abundance kept stable in all dilutions (between 13.414 and 58.062 reads) due to equimolar pooling of PCR products during library preparation, read abundances of individual OTUs increased or decreased with progressing dilution. OTUs with low initial read abundance (< 50 reads) failed to be detected in at least one of the dilutions, while OTUs with high initial read abundance (e.g., > 700 reads) were almost constantly detected throughout the dilutions.

Figure 5.

Effect of DNA input dilutions on detection rates in artificial communities. Mean number of detected OTUs (± range, n = 3) in a chironomid mock community with 16 individuals, using different metabarcoding methods and dilutions (logarithmic scale) of purified DNA extracts (orange) or tissue-water mixes (blue). We tested both the common and the direct PCR approach with (triangles, dark colours) and without (squares, light colours) prior size-sorting.

Experiment 4: Comparing standard versus dPCR metabarcoding approach for natural communities

Chironomid communities of eight different ponds in a mesocosm study were analysed using two different metabarcoding approaches, i.e., with and without DNA extraction. Ponds had been treated with Bti or left as a control and chironomid emergence was sampled on 5 sampling dates over a period of 17 days. In total, we detected 114 chironomid OTUs, corresponding to 36 different genera. Pairwise comparing each OTU in each sample, there were deviations between the two metabarcoding approaches in on average 8.3 (± 5.2 SD) percent of the OTUs per sample. In 6 out of 136 OTUs we saw a deviation between the two metabarcoding approaches in more than 4 samples (> 10%). These OTUs showed a relatively low mean read abundance (10 reads and less) compared to the overall mean read abundance (447 reads). Correspondence analysis showed high overlap in metabarcoding results from both approaches (Fig. 6). PERMANOVA indicated that a statistically significant effect on the community composition resulted from treatment (Bti vs. control) and sampling date, but not from the applied metabarcoding method (with or without DNA extraction; Table 1).

Table 1.

Results of the three-way PERMANOVA on natural communities. Assessed were the effects of Bti-treatment, sampling date, and DNA preparation method as well as their interactions on presence-absence chironomid community composition.

Source of variation DF sum of squares F statistics R² p
Bti-Treatment 1 0.870 3.730 0.046 0.001 ***
Sampling Date 4 2.534 2.715 0.135 0.001 ***
DNA Preparation Method 1 0.080 0.341 0.004 0.999
Bti-Treatment: Sampling Date 4 1.177 1.261 0.063 0.058
Bti-Treatment: DNA Preparation Method 1 0.004 0.018 < 0.001 1
Sampling Date: DNA Preparation Method 4 0.058 0.063 0.003 1
Bti-Treatment: Sampling Date: DNA Preparation Method 4 0.007 0.007 < 0.001 1
Residuals 60 14.001 0.747
Total 79 18.731
Figure 6.

Correspondence analysis (CA) plot comparing standard versus dPCR metabarcoding approach for natural communities. Indicated are differences of chironomid communities from four Bti-treated (triangles, dark colours) and four control (squares, light colours) ponds of a mesocosm study over five sampling dates (N = 40). Two different metabarcoding approaches, i.e., with (orange) or without (blue) DNA extraction, were used for each community. Proportion explained per axis: CA1–7.8%, CA2–7.5%. Ellipsoids indicate standard deviations of points. Two very distant points are not visible in this figure.

Discussion

In the current study, we showed for the first time that direct PCR protocols can be combined with insect community metabarcoding approaches. Metabarcoding of chironomid communities with and without DNA extraction, produced highly comparable results. Both methods led to similar detection rates in artificially created communities with known chironomid composition. Furthermore, samples of natural communities from an ongoing mesocosm study (Kolbenschlag et al. 2023), which had been processed with and without DNA extraction, showed to be highly similar in their species composition, allowing identical conclusions about the community response to an environmental stressor.

We observed high similarity between the presence-absence chironomid community composition detected by metabarcoding with and without DNA extraction. We explain the faster decrease of specimen detection rates in sequential dilutions of the dPCR approach as compared to the common approach with the stochastic effect that arises, when very small amounts of input material are pipetted (i.e., 1 µl of purified DNA extract or 5 µl of tissue-water mix). Stochasticity probably also explains the observed higher read abundance variation in the dPCR samples as compared to the common approach. As DNA is purified and concentrated during DNA extraction, the probability to represent the whole community in a small amount of purified DNA extract is much higher than in a similar amount of tissue-water mix. This lower representativeness could indeed raise concern about the reliability of direct PCR metabarcoding when studying natural communities that are typically composed of few abundant and many rare species (e.g., Fisher et al. 1943; Hubbell 2001). However, we found only minor differences between the dPCR and DNA extraction metabarcoding approach when studying natural chironomid communities (Fig. 6). In addition, the differences did not alter the conclusions about the outcome of the ecotoxicological case study. Therefore, we conclude that direct PCR metabarcoding is an equally valid tool to characterize chironomid community composition as compared to the common metabarcoding approach that includes DNA extraction.

While detection rates of the two approaches were largely similar, read abundance was influenced by the method used (Fig. 2). The individual read abundance was lower in direct PCR approaches for most of the specimens, while the variation between specimens was several factors higher. Elbrecht and Leese (2015) argued that the different body mass of specimens could have an effect on read abundance and recommended size-sorting of individuals prior to analysis. As our specimens varied in their individual mass up to a factor of more than 20, we expected an alignment of individual read abundances after equalizing the available tissue or DNA per specimen. Surprisingly, even though this treatment substantially improved the detection rates of very small specimens, individual read abundances still varied greatly between specimens (Fig. 2).

Investigating the reasons for differential individual read abundances, we found out that adjusting the availability of purified DNA prior to PCR did not eliminate variation. At the same time, varying the availably of DNA target fragments after PCR had a much larger influence on read abundance (Figs 3, 4). We therefore conclude that most of the variation arose during PCR probably due to primer bias, even though our specimens are closely related taxa and primers were highly degenerated. If and to what extent interspecific mitochondrial copy number variation (as described by Stefano et al. 2017) plays a role here is beyond the scope of this study and further studies are required to evaluate the potential impact of heteroplasmy and nuclear pseudogenes of mitochondrial origin on read abundance in insect metabarcoding studies (cf. D’Errico et al. 2004 and Guo et al. 2022). Given the weak relationship between individual tissue mass and read abundance, we agree with previous studies (e.g., Amend et al. 2010, Sun et al. 2015) that metabarcoding of arthropod communities should be evaluated based on presence-absence data, rather than assuming read abundance would reflect biomass proportions.

We tested the comparability of the two metabarcoding approaches using specimens from the family Chironomidae. However, dPCR protocols have been developed and successfully applied to several other invertebrate taxa, e.g., mosquitoes: Werblow et al. (2016), spider mites: Sakamoto and Gotoh (2017), ants: Wang et al. (2018). Wong et al. (2014) tested the suitability of dPCR protocols for a range of different invertebrates. While they were able to successfully optimize dPCR procedures for different flies and sea stars, the authors reported low success rates for beetles, copepods, ants and odonats. Wong et al. (2014) assumed that their dPCR protocol was less successful for heavily sclerotized taxa and/or those who carry a lot of PCR inhibitors (e.g., melanin, haemocyanin or secretions from exocrine glands). In contrast, the here presented dPCR approach, using finely ground tissue in water, produced COI sequences from specimens of the arthropod orders Araneae (Pirata and Tetragnatha), Coleoptera (Dasytes and Neocrepidodera), Diptera (Bezzia and Helina), Ephemeroptera (Caenis and Cloeon), Hemiptera (Pyrrhocoris), Hymenoptera (Chelonus, Formica and Lasius), Lepidoptera (Eudonia) and Odonata (Coenagrion). We assume that our thorough mechanical tissue-breakdown step (cf., Elbrecht and Leese 2015; Buchner et al. 2021), supported by the dilution of inhibitors when tissue is mixed with water, resulted in a higher PCR success for several difficult taxa. Additionally, longevity of samples could be easily improved by using buffer instead of water to dissolve ground tissue (see Anchordoquy and Molina 2007). Our results are in accordance with Thongjued et al. (2019) who concluded that a thorough cell lysis, dilution of inhibitors and the use of an inhibitor tolerant DNA polymerase contribute to successful direct PCR even in challenging taxa. Although our initial positive results on Sanger sequencing of dPCR products of various arthropods suggest a successful extension to metabarcoding approaches, applicability of direct PCR metabarcoding to other arthropod groups needs to be verified. As shown for chironomids, representative mock communities need to be established to ensure reliable detection of relevant taxa in metabarcoding of mixed arthropod samples and to rule out dPCR-related biases between taxa. Further, upscaling the approach for very large sample sizes, such as from Malaise traps, could be tested. Given these initial tests are successful, dPCR metabarcoding might become available in mixed arthropod community samples which do not require laborious morphological sorting and therefore contribute to time and cost-efficient monitoring of natural, environmentally stressed or experimental communities.

In conclusion, we showed that DNA extraction can be omitted in chironomid community metabarcoding while preserving the informative value of the presence-absence community composition. As the adoption of our proposed direct PCR protocols is relatively easy and inexpensive, direct PCR metabarcoding has the potential to become a standard procedure in chironomid community analysis. With modern PCR reagents being more robust to contamination by inhibitors, we assume that direct PCR metabarcoding, where DNA extraction is replaced by solely a mechanical tissue-breakdown step, is applicable for a wide range of arthropod taxa and encourage further comparative studies. The opportunity to avoid DNA extraction steps might even aid the ongoing development to miniaturize the instrumental requirements for PCR and metabarcoding (e.g., for application in the field or for live on-site monitoring). In addition, laboratories implement more and more automation processes, in order to reduce the costs for high-throughput application of metabarcoding. Eliminating the DNA extraction step entirely might contribute a substantial improvement in this development. In the end, faster and cheaper metabarcoding procedures will boost our ability to monitor the diversity of arthropod communities and appropriately target conservation actions to combat the ongoing global biodiversity loss.

Acknowledgements

The authors would like to cordially thank Sara Kolbenschlag and the Eußerthal Ecosystem Research Station Team for providing preserved chironomids. We are grateful to Sophie Stoll for fabulous work in the laboratory and to our many student helpers for size-sorting insects. We would like to thank Carola Greve, Damian Baranski und Alexander Ben Hamadou of TBG for laboratory support during the measurements of target fragment concentrations. Thanks to Pavel Bystřický for help in the lab during initial tests of the direct PCR approach. Special thanks to Thomas Mehner for constructive advice on earlier versions of the manuscript and continuous support throughout the writing process. We further thank Verena Gerstle, Sara Kolbenschlag, Sebastian Pietz, Alexis Roodt and Caroline Ganglo for a stimulating atmosphere while discussing scientific writing. Thank you very much to Dominik Buchner for his effort in providing a detailed and constructive review.

Additional information

Conflict of interest

The authors declare they have no conflict of interest.

Ethical statement

No ethical statement was reported.

Funding

This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Research Training Group SystemLink 326210499/GRK2360.

Author contributions

NR and KS conceived the ideas and designed methodology; NR collected the data; NR and KS analysed the data; NR led the writing of the manuscript. Both authors contributed critically to the drafts and gave final approval for publication.

Author ORCIDs

Nina Röder https://orcid.org/0000-0002-9681-7538

Klaus Schwenk https://orcid.org/0000-0003-2427-4332

Data availability

Raw sequences were deposited in GenBank SRA archive and are available with the BioProject accession number PRJNA989176. Data and R scripts are available from Zenodo https://doi.org/10.5281/zenodo.8074454 (Röder and Schwenk 2023).

References

  • Aljanabi S (1997) Universal and rapid salt-extraction of high quality genomic DNA for PCR- based techniques. Nucleic Acids Research 25(22): 4692–4693. https://doi.org/10.1093/nar/25.22.4692
  • Baloğlu B, Clews E, Meier R (2018) NGS barcoding reveals high resistance of a hyperdiverse chironomid (Diptera) swamp fauna against invasion from adjacent freshwater reservoirs. Frontiers in Zoology 15(1): 1–31. https://doi.org/10.1186/s12983-018-0276-7
  • Buchner D, Haase P, Leese F (2021) Wet grinding of invertebrate bulk samples – a scalable and cost-efficient protocol for metabarcoding and metagenomics. Metabarcoding and Metagenomics 5: e67533. https://doi.org/10.3897/mbmg.5.67533
  • Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP (2016) DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods 13(7): 581–583. https://doi.org/10.1038/nmeth.3869
  • Creer S, Deiner K, Frey S, Porazinska D, Taberlet P, Thomas WK, Potter C, Bik HM (2016) The ecologist’s field guide to sequence‐based identification of biodiversity. Methods in Ecology and Evolution 7(9): 1008–1018. https://doi.org/10.1111/2041-210X.12574
  • D’Errico I, Gadaleta G, Saccone C (2004) Pseudogenes in metazoa: Origin and features. Briefings in Functional Genomics & Proteomics 3(2): 157–167. https://doi.org/10.1093/bfgp/3.2.157
  • Drake LE, Cuff JP, Young RE, Marchbank A, Chadwick EA, Symondson WOC (2021) An assessment of minimum sequence copy thresholds for identifying and reducing the prevalence of artefacts in dietary metabarcoding data. Methods in Ecology and Evolution 13(3): 694–710. https://doi.org/10.1111/2041-210X.13780
  • Elbrecht V, Leese F (2015) Can DNA-based ecosystem assessments quantify species abundance? Testing primer bias and biomass—sequence relationships with an innovative metabarcoding protocol. PLoS ONE 10(7): e0130324. https://doi.org/10.1371/journal.pone.0130324
  • Elbrecht V, Leese F (2017) Validation and Development of COI Metabarcoding Primers for Freshwater Macroinvertebrate Bioassessment. Frontiers in Environmental Science 5. https://doi.org/10.3389/fenvs.2017.00011
  • Elbrecht V, Steinke D (2018) Scaling up DNA metabarcoding for freshwater macrozoobenthos monitoring. Freshwater Biology 64: 380–387. https://doi.org/10.1111/fwb.13220
  • Fisher RA, Corbet AS, Williams CB (1943) The relation between the number of species and the number of individuals in a random sample of an animal population. Journal of Animal Ecology 12(1): 42–58. https://doi.org/10.2307/1411
  • Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R (1994) DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology 3: 294–299. https://pubmed.ncbi.nlm.nih.gov/7881515/
  • Gostel MR, Kress WJ (2022) The Expanding role of DNA barcodes: Indispensable tools for ecology, evolution, and conservation. Diversity (Basel) 14(3): e213. https://doi.org/10.3390/d14030213
  • Hubbell SP (2001) The unified neutral theory of biodiversity and biogeography. Princeton University Press.
  • Kolbenschlag S, Gerstle V, Eberhardt J, Bollinger E, Schulz R, Brühl CA, Bundschuh M (2023) A temporal perspective on aquatic subsidy: Bti affects emergence of Chironomidae. Ecotoxicology and Environmental Safety 250: e114503. https://doi.org/10.1016/j.ecoenv.2023.114503
  • Lewis RL, Beckenbach AT, Mooers AØ (2005) The phylogeny of the subgroups within the melanogaster species group: Likelihood tests on COI and COII sequences and a Bayesian estimate of phylogeny. Molecular Phylogenetics and Evolution 37(1): 15–24. https://doi.org/10.1016/j.ympev.2005.02.018
  • Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Szoecs E, Wagner H (2019) vegan: Community Ecology Package. https://CRAN.R-project.org/package=vegan
  • R Core Team (2019) R: A language and environment for statistical computing (version 3.6.2). R Foundation for Statistical Computing, Vienna. https://www.R-project.org/
  • Röder N, Schwenk K (2023) Data from: Direct PCR meets high-throughput sequencing – metabarcoding of chironomid communities without DNA extraction, Zenodo, Dataset. https://doi.org/10.5281/zenodo.8074454
  • RStudio Team (2019) RStudio: Integrated Development Environment for R. RStudio Inc, Boston, MA.http://www.rstudio.com/
  • Sabo JL, Bastow JL, Power ME (2002) Length-mass relationships for adult aquatic and terrestrial invertebrates in a California watershed. Journal of the North American Benthological Society 21(2): 336–343. https://doi.org/10.2307/1468420
  • Sakamoto H, Gotoh T (2017) Non-destructive direct polymerase chain reaction (direct PCR) greatly facilitates molecular identification of spider mites (Acari: Tetranychidae). Applied Entomology and Zoology 52(4): 661–665. https://doi.org/10.1007/s13355-017-0512-1
  • Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P (1994) Evolution, Weighting, and Phylogenetic Utility of Mitochondrial Gene Sequences and a Compilation of Conserved Polymerase Chain Reaction Primers. Annals of the Entomological Society of America 87(6): 651–701. https://doi.org/10.1093/aesa/87.6.651
  • Sun C, Zhao Y, Li H, Dong Y, MacIsaac HJ, Zhan A (2015) Unreliable quantitation of species abundance based on high-throughput sequencing data of zooplankton communities. Aquatic Biology 24(1): 9–15. https://doi.org/10.3354/ab00629
  • Thongjued K, Chotigeat W, Bumrungsri S, Thanakiatkrai P, Kitpipit T (2019) A new cost-effective and fast direct PCR protocol for insects based on PBS buffer. Molecular Ecology Resources 19(3): 691–701. https://doi.org/10.1111/1755-0998.13005
  • Wang WY, Srivathsan A, Foo M, Yamane SK, Meier R (2018) Sorting specimen-rich invertebrate samples with cost-effective NGS barcodes: Validating a reverse workflow for specimen processing. Molecular Ecology Resources 18(3): 490–501. https://doi.org/10.1111/1755-0998.12751
  • Werblow A, Flechl E, Klimpel S, Zittra C, Lebl K, Kieser K, Laciny A, Silbermayr K, Melaun C, Fuehrer H-P (2016) Direct PCR of indigenous and invasive mosquito species: A time- and cost-effective technique of mosquito barcoding. Medical and Veterinary Entomology 30(1): 8–13. https://doi.org/10.1111/mve.12154
  • Wong WH, Tay YC, Puniamoorthy J, Balke M, Cranston PS, Meier R (2014) ‘Direct PCR’ optimization yields a rapid, cost-effective, nondestructive and efficient method for obtaining DNA barcodes without DNA extraction. Molecular Ecology Resources 14(6): 1271–1280. https://doi.org/10.1111/1755-0998.12275
  • Yu DW, Ji Y, Emerson BC, Wang X, Ye C, Yang C, Ding Z (2012) Biodiversity soup: Metabarcoding of arthropods for rapid biodiversity assessment and biomonitoring. Methods in Ecology and Evolution 3(4): 613–623. https://doi.org/10.1111/j.2041-210X.2012.00198.x

Supplementary materials

Supplementary material 1 

Overview of chironomid size classes

Nina Röder, Klaus Schwenk

Data type: additional methodological information

Explanation note: Chironomid individuals were sorted into four different size categories according to their body length (from the anterior margin of the head between the antennae to the end of the posterior abdominal segment) and shape (thin – usually males or thick – usually females). They were counted and average weight per specimen was determined.

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (18.71 kb)
Supplementary material 2 

Composition of the two artificial chironomid communities

Nina Röder, Klaus Schwenk

Data type: additional methodological information

Explanation note: Sanger sequences of specimens were compared to BOLD database.

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (26.54 kb)
login to comment