Pooling size sorted Malaise trap fractions to maximize taxon recovery with metabarcoding

Introduction

DNA metabarcoding is a useful approach to characterize arthropod communities. Instead of DNA barcoding individual specimens (Hebert et al., 2003), DNA is usually extracted from homogenized bulk samples (Zizka, Geiger & Leese, 2020; Mata et al., 2020; Turunen et al., 2021; Hardulak et al., 2020). A barcoding marker is amplified for the whole community and sequenced with high throughput sequencing. This allows for the identification of entire communities (Ji et al., 2013), often at species level resolution, depending on marker choice and reference database completeness (Leray et al., 2019; Weigand et al., 2019). However, DNA metabarcoding is affected by a large number of methodological biases preventing accurate specimen counts and often not all taxa in the sample are detected (Elbrecht et al., 2017; Piper et al., 2019). While there is a loose correlation between biomass and read abundance (Lamb et al., 2019), metabarcoding results are often skewed by primer bias (Elbrecht & Leese, 2015; Piñol et al., 2015; Krehenwinkel et al., 2017), mitochondrial copy number variation (Braukmann et al., 2019; Deagle et al., 2019), sequencing errors (Singer et al., 2019) and variation in specimen size, which can result in smaller taxa remaining undetected (Elbrecht, Peinert & Leese, 2017).

Arthropod bulk samples often show a wide range of specimen sizes (Aylagas et al., 2016). While imaginal stages of the same species have similar biomass, there can be substantial size variation between the life stages within species and of course between different species (Stork & Blackburn, 1993; Siemann, Tilman & Haarstad, 1996). Thus, when extracting DNA from bulk samples containing rare or small taxa with low biomass, those might be overshadowed by biomass rich taxa that contribute a large amount of DNA to the extraction. This natural variation can jeopardize the detection of small and rare specimens (Deagle et al., 2018), especially when taking PCR effects into account (Kelly, Shelton & Gallego, 2019). To maximize taxon detection in bulk samples showing substantial specimen biomass variation (Aylagas et al., 2016), size sorting might be required (Liu et al., 2020). Metabarcoding of marine arthropod bulk samples is commonly preceded by size fractioning of samples through sieving (Cowart et al., 2015; Leray & Knowlton, 2015; Wangensteen & Turon, 2017; Wangensteen et al., 2018). Size sorting has also been applied to freshwater macroinvertebrates (Elbrecht, Peinert & Leese, 2017; Carew, Coleman & Hoffmann, 2018; Zizka et al., 2019) and terrestrial arthropods (Ji et al., 2013; Creedy, Ng & Vogler, 2019; Hausmann et al., 2020). While the advantages of size sorting seem logical, the process can take a substantial amount of time if done manually (Elbrecht, Peinert & Leese, 2017) or even when using semi-automated methods like sieving. Additionally, while size sorting likely increases taxon recovery, it also generates sub samples that can either be pooled or sequenced individually. Regardless of the chosen size sorting and metabarcoding strategy, the amount of laboratory work and cost increases compared to the processing of a single metabarcoding bulk sample without size sorting. Thus, some authors suggest increasing sequencing depth instead of applying size sorting strategies, especially when size variation is limited (Creedy, Ng & Vogler, 2019; Elbrecht, Peinert & Leese, 2017). In fact, most non-marine studies do metabarcode complete bulk samples without prior size sorting (e.g. Hardulak et al., 2020; Hausmann et al., 2020; Steinke et al., 2020; Bush et al., 2020).

In this study, we explore the effects of size sorting Malaise trap samples through sieving into four size fractions to maximize taxon recovery. We thereby want to circumvent time consuming and therefore expensive manual sorting, which is not feasible for large-scale approaches. Furthermore, manual sorting and determination of samples is difficult to standardize and depends on the researcher performing the sorting. We therefore focus on reproducible sieving approaches including a standardized mesh size and a defined sieving time and amplitude. Prior to DNA extraction, we specifically test the effect of pooling lysates of each size fraction in different proportions. Each size fraction is also extracted and metabarcoded individually, giving us the opportunity to further validate different size sorting scenarios. The goal of this study is to find a standardized and scalable DNA metabarcoding strategy for Malaise trap samples, which also captures the diversity of small and rare specimens.

Materials and Methods

Tissue homogenization, lysate pooling and DNA extraction

To reduce electrostatic charging after dry sieving, sieved insects were moistened with 96% ethanol before transferring them into 30 ml grinding tubes (Nalgene; Thermo Fisher Scientific, Waltham, MA, USA). Specimens were dried again for up to 5 days at 50 °C in the oven (INCU Line ILS 6, VWR, Radnor, PA, USA). Homogenization was performed in grinding tubes with five mm diameter steel beads (Retsch, Haan, Germany), using a mixer mill (MM400; Retsch, Haan, Germany) for 2 min at 30 hz. Ground tissue from each size fraction was transferred to either 50 ml falcon tubes or 500 ml bottles and weighed. For tissue lysis, 10 ml ATL lysis buffer (Qiagen, Hilden, Germany) supplemented with 2% Proteinase K (20 mg/ml) was added per g of tissue. Ground tissue with lysis buffer was incubated overnight at 56 °C in a shaking incubator (INCU Line ILS 6; VWR, Radnor, PA, USA) and centrifuged for 5 min at 4,149 g (Mega Star 1.6; VWR, Radnor, PA, USA). Samples were centrifuged to pelletize solid tissue at the bottom of the tubes and to facilitate the pooling of lysates by pipetting the respective proportions of supernatant. Thereby, lysates of the four size fractions (S < 2 mm, M > 2 mm, L > 4 mm, XL > 8 mm) were pooled in volumes ranging from 0.18 to 200 µl of supernatant (see Fig. 1 & Table S1 for details of lysate volumes pooled for the different ratios). In total, 200 µl of pooled (or single fraction) lysate was used for each DNA extraction. Since the amount of extraction buffer added was proportional to the weight of tissue in the respective size fraction, we assumed the DNA concentration in each lysate to be the same. Thus, pooling lysate volumes of the four size fractions proportional to the tissue weight simulates an unsorted Malaise trap sample (L1_gA, L2_gA, L3_gA). The other schemes (shown in Fig. 1, Table S1) simulate different pooling strategies of size fractions to test the effects on taxon recovery. In addition, individual size fractions (S, M, L, XL) of all three samples were also processed and sequenced individually and used for later in silico analysis.

DNA was extracted from the lysate using a Qiagen DNeasy 96 well Blood & Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol with only one adjustment (2% proteinase K). DNA was eluted with 200 µl buffer AE (Qiagen, Hilden, Germany). The 96 well plate map with sample locations is available as supporting information (Fig. S2), with 27 wells being used for primer tests as part of a separate project. In nine negative controls only ATL lysis buffer (Qiagen, Hilden, Germany) and Proteinase K without tissue were used for extraction (DNeasy 96 well Blood & Tissue Kit; Qiagen, Hilden, Germany). The eluates obtained from these negative controls were also included in downstream amplicon and indexing PCR. DNA quality was examined via gel electrophoresis on a 0.5% agarose gel with 0.1% GelRed (10,000×; Biotium Inc., Fremont, CA, USA) and 1× TBE. One µl of GeneRuler 100 bp plus DNA ladder (Thermo Fisher Scientific, Waltham, MA, USA) and three µl of each DNA extract were loaded onto the gel.

In silico analysis

To test if the partitioning of a sample into two size fractions would be sufficient to achieve optimal taxon recovery, the individually sequenced size fractions S+M (S < 2 mm, M > 2 mm) and L+XL (L > 4 mm, XL > 8 mm) were combined and pooled in silico in different ratios (Fig. 2, Fig. S3). First, based on data obtained from the separately sequenced size fractions, sequences of class S and M were pooled in silico for each sample after quality filtering (sequences present in both replicates, representing at least 0.01% of total observed reads per sample). Pooling was conducted to simulate unsorted samples, according to the dry weight of each size fraction. An identical process was carried out for the sequences from size fractions L and XL, in which they were pooled to simulate unsorted samples (Fig. S3). In the next step, the two reconstructed fractions (S+M and L+XL) were pooled in silico in different proportions, based on a sequencing depth of 10,000 reads. Starting with 100% S+M and 0% L+XL, the ratio was altered in 5% steps to 0% S+M and 100% L+XL, taking into account read numbers (e.g. 10 × read abundance S+M, 0 × read abundance L+XL). Relative read numbers were calculated for merged samples and a filtering threshold of 0.01% was applied. Thereby OTUs with a read abundance <0.01% per sample were discarded. Detailed information about in silico pooling is available in Supplementary Material Script S2 (“silico_pool_v1.R”).

Rarefaction was applied to the unfiltered dataset (sequences present in both replicates, representing at least 0.01% of total observed reads per sample). The following examples detail a simulated sequencing depth of 500 reads (lowest sequencing depth) and 50,000 reads (highest sequencing depth (Fig. 3)). First, the sequencing depth was calculated using a filtering threshold defining the minimum abundance an OTU should have to be included in the analysis (for a sequencing depth of 500 reads: 100/500 = 0.2%; for a sequencing depth of 50,000 reads: 100/50,000 = 0.002%). In the next step, OTUs were excluded from the samples (g, equ, invg, fibo, log, 4×) if the minimum abundance of 0.2% (for 500 reads) or 0.002% (for 50,000 reads) was not reached respectively. Detailed information about in silico pooling is available in Supplementary Material Script S3 (“rarefaction.R”).

Results

Sorting dried specimens into four size fractions using the “sieving tower” (Fig. S1A) took 2–3 h per sample in total as statically charged specimens and broken off tissue pieces had to be collected from the sieve surface. The alternative “shaker” design was also affected by electrostatic charging. This problem could be overcome by moistening the dried samples before transferring them, which however again results in a higher total processing time. Wet sieving was not affected by electrostatic charging of specimens and enabled size sorting in approximately 15 min including sample transfer. Alternating between gentle shaking and letting small specimens sink through the four mm mesh resulted in clear size fractionation with one order of magnitude difference in average specimen weight between the small and large size fractions (Fig. S1B).

The MiSeq run generated 21,911,822 PE reads which are available as raw data as well as demultiplexed samples under the NCBI SRA accession PRJNA625408. No-call events were observed on several tiles in read 2, leading to N being inserted at positions 16, 17 and 25. A total of 20 out of 38 tiles were affected by one or more No-call events, but as respective areas were mostly removed by primer trimming, analysis was not impaired. Average sequencing depth per sample was 249,423 reads (SD = 92.075) with the lowest sample containing 146,344 reads. In bioinformatic processing on average 23.40% (SD = 4.16%) of reads were discarded. After complete quality filtering, 99% of the 2,197 OTUs were assigned to the phylum Arthropoda (Table S2).

To assign a specimen size to each OTU, relative read abundance of all three samples was combined for each size class. Then, the size class with the highest abundance was selected as a specimen size estimate. Three quarters of OTUs were present in only one size fraction (Figs. S3, S4, S5). Thus, most OTUs could be clearly assigned to a size class using this strategy, with a few exceptions in the large and medium size category (e.g. OTU1 Myrmica ruginodis and OTU11 Sphaerophoria). Here, the size category with the highest read number was assigned, even when an OTU was found in high abundance in two fractions. For some low-abundance OTUs, no reads were assigned in any of the separately sequenced size classes, but only in the pooled samples (this was the case for 9.38% of OTUs). Thus, the respective OTUs had no size class assigned. Most OTUs were assigned to smaller size classes, with 41.51% in S, 24.53% in M, 19.21% in L and 5.37% in XL (Table S1). The individual size fractions (S, M, L and XL) recovered mostly reads from their respective size class (Fig. 2), with only a few exceptions (e.g. OUT1, Myrmica ruginodis). Here a high number of reads was assigned to size fraction L and M (43.2% and 55.7%). Thus, for OTU1 the size class M was assigned, despite the large size fraction also containing a substantial number of reads. The total number of taxa recovered was the highest in the smaller size fractions, and no more than 100 OTUs were recovered in the extra-large size fraction (Fig. 2).

When lysates were pooled based on dry size fraction weight (g) (simulating an unsorted Malaise trap sample), results were dominated by large specimens. Compared to the other pooling strategies, the simulated unsorted samples showed the weakest OTU recovery, followed by the pooling of all four lysates in equal proportions (equ). The other four pooling methods (invg, fibo, 4× and log) recovered more OTUs, with equally high numbers detected across the three Malaise trap samples. The logarithmic pooling (log) recovered the highest proportion of small specimens, while showing reduced recovery of large ones compared to the other methods. This was also supported by Jaccard and Bray–Curtis dissimilarity analysis (Fig. S6). Here, especially samples pooled using the log and 4× strategies cluster closely with samples where size class S had been sequenced individually. Samples pooled using strategies equ or g showed high similarity to size fraction L or even XL when read numbers were taken into account.

In silico pooling of the four individually sequenced size fractions at 1/4 sequencing depth of the lysates pooled in equal proportions (equ), showed that the individually sequenced size fractions recovered an average of 26.56% more OTUs than the equally pooled lysates (Fig. 2). However, the numbers of OTUs recovered from the in silico pool of individually sequenced size fractions were comparable to OTU counts recovered by lysates pooled in favour of small size fractions.

Rarefaction analysis showed that samples pooled in favour of small size fractions (invg, fibo, 4×, log) had a similar number of OTUs at around 15,000 reads sequencing depth compared to what was recovered from the unsorted sample (g) at 50,000 reads sequencing depth (Fig. 3). This represents a factor of approximately 3–4× reduction of sequencing depth compared to unsorted samples. In this comparison, unsorted samples (g) were more biased towards large specimens, while the samples pooled in favour of small specimens showed a higher detection rate of taxa assigned to the small size class (Fig. 2). Thus, while the number of OTUs recovered from the unsorted samples at 3–4× sequencing depth might be similar to that of size sorted samples, there are still differences in OTU composition between different pooling strategies (Fig. 2, Fig. S6).

To test if separation of a sample in only two size fractions would be sufficient to obtain optimal taxon recovery, individual size fractions S+M and L+XL were combined and pooled in silico in different ratios (Fig. S3). The recovery of OTUs was highest when adding 10–20% of the large fraction to the small fraction. The overall OTU recovery in the in silico pooled samples did in fact exceed the number of OTUs recovered by the processed lysate pools (invg, fibo, 4×, log) by up to 10–20%.

Discussion

Previous studies have demonstrated that size sorting before metabarcoding can improve taxon recovery (Elbrecht, Peinert & Leese, 2017; Creedy, Ng & Vogler, 2019). This is particularly relevant when dealing with bulk samples containing specimens with a wide variation in biomass. We were able to confirm these trends for Malaise trap samples. Experiments were done on three Malaise samples only, which represents low replication effort and thus does not allow for statistical testing of results. However, the results presented here clearly support expectations and findings of previous studies (Elbrecht, Peinert & Leese, 2017; Creedy, Ng & Vogler, 2019). As illustrated in Figs. 2 and 3, our results present distinct patterns clearly showing differences between the applied strategies and therefore justifying the limited number of replicates and the absence of statistical testing.

Small and medium sized fractions contained the majority of OTUs in each of the three samples, consistent with studies from marine benthos (Wangensteen et al., 2018). When the four size fractions were pooled in different proportions, the fractions simulating a completely unsorted Malaise trap were dominated by large specimens, showing consistently the lowest recovery of OTUs. Interestingly, all four ways of pooling in favour of the smaller fractions (inverse weight, logarithmically, fibonacci based and 4-fold increase) recovered a similar number of OTUs. This result was unexpected, as for example logarithmic pooling strongly favours small fractions while severely underrepresenting larger ones. A reason for this effect could be OTU sharing across the individually sequenced size classes due to small specimens attached to larger ones, consequently also detected in larger fractions. This would also indicate that sieving does not perfectly separate individuals of different size classes. However, for the dry sieving approach, the aggregation of specimens of different size classes was rarely observed and counteracted through the shaking amplitude. The phenomenon occurred more often with wet sieving but was counteracted with an increased mesh surface, the use of higher ethanol volumes during sieving and a higher sieving amplitude. Furthermore, smaller specimens may be detected as secondary prey items in the guts of larger specimens and would therefore also be detected in larger size fractions. Additionally, legs and other body parts can easily break off the brittle specimens during dry sieving, and small size fractions therefore contain tissue from large specimens. Because small pieces of tissue will not break off from all specimens equally and large taxa found in the small fractions might not be representative of the real diversity in the samples, specimen damage during sieving and the subsequent detection of larger specimens in smaller size classes should be avoided. As shown here, wet sieving is a promising approach which is more gentle than dry sieving. In addition, different size fractions from different samples should be pooled in the same proportions for each sample to ensure comparability. Choice of pooling strategy might depend on study goals, and on size class and taxonomic composition of the samples, which can vary greatly depending on location. If taxa of interest fall within a particular size range, pooling should be adjusted accordingly, given that the number of reads and OTUs of differently sized specimens depends on the pooling process.

While we have shown that size sorting can increase taxon recovery in Malaise trap samples, it remains to be discussed how this can be done most efficiently. This is of particular importance for upscaled metabarcoding as required for biodiversity monitoring. Size sorting can only be used in larger scale metabarcoding projects if the process is time and cost efficient and does not substantially increase laboratory work. When the number of samples is limited and specimen biomass does not vary substantially, biomass biases can be mitigated (Creedy, Ng & Vogler, 2019) or minimized by only including a small part of the tissue (e.g. one leg) from large specimens (Ji et al., 2013; Carew, Coleman & Hoffmann, 2018). While manual sorting and processing of specimens can be a good strategy to reduce biomass bias, it is rather work intensive and difficult to scale or standardize. An automated approach would be desirable, but automated sieving techniques come with their own challenges. With the dry sieving approach used in this study, major methodological limitations arose, as it took several hours per sample to pick electrostatically charged individual specimens from the equipment. Consequently, we recommend using the alternative approach of wet sieving, completely immersing the specimens and sieve in ethanol. This approach can be automated more easily by stacking the sieve containers on shakers and thereby reducing hands-on time and increasing scalability. However, it needs to be considered that intensive cleaning of the equipment is required between sample processing steps to prevent cross contamination. Here, at least 15 min per sample should be calculated. Furthermore, wet sieving is conducted in pure 96% ethanol and the costs for additional amounts of ethanol need to be considered. Wet sieved specimens are less prone to damage and might be usable for morphological identification if non-destructive metabarcoding methods are applied, using for example preservative ethanol (Hajibabaei et al., 2012) or lysis buffer protocols (Nielsen et al., 2019).

Our in silico evaluations suggest that size sorting with a single sieve of four mm mesh size represents a good trade-off between taxon recovery and laboratory workload. Pooling the large and small size fractions in a proportion of 1:10 to 2:10 showed the highest taxon recovery in our in silico test, which however should be confirmed experimentally, since in silico comparisons can have limited reliability and due to the low replicate number processed in this study. For example, when in silico pooling 1/4 of each of the individually sequenced size fractions, more taxa were recovered than sequencing equally pooled lysate proportions. Thus, taxon recovery for the in silico estimates is slightly overrepresented, probably due to primer bias effects (Elbrecht & Leese, 2015; Piñol et al., 2015). For example, if a specific taxon in the large fraction is very well amplified, it will make up a large proportion of the reads in the pooled sample, decreasing the sequencing depth of all other specimens. When individual size fractions are sequenced, primer biases are limited to each individual fraction, increasing overall taxon recovery which is why some authors recommend the taxonomic sorting of bulk samples and sequencing each taxonomic group individually (Morinière et al., 2016; Beentjes et al., 2019). A taxonomic sorting step is however impractical for metabarcoding, as it substantially prolongs processing time per sample and does not necessarily increase taxon recovery as most arthropod primer binding sites show the same variability pattern (Elbrecht et al., 2019).

Other studies implement individual sequencing of each size fraction (Cowart et al., 2015; Wangensteen et al., 2018) and provided information about community composition for each size class to increase the detection of small and rare taxa. Sequencing of separate size classes doubles the costs for sample preparation and sequencing per sample. As an alternative, lysates or extracted DNA can be pooled proportionally to enrich the small size fractions sufficiently. Thereby, size sorting and subsequent pooling of lysates of each fraction before the DNA extraction is particularly advantageous as it reduces costs and time required per sample compared to pooling strategies after DNA extraction or PCR.

In our study, we confirm that size sorting can be applied to increase taxon recovery in metabarcoding of Malaise traps samples. However, the implementation of techniques, e.g. number of size fractions or pooling proportions, needs to be adjusted to the aim of the study. In some cases, size sorting may be skipped, as a sufficient signal can be achieved with a high sequencing depth for samples with limited specimen size variability (Creedy, Ng & Vogler, 2019). While Creedy, Ng & Vogler (2019) do not give exact numbers on sequencing depth effects, we simulate that a 3- to 4-fold increase in sequencing depth recovers similar numbers of OTUs from an unsorted sample compared to size sorted and proportionally pooled samples. Calculations are based on rarefaction analysis, which is explained in the materials and methods section (“In silico analysis”) and illustrated in Fig. 3. However, by increasing the sequencing depth, metabarcoding results will still be dominated by biomass rich taxa and it should be further experimentally tested how a higher sequencing depth can adequately replace size sorting.

Given that size sorting can be done in around 15 min using the wet sieving method demonstrated here, it seems like a worthwhile step to prepare Malaise trap samples for metabarcoding. Further Malaise trap sampling and additional size sorting experiments could be useful to quantify effects of large and abundant species on the detection of small or rare ones. Examination of Malaise trap samples or mock communities with known taxonomic composition and assigned size classes would be ideal to evaluate the results presented here. In addition to in silico assays, experimental approaches testing how a higher sequencing depth can compensate for size sorting and how different pooling proportions affect detected diversity patterns would be helpful to develop standardized monitoring protocols for flying arthropods.

Conclusions

We demonstrate that size sorting of Malaise trap samples can increase recovery of small specimens. While our analyses show distinct patterns clearly supporting previous studies, the limited number of samples processed and the resulting lack of statistical support needs to be carefully considered. Pooling lysates in the extraction process can be an efficient way of improving the representation of small taxa. The methods presented here reduce laboratory work and costs, in comparison to sequencing each size fraction individually. In silico tests suggest that sieving into two size fractions (four mm mesh size) and subsequent pooling in a proportion of 1:10 (one part size class L and nine parts size class S) or 2:10 can substantially increase taxon recovery. Increasing sequencing depth by 3 to 4-fold can also be an alternative to size sorting, but might still underrepresent small and rare taxa due to biomass biases. Wet sieving of specimens in ethanol may be the best compromise for optimizing taxon recovery with low effort, compared to other approaches. However, this method remains to be thoroughly tested. Ultimately, the decision of whether or not samples should be size sorted depends on the research goals, budget and timeline as well as on the sample composition and size heterogeneity of specimens in the samples.

Supplemental Information