A total of 23 terrestrial arthropod species (including 21 insects, a collembolan and a crustacean) were obtained as live specimens from either standardised cultures of commercial suppliers, donations by other laboratories and the v (NRM, Stockholm) vertebrate collection, or collected from the surroundings of the NRM (Suppl. material 7: Table S1). All specimens were killed by submersion in 99% ethanol. Ten types of mock communities (A-J), with four replicates each (a total of 40 tubes), were assembled in 50 mL Falcon tubes (see Suppl. material 1: “Community types” for details). Each community was composed of 22 species, so each of them had a species missing that was present in the rest. In other words, our study explicitly addresses the effect of the extraction protocol on the detectability of species against a community background of standardised complexity.

The number of individuals of each species was always the same regardless of the community type. For example, Drosophila melanogaster was represented by six individuals in all ten community types, while D. yakuba was represented by three individuals in all community types, except for Community H from which D. yakuba was excluded). The total number of insects per community ranged from 70 to 74. Since our specific interest was in the impact of species properties on the detectability of species, rather than of individual variation in body size, the average weight per specimen was computed by recording the dry weight of ten individuals per species. We then selected all individuals of a species from the cultures to be of the same size. As a result of this rationale, the study is well aimed to detect effects of species averages, whereas it provides little information on the effects of added variation in individual size. The proportions of each species in the communities in terms of weight and numbers was registered as abundance in biomass or number, respectively.

Reference DNA barcode sequences were constructed using one additional individual of each species and one individual of the bumblebee Bombus pascuorum (used as a control for quantifying index swapping during the sequencing run, see below) as follows. DNA was extracted from all individuals using the KingFisher Cell and Tissue DNA kit on a KingFisher Duo robot (Thermo Fisher Scientific), except for Encarsia formosa, Folsomia candida and Tuberculatus annulatus, which were processed using QIAamp DNA Micro kit (Qiagen) following the manufacturer’s protocol. The entire barcoding region of COI (658 bp) and a fragment of 450–490 bp (depending on the species) close to the 5’ end of the mitochondrial 16S rRNA gene were amplified and sequenced. COI was amplified using the primers jgLCO1490-jgHCO2198 (Geller et al. 2013), with the exception of Formica rufa, which was amplified with the primers LepF1-LepR1 (Hebert et al. 2004) following failed amplification attempts with the previous pair. 16S was amplified with the primers 16Sar-16Sb2 (Simon et al. 1994; Cognato and Vogler 2001) in all species. The PCR mix consisted of one Illustra Hot Start Taq RTG bead (GE Healthcare Life Sciences), 1 µL (10 pmoles) of each primer, 2 µL of DNA template and 21 µL of biology-grade water (final volume: 25 µL); the temperature protocol consisted of an initial phase of denaturation at 95 °C for 5 min followed by 40 cycles of 30 s of denaturation at 95 °C, 45 s of annealing at 50 °C (COI-jg) / 45 °C (COI-Lep) / 48 °C (16S), and 60 s of extension at 68 °C, finishing with a final phase of extension of 10 min at 72 °C. PCR success was checked in an agarose gel and those reactions with positive bands were cleaned of single-strand DNA molecules with ExoSAP-IT (Thermo Fisher Scientific) following the manufacturer’s protocol and sent to Macrogen Europe B. V. (Amsterdam, Netherlands) for two-strand Sanger sequencing. The sequences were merged, edited and trimmed of primers using GENEIOUS v8.1.7 (Kearse et al. 2012). The reference sequences of both genes of all 24 species are given in the Suppl. material 2, 3: “COI reference library” and “16S reference library”, respectively.

The mock communities were subjected to four different lysis treatments, resulting from the combination of two digestion buffers (referred to as B1 or Gentle, and B2 or Aggressive) and two incubation times (Fig. 1). Buffer B1 was taken from Vesterinen et al. (2016) (which was, in turn, modified from Aljanabi and Martinez (1997)) and consisted of 400 mM NaCl, 10 mM Tris-HCl pH 8.0, 2 mM EDTA pH 8.0, 2% SDS, 0.1% proteinase K 20 mg/mL (which corresponds to 0.069 mM) and molecular biology grade water. Buffer B2 was modified from Gilbert et al. (2007) (see also Ji et al. 2020) and consisted of 115 mM NaCl, 47 mM Tris-HCl pH 8.0, 3 mM CaCl₂, 2.4% SDS, 40 mM Dithiothreitol (DTT), 1% proteinase K 20 mg/mL (which corresponds to 0.69 mM) and molecular biology grade water (see Suppl. material 8: “Lysis buffers and purification protocols” for details).

Figure 1.

Schematic overview of the experimental design and treatments. Each community (A–J) was represented by four initial replicates (four tubes containing the same mix of species, each of them with the same number of specimens). Each replicate was incubated in one buffer (B1 or B2) for 2 h 30’ (LT1) or 5 h (LT2). Then, DNA from each lysate (B1–LT1, B1–LT2, etc.) was extracted by using both a manual salt saturation-salt protocol (P1) or silica-coated magnetic beads in a robot (P2). Thus, for each original community (A–J), eight DNA extracts were obtained, one from each combination of purification method, buffer and lysis time.

Of the reagents used, sodium dodecyl sulphate (SDS) is a surfactant that breaks cells by disrupting the membrane, and both proteinase K and DTT have proteolytic activity. Thus, having a higher concentration of these compounds, buffer B2 is expected to produce a more aggressive digestion than buffer B1. In addition, buffer B1 contains EDTA, which inactivates proteolytic enzymes, such as nucleases and proteases, while buffer B2 contains Ca²⁺, which activates these enzymes.

In terms of incubation times, samples were split into two different treatments: 2 hours and 30 minutes vs. 5 hours (referred to as LT1 and LT2). When combined with the different lysis buffers, this resulted in four different treatments, each representing a unique combination of buffer and lysis time (B1–LT1, B1–LT2, B2–LT1, B2–LT2). In each case, the samples were incubated in 20 mL of buffer at 56 °C with a slight agitation in an orbital shaker, for each community replicate. Once the incubation time was over, the lysate was decanted out and collected for DNA extraction. The insects were rinsed with molecular biology-grade water, then with clean 70% ethanol, and finally stored in 80% ethanol. The insects remained at all times inside the tubes.

DNA from each lysate was extracted using two purification methods (referred to as P1 and P2). In short, the protocols differed as follows (see Suppl. material 8: “Lysis buffers and purification protocols” for full details). Purification P1 followed a manual salt saturation protocol, modified from Vesterinen et al. (2016) and Aljanabi and Martinez (1997). Taking 7.5 mL of lysate as input material, the proteins and cell membranes were first precipitated using a saturated salt solution, and the DNA was then precipitated using isopropanol, both precipitation steps being aided by centrifugation, with an elution volume of 150 µL. Purification P2 was conducted using silica-coated magnetic beads from the KingFisher Cell and Tissue DNA kit on a KingFisher Duo robot, following the manufacturer’s protocol, with an input volume of 225 µL of lysate and an elution volume of 150 µL as well. Summarising, for each community, eight DNA extracts were obtained, making a total of 80 DNA samples (10 mock communities × 2 lysis buffers × 2 incubation times × 2 purification methods). The DNA extracts were measured for concentration and purity (ratio of absorbance at 260 and 280 nm, henceforth referred to as the A260/A280 ratio) using a NanoVue instrument (version 4282 v1.7.3, GE Healthcare Life Sciences). A260/A280 ratios between 1.8 and 2 are considered optimal.

A 321 bp fragment of COI was amplified with a modified version of the primer pair BF2-BR1 (Elbrecht and Leese 2017a, modified in Marquina et al. 2019a) and a 345 bp fragment of mitochondrial 16S rRNA was amplified with the primer pair Chiar16SF-Chiar16SR (Marquina et al. 2019a). A two-step PCR protocol was followed for library preparation, with a marker-specific first-round PCR with primers with Illumina overhangs attached at the 5’ end and an indexing second-round PCR with indexed Illumina adapters. The PCR mix consisted of one Illustra Hot Start Taq RTG bead (GE Healthcare Life Sciences), 10 pmoles of each primer, 2 µL of DNA template and 21 µL of biology‐grade water (final volume: 25 µL). The temperature protocol of the two rounds of PCR consisted of an initial phase of denaturation at 95 °C for 5 min followed by 25 cycles of 30 s of denaturation at 95 °C, 45 s of annealing at 48 °C (COI) / 50 °C (16S) and 45 s of extension at 68 °C (first round), or 15 cycles of 30 s of denaturation at 95 °C, 30 s of annealing at 62 °C and 60 s of extension at 72 °C (second round), and a final extension phase of 10 min at 72 °C. Four libraries containing only Bombus pascuorum amplicons were added to control for index swapping, and extraction and PCR blanks to control for other sources of error. Without separate index swapping controls, it would be impossible to differentiate between “stray” sequences resulting from tag-switching amongst sequences emanating from the original samples and those resulting from contamination. PCR products were quantified with a Qubit Fluorometer (Thermo Fisher Scientific), pooled in equimolar concentration and purified with QIAquick gel extraction kit (Qiagen) after cutting the bands of the desired length from an agarose gel. They were sequenced on an Illumina MiSeq using v3 chemistry and a 2 × 300 bp paired‐end run at SciLifeLab facilities (Stockholm).

The detailed bioinformatic pipeline with commands and options can be found in Suppl. material 4: “Bioinformatic pipeline”. Sequences were processed with the OBITools pipeline (Boyer et al. 2016), complemented by other programmes and scripts (mainly from https://github.com/metabarpark/R_scripts_metabarpark). An attribute containing sample information was added to all sequence headers before pooling together all the libraries in a single file (COI and 16S separately). Reads were quality checked with FastQC (Andrews 2010), the ends trimmed when average quality dropped below a Phred score of 28, and finally paired-end-merged, discarding sequences with an alignment score lower than 40. Primers were trimmed away using CUTADAPT v1.8.0 (Martin 2011), and only the reads with the desired length were kept for downstream analysis (310–330 bp for COI, 290–370 bp for 16S).

Subsequently, the reads were dereplicated and chimeras were filtered out using the uchime_denovo function in VSEARCH v2.7.1 (Rognes et al. 2016) and clustered into Molecular Operational Taxonomic Units (MOTUs) using SWARM v2.1.13 (Mahé et al. 2015). The maximum distance d allowed during clustering was 13 for COI and 6 for 16S, which corresponds to a sequence divergence of 3–4% and 1–2% for COI and 16S, respectively. MOTU occurrence tables were curated with LULU v0.1.0 (Frøslev et al. 2017) to collapse NUMT-derived erroneously generated MOTUs into their parent MOTUs. The centroid sequences of every MOTU generated by SWARM were compared against the previously constructed reference libraries using the ecotag script from OBITools, and the resulting file was merged with the abundance table. The resulting file was curated with the refine_MOTU_table script (https://github.com/metagusano/metabarcoding_scripts) to remove MOTUs with a relative abundance per sample lower or equal to that generated by index swapping (i.e. abundance of reads assigned to Bombus pascuorum in the mock community samples), as well as with less than 10 reads in total, and to collapse MOTUs with coincident species identification.

All data analyses were performed in R v.3.3.3 (R Core Team 2017). We first visualised the differences in the estimated community composition produced by different methods, using non-metric multidimensional scaling (NMDS), based on Bray-Curtis dissimilarity (functions vegdist and metaMDS from package ‘vegan’ (Oksanen et al. 2013)). To then pinpoint the effects of individual methodological choices, we ran a series of tests. First, we tested whether differences in digestion buffer and incubation time during lysis, and in purification method post-lysis, had effects on the concentration and purity of the DNA extract. Specifically, we applied a split-split-plot analysis of variance (ANOVA) to values of DNA concentration as a function of Buffer (main plot), Incubation time (subplot) and Purification (sub-subplot), with Community type (A-J) as replicate (function ssp.plot from package ‘agricolae’ (de Mendiburu 2019)). For the purity of the DNA extract, another split-split-plot ANOVA was fitted to the A260/A280 values with the same formula as for the concentration. Subsequently, a post-hoc least significant difference test (function LSD.test from package ‘agricolae’) was used to conduct pairwise comparisons between groups. Then, to investigate if the differences in concentration and purity of the DNA extract had any effects on the recovery of species after sequencing, we repeated the split-split-plot analysis with the number of species recovered as the response variable. This analysis was performed for each marker separately.

As the experimental set-up did not allow for an analysis of correlation of real abundance to read abundances, we then investigated how the lysis buffers and the digestion times performed in recovering compositional-related metrics from the samples. We compared the differential in alpha diversity (Shannon index, H’) from the original mock communities to estimates obtained through metabarcoding (function diversity from package ‘vegan’). We also calculated the Kullback-Leibler Divergence between the observed community and the original, known composition of the mock community (function KLD from package ‘LaplacesDemon’ (Statisticat 2018)). In brief, the Kullback-Leibler Divergence takes a distribution as its reference (in this case, the true proportions of the species in the mock communities) and calculates how much information must be added to a second distribution (the proportion of reads of each species in the metabarcoding sample) to make it equal to the reference. To model the impact of methodological choices on the Shannon Index and the Kullback-Leibler Divergence between the true and observed sample, we applied another split-split-plot analysis with the values as a function of the three factors, as previously done for purity and concentration. The purification term was not significant in any of the split-split-plot analyses of the Shannon Index and the Kullback-Leibler Divergence and was thus dropped. These analyses were done separately for COI and 16S and using both reference distributions based on relative abundance and based on relative biomass, and with both datasets rarefied to the number of reads of the sample with the lowest number of reads.

The concentration of the DNA extracts ranged from 4.6 to 371.5 ng/µL (Suppl. material 7: Table S2). The lowest mean values corresponded to the replicates incubated in buffer B1 (for 2 h 30’ or 5 h) and subsequently purified using the extraction robot. By comparison, the highest mean values corresponded to the replicates incubated in buffer B2 for 5 hours, then purified with the manual salt saturation method. The ratio A260/A280 ranged from 1.05 to 2.00. The lowest mean values corresponded to the replicates incubated in buffer B1 for 2 h 30’ and purified with the robot and those incubated in buffer B2 also for 2 h 30’ and purified manually. The highest mean values corresponded to the replicates incubated in buffer B1 for 5 hours and purified manually.

All three factors and their interaction had significant effects on the DNA concentration of the extracts (Suppl. material 7: Table S3). Lysis with buffer B2 increased the DNA concentration, same as longer (LT2) incubation (Fig. 3, upper panel). In addition, extracting the DNA from the lysate with the salt saturation method generated extracts with a much higher concentration than those extracted with the automated robot extraction using silica beads. Undoubtedly, this was partly influenced by the starting volume being 7.5 mL in the first case and only 0.225 mL in the latter.

Several factors had a significant effect on the purity of the extract (ratio A260/A280; Suppl. material 7: Table S4). However, neither Buffer nor the three-way interaction showed significant effects. When buffer B1 was used for lysis, the purification method based on salt saturation produced extracts of significantly higher purity than did the method based on silica beads, but the incubation time did not have a strong effect. In contrast, when using buffer B2, LT2 improved the purity of the extract (Fig. 3, lower panel).

Figure 3.

Concentration and purity of the DNA extracts from different extraction methods. DNA concentration (upper panel) clearly increases with buffer aggressiveness and incubation time using the manual salt saturation purification protocol, while the increase due to incubation time is less clear, but the effect of lysis buffer can still be appreciated when using the robot purification protocol. Note that the starting input volume of lysate is 7.5 mL for the manual purification method and 225 µL (30 times smaller approximately) for the robot, while elution volume is 150 µL in both cases. Purity of the DNA extract (lower panel) is higher for the manual purification and the longer incubation times, regardless of the lysis buffer.

These differences had no significant effect on the number of species recovered with 16S. However, for COI, the Buffer and Purification effects were both significant, albeit small. Specifically, the number of species recovered was slightly higher for buffer B2 and for the salt saturation protocol. For COI, the buffer affected the mean number of species recovered (Suppl. material 7: Table S5). Against this background, the lysis time did not modify the mean (Suppl. material 7: Table S3), but the purification method did (Suppl. material 7: Table S3). Given significant two-way interactions between Buffer and Purification, we arrive at a scenario where combination B2–P1 provides the highest response in terms of species recovered, with B2–P2 second, then B1–P1 and B1–P2 the lowest (Suppl. material 7: Fig. S1). For 16S, there was no significant effect of any of the factors in the number of species recovered nor of their interactions (Suppl. material 7: Table S6). The average number of species recovered ranged between 20.5–21.5 for COI and 18.2–19 for 16S.

In terms of community composition, the metabarcoding results did not resemble the mock communities in neither specimen abundance nor biomass and they all clustered together in the NMDS plot (Suppl. material 7: Figs S2–S4). Although more similar amongst themselves than to the mock communities, replicates treated with different lysis buffers and incubation times showed great variability. Purification replicates, on the other hand, were very similar to each other (see bars distribution in Suppl. material 7: Fig. S3,S4). A few large species dominated most samples, namely Calliphora vomitoria (Diptera), Acheta domesticus (Orthoptera), and Locusta migratoria (Orthoptera), but A. domesticus was over-represented in the metabarcoding datasets compared to the real abundance in terms of biomass. Fairly small insects like Aphidius colemani, A. ervi, Encarsia formosa (Hymenoptera), Aphidoletes aphidimyza (Diptera) and Folsomia candida (Collembola), were more abundant with respect to read numbers for both markers than in specimen number or biomass in the communities, as well as Drosophila simulans and D. yakuba (Diptera). On the other hand, larger species like Blaptica dubia (Blattodea), Dermestes haemorrhoidalis (Coleoptera), and Polygonia c-album (Lepidoptera), were also under-represented in the dataset. For other species, the bias depended on the marker. For instance, Macrolophus pygmaeus (Hemiptera) was over-represented in the 16S dataset, as well as some Drosophila species (Diptera) in the COI dataset, while Cryptolaemus montrouzieri (Coleoptera) was usually under-represented in the 16S dataset (Fig. 4, and Suppl. material 7: Figs S2, S3).

Figure 4.

Representation of sequencing reads relative to biomass per sample. Relative representation is calculated as the log-ratio between relative read abundance and relative abundance in biomass of each species in each replicate. A higher log-ratio indicates that the species is over-represented in the metabarcoding dataset, while a lower value indicates that the species is under-represented in relation to its relative abundance in biomass in the mock community. Each community’s replicates are indicated in the following order: B1.LT1.P1, B1.LT1.P2, B1.LT2.P1, B1.LT2.P2, B2.LT1.P1, B2.LT1.P2, B2.LT2.P1, and B2.LT2.P2. No reads were recovered for COI from sample C-B2.LT2.P1.

Regarding alpha diversity, replicates incubated in buffer B1 returned values of the Shannon Index (H’) with a smaller decrease compared to those of the mock community based both on biomass and specimen numbers, irrespective of whether they had been incubated for LT1 or LT2 (Fig. 5, Suppl. material 7: Tables S7–S10). LT2 produced samples with a slightly lower value of H’ in the 16S dataset, but the differences were not significant. In contrast, lysis in buffer B2 generated metabarcoding community estimates with a significantly lower diversity than the mock communities they originated from. This was true regardless of marker or lysis time.

Figure 5.

Estimated decrease in alpha diversity, measured as Shannon Index (H’), for different incubation treatments and markers. A short and gentle lysis (B1, LT1) recovers diversity values closer to the actual values of the original sample measured in biomass with the 16S marker (blue) and an increase in lysis time and chemical aggressiveness of the buffer returns more distant values. With COI (orange), this effect is dependent only on the lysis buffer. The decrease in alpha diversity compared to the mock samples based on number of individuals is greater than based on biomass, but they reproduce the same pattern. The black line indicates H’(mock community) = H’(metabarcoding sample).

The values of the Kullback-Leibler Divergences (i.e. the amount of information that is needed to transform the relative abundance distribution of species obtained with metabarcoding data into the original distribution of the mock communities) for the four treatments (two buffers, two incubation times) were quite similar regardless of whether community composition was based on biomass or specimen number (Fig. 6). For both 16S and COI, only Buffer had significant effects in the two cases (Suppl. material 7: Tables S11–S14). In the case of 16S, those replicates incubated in buffer B1 for LT1 had the lowest divergence values to the mock communities, although the interaction between Buffer and Lysis time was not significant, and a LT2 in B1 or incubation in buffer B2, increased the divergence. For COI, incubation time did not induce any difference in the value of the Kullback-Leibler Divergences, whereas the buffer in which the replicates were incubated did.

Figure 6.

Kullback-Leibler Divergences between the true community composition and the metabarcoding estimates of it. Community composition is measured in terms of biomass (left) or the number of specimens (right). Data are shown both for the 16S marker (blue) and the COI marker (orange). For the 16S marker, the divergence between the metabarcoding and the original sample increases with buffer aggressiveness and incubation time, while for the COI marker, the divergence is only affected by an increase in buffer aggressiveness.

As far as we could judge, all lysis protocols applied here essentially left insect morphology intact. Much of this beneficial outcome may be due to the limited lysis time used, since the digestion step of each protocol here examined was less than 5 hours. The incubation times are, thus, much shorter than those used in other protocols for terrestrial insect samples, which range from around 14 hours of lysis to up to 72 hours (Vesterinen et al. 2016; Nielsen et al. 2019; Ji et al. 2020). The lysis times used here are more similar to those previously applied to aquatic invertebrates (Carew et al. 2018). In previous pilot tests (unpublished), we applied much longer digestion times of 48 – 72 hours, but then found them to be highly damaging for the insects, to an extent where identification even to order proved difficult and, in some cases, impossible.

The high level of morphological preservation here achieved is hope-inspiring. From a taxonomist perspective, it allows the later description of new species from the material treated. As an exciting scenario, we may then apply bulk metabarcoding to generate taxonomic lists of contents for large sets of bulk samples. Such lists may then be offered to expert taxonomists, allowing them to direct their input to those samples offering the highest reward in terms of new and interesting species to examine. This is a quantum leap from the tedious manual sorting of mass samples, where the main effort typically goes into dealing with the most abundant and typically less interesting taxa. Such tasks represent the poorest possible use of skilled taxonomists, whose availability tends to be in short supply.

All methods of DNA extraction tested here yielded DNA of sufficient concentration and quality for successful PCR and sequencing. DNA concentration was consistently higher for the replicates of which lysates had been purified with the salt saturation method, compared to those in which the extraction was done using silica-coated magnetic beads in a robot. This is not unexpected, as the salt saturation method started from a lysate volume of 7.5 mL, whereas the robot-based method started from only 225 µL, with both methods ending at a final elution volume of 150 µL. Thus, it is likely that the difference in initial quantity of DNA is reflected in the final concentration. In addition, it is important to note that the amount of beads in the reaction was the same for all four treatments, which might explain why the longer lysis times did not increase the DNA concentration more. In any case, the results are intuitive: a more chemically aggressive buffer, a longer digestion time and a larger input volume will all increase the concentration of DNA in the extract.

In terms of sample purity, the overall values achieved were high. Compared to the value of the A260/280 ratio considered ideal (1.8–2.0), we found highly adequate readings (1.5–1.9). In three of the four lysis treatments, the replicates purified with the salt saturation method produced a higher A260/A280 ratio, but those corresponding to buffer B2 and short incubation time showed the opposite relation. This could possibly be explained if one assumes that the higher concentration of proteinase K and the presence of DTT and Ca²⁺ in buffer B2 released more proteins to the lysate than buffer B1, but that the short incubation time was not enough to hydrolyse these proteins completely. However, in general, a longer incubation time produced DNA extracts with higher purity, same as the manual purification with the salt saturation method. Although significant, these differences had only a small effect on the species recovery. This differs from a previous study, in which the salt saturation method was shown to provide metabarcoding data with higher species richness than those provided by commercial kits(Kaunisto et al. 2017). However, the latter study was based on faecal material, and the results must, thus, be taken with caution. Since the contents of faeces depart from that of bulk insect sample, several other factors may bias the output (see, for example, Nielsen et al. 2018). Our results are in line with a more recent study (Nielsen et al. 2019) that used mock insect community samples and detected no differences in species recovery between the methods. As long as the lysate is well mixed, there is no indication that diversity recovery is affected by volume subsampling.

In our study, all analyses were based on communities of known composition. In terms of species recovery, the COI marker was able to detect all 23 species we used in the communities, albeit not all species were detected in all the samples where they were present. For instance, small and delicate insects like Tuberculatus annulatus (Hemiptera) and Chrysoperla carnea larvae (Neuroptera) yielded low read abundances in most samples and appeared missing from many of those that were subjected to lysis with buffer B2. The 16S marker failed to detect the isopod Porcellionides pruinosus and the two species of Formica (Hymenoptera). The absence of P. pruinosus is not surprising, as the 16S primers used in this study had low degeneracy and were designed to target only insects (Marquina et al. 2019a). The case of Formica seems to be attributable to the resistance of the cuticle to mild lysis, as it is missing from the 16S dataset, but also very seldom recovered in the COI dataset, in combination with low binding affinity to both 16S and COI primers.

Importantly, the current study was explicitly aimed at evaluating the effect of the extraction protocol on the detectability of species against a community background of standardised complexity. Our communities were varied by excluding a single species amongst 23, whereas we did not vary the background complexity from highly species-poor to highly species-rich samples. As variation in the latter dimension provides an important aspect of natural communities, its effect should be the target of future studies. What we do see is that species detection rates, even against a standardised background, will never reach 100%. This pattern matches that reported by other studies. When using communities of known composition, both Carew et al. (2018) and Nielsen et al. (2019) reported significant variation in detection success, which they, too, attributed to a mix of effects in unknown proportions: species-specific characteristics, PCR biases and lack of homogeneity of the DNA extract. Thus, species-to-species variation in detectability remains an important challenge in the metabarcoding-based characterisation of insect communities, no matter which extraction method is used.

Accurate abundance estimation is currently one of the main research fronts in metabarcoding. Early studies suggested that metabarcoding was unsuitable for quantification, and that accurate abundance estimates might only be achieved through shotgun sequencing using mitochondrial metagenomics (e.g. Crampton-Platt et al. 2016; Bista et al. 2018; reviewed in Lamb et al. 2019). Yet, recent developments using correction factors or spike-ins of known concentration (Krehenwinkel et al. 2017; Ji et al. 2020; Ershova et al. 2021) have yielded promising results. A complicating factor that will undoubtedly introduce additional noise when using mitochondrial markers to estimate abundance is that the mitogenome copy number may vary considerably between specimens due to size differences and to differences in physiological state (ratio of mitochondrion-rich to mitochondrion-poor tissue). However, this would affect both metagenomics and metabarcoding.

Importantly, the current study focuses on the impact of a single step in sample processing: that of the lysis phase of DNA extraction. What we find is some factors that clearly contribute to a poor overall relationship between species abundance and read abundance For instance, the failure of the metabarcoding data to estimate the specimen counts of the different species appears to be due to a large extent to the over-representation, at least in some treatments, of DNA from large species that were represented by only one or a few specimens (see, for example, Calliphora, Acheta and Locusta in Suppl. material 7: Figs S3, S4). It thus seems possible that the biases are consistent and that, given enough training data, a machine-learning model could compensate for the variance in read abundance for such taxa. If this proves the case, then it will be possible to improve the estimation of specimen counts considerably for other taxa.

In terms of other future improvements, it is quite plausible that mild lysis protocols can be further optimised to more accurately represent the contents of the sample processed. In our experiment, a moderately chemically aggressive lysis buffer and a short incubation time tended to reduce the difference in estimates of alpha diversity compared to the actual communities significantly more than a more destructive buffer or longer incubation times, showing that even when no precise estimates about species abundance distributions can be obtained, still some ecological insight can be drawn using this method. This likely illustrates simple considerations based on the relation between body volume and surface. During the early part of the lysis, both large and small insects presumably release DNA from their tissues in contact with the buffer, at a rate proportional to the exposed surface (roughly equivalent to the square of the body size). As the incubation time increases, the digestion will continue towards the internal tissues and, thus, the released DNA will be proportional to the volume of the individual (roughly equivalent to the cube of the body size) (Nielsen et al. 2019). Thus, it seems likely that larger insects will contribute proportionally more to the DNA pool the longer the lysis period or the more invasive the digestion buffer. These predicted patterns fit the observed Kullback-Leibler Divergences well, in that those replicates incubated in the gentler buffer and for a shorter period of time were more similar to the mock samples than those incubated in a more aggressive buffer or for a longer time. In this sense, longer incubation times would approximate the proportions of DNA from small versus large individuals obtained in homogenisation protocols, in which larger individuals contribute proportionally more to the DNA pool than smaller ones (Elbrecht et al. 2017b).

As a final caveat, we would like to re-emphasise that our results are based on a series of mock communities, the complexity of which is drastically lower than many real samples from Malaise traps or other efficient insect traps. In the future, experiments similar to the ones here conducted should thus be aimed at varying other aspects of community context, including significantly more diverse samples (Creedy et al. 2019), to verify that our results are consistent and scalable. Of course, the costs for the lysis will be higher for real Malaise trap samples (especially those from habitats with specimen-abundant insect faunas), but still affordable. We estimate that the current cost of lysis and DNA purification would be ~ 4.5 €/sample for a lysis volume of 200 mL, and ~6 €/sample for a volume of 500 mL. The former volume will generally suffice for a typical Malaise trap sample, in our experience. For destructive sampling, extraction can be done for a smaller volume of biomass after homogenisation, thus lowering the reagent cost; on the other hand, the costs of equipment and tubes required for the homogenisation, or additional labour required, may not be trivial. However, we feel that any additional cost of mild lysis is more than outweighed by the benefits of being able to retain the specimens for later morphology-based work or single-specimen sequencing.