Received: 23 September 2020 | Revised: 21 June 2021 | Accepted: 7 July 2021 DOI: 10.1111/1755-0998.13464 RESOURCE ARTICLE Relative species abundance estimation in artificial mixtures of insects using mito-­metagenomics and a correction factor for the mitochondrial DNA copy number Lidia Garrido-­Sanz1 | Miquel Àngel Senar1 | Josep Piñol1,2 1 Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Spain Abstract 2 CREAF, Cerdanyola del Vallès, Spain Mito-­metagenomics (MMG) is becoming an alternative to amplicon metabarcoding for Correspondence the assessment of biodiversity in complex biological samples using high-­throughput Lidia Garrido-­Sanz, Universitat Autònoma sequencing. Whereas MMG overcomes the biases introduced by the PCR step in the de Barcelona, Cerdanyola del Vallès 08193, Spain. generation of amplicons, it is not yet a technique free of shortcomings. First, as the Email: Lidia.Garrido@uab.cat reads are obtained from shotgun sequencing, a very low proportion of reads map into Funding information the mitogenomes, so a high sequencing effort is needed. Second, as the number of Spanish Government grant, Grant/Award mitogenomes per cell can vary among species, the relative species abundance (RSA) Number: TIN2017-­8 4553-­C2-­1-­R; Agència de Gestió d'Ajuts Universitaris i de in a mixture could be wrongly estimated. Here, we challenge the MMG method to es- Recerca (AGAUR), Grant/Award Number: timate the RSA using artificial libraries of 17 insect species whose complete genomes 2017-­SGR-­1001 are available on public repositories. With fresh specimens of these species, we cre- ated single-­species libraries to calibrate the bioinformatic pipeline and mixed-­species libraries to estimate the RSA. Our results showed that the MMG approach confidently recovers the species list of the mixtures, even when they contain congeneric species. The method was also able to estimate the abundance of a species across different samples (within-­species estimation) but failed to estimate the RSA within a single sam- ple (across-­species estimation) unless a correction factor accounting for the variable number of mitogenomes per cell was used. To estimate this correction factor, we used the proportion of reads mapping into mitogenomes in the single-­species libraries and the lengths of the whole genomes and mitogenomes. KEYWORDS Metazoa, mitochondrial genomes, mitogenome skimming, mock sample, next-­generation sequencing, PCR-­free 1 | I NTRO D U C TI O N MMG consists in the shotgun sequencing of a DNA sample followed by the mapping of the reads to mitochondrial genomes (hereafter, Mitochondrial metagenomics or mito-­ metagenomics (hereafter, mitogenomes) obtained from online repositories or ad hoc assem- MMG) is becoming an alternative to the classical amplicon metabar- blages (Crampton-­Platt et al., 2016). On the positive side, MMG coding for the large-­scale assessment of biodiversity of Metazoa avoids the amplification biases caused by the PCR step (Elbrecht & (Crampton-­Platt et al., 2015; Tang et al., 2014; Zhou et al., 2013). Leese, 2015; Piñol et al., 2015; Taberlet et al., 2012). On the negative This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Molecular Ecology Resources published by John Wiley & Sons Ltd. Mol Ecol Resour. 2021;00:1–15.  wileyonlinelibrary.com/journal/men | 1 2 | GARRIDO-­SANZ et al. side, the tiny size of the mitogenome compared to the nuclear ge- (there is a database of genomes sizes with 1344 insect species on nome produces a high number of non-­informative reads (Tang et al., it; Gregory (2020), accessed on 25 March 2020). Both the variation 2014), so a great sequencing depth is needed. across-­species of mitochondrial copy number and genome size affect It is generally assumed that MMG quantifies satisfactorily the MMG, but also any amplicon metabarcoding method that targets ge- relative species abundance (hereafter, RSA) of complex mixtures nomic regions with a variable copy number, such as COI in animals (Gómez-­Rodríguez et al., 2015; Zhou et al., 2013). However, as far as (Hebert et al., 2003), ITS in fungi (Schoch et al., 2012), or rbcL + we know, there are only five studies that tested the MMG method matK in plants (CBoL Plant Working Group, 2009). (plus one using chloroplast metagenomics in plants) using shotgun Here we explore the quantitative capabilities of MMG for the samples of known composition (Bista et al., 2018; Gómez-­Rodríguez estimation of RSA of heterogeneous mixtures of insects. For this et al., 2015; Gueuning et al., 2019; Ji et al., 2020; Lang et al., 2019; purpose, we prepared single-­species and artificial mixed-­species Tang et al., 2015). In general, the relationship between the expected libraries with several species of insects whose entire genome has and estimated RSA is statistically significant, but with high variability already been sequenced. We included four species of Drosophila in the goodness of fit. to assess the ability of the method to set apart closely related How the RSA of complex mixtures is presented in the literature species. The single-­species libraries allowed the calculation of a requires some clarification. First, the RSA can be expressed as a pro- reliable mitochondrial DNA copy number (NM) for each species portion of the species biomass (e.g., Gueuning et al., 2019) or indi- that was further used as a correction factor for the across-­species vidual counts (e.g., Lang et al., 2019); alternatively, the RSA can refer estimation of RSA of the mixed-­species libraries. In particular, we to the proportion of the DNA amount of each species in the mixture. addressed the following questions: (1) Is the MMG method able to Whilst the former approach is more meaningful for most ecological identify species in complex mixtures, even when they are of the studies, we adopt the latter approach here because it allows the in- same genus? As an approach to real samples, we investigated the dependent evaluation of different sources of bias on the RSA esti- robustness of the method in the absence of the mitogenome of mation. Second, some studies provide the relative abundance of one the focal species. (2) Can MMG estimate the RSA of complex sam- species in different samples (e.g., Bista et al., 2018), whereas oth- ples? Is it necessary the use of the N M correction factor for the ers report the abundance of several species in a single sample (e.g., across-­species estimation of RSA? (3) Finally, can the number of Saitoh et al., 2016). Ji et al. (2020) named within-­species estimation sequenced reads be reduced and still recover all species in a com- the former (is species i more abundant in sample s than in sample r?) plex sample of insects? and across-­species the latter (is species i more abundant than spe- cies j in sample s?). This distinction is important because there are species-­specific characteristics that influence the across-­species es- 2 | M ATE R I A L S A N D M E TH O DS timation but not the within-­species estimation. The most important of these characteristics is the variable number of mitogenomes per 2.1 | Selection of species, preparation of the DNA nuclear genome (mitochondrial DNA copy number). Thus, a species i libraries, sequencing, and quality control with twice the number of mitogenomes per nuclear genome than an- other species j will produce twice many mitochondrial reads as well; We selected 17 species of insects whose complete genome is al- without a proper correcting factor, species i would, apparently, be ready sequenced and available on the RefSeq repository (Table 1). twice more abundant in the mixture than species j. This fact is known Individuals of these species were captured alive or in fly traps from (Bista et al., 2018; Piñol et al., 2015; Tang et al., 2014), but there are various locations and DNA was extracted using DNeasy Blood & not reliable solutions to the problem because little is known about Tissue Kit (Qiagen) from c. 20 mg of fresh material. With the DNA the causes of the variation of the mitochondrial copy number across-­ extracts, we prepared two kinds of libraries: 21 single-­species librar- species (but see Liu et al., 2018, that reported a higher mitochondrial ies and six mixed-­species libraries. Four species were sequenced copy number in organs with a high metabolic rate and in species liv- twice in different runs using different extractions to test the repeat- ing at low altitude than in their counterparts at high altitude in the ability of the method (Table 1). The same extracts used for the first Tibetan Plateau). run of single-­species libraries were also used to create six artificial The size of the nuclear genome of the species also affects the mixed-­species libraries of 7–­8 species at known relative DNA con- across-­species RSA estimation. Being all other things equal, a spe- centrations to test the ability of the method to estimate the RSA cies r with a nuclear genome half as big as that of another species s (Table S1). In libraries no. 1 and no. 2, the RSA was highly variable will produce twice many mitochondrial reads because the mitochon- (from ~50%, the most abundant, to ~0.4%, the least abundant); in drial DNA is diluted in a smaller amount of nuclear DNA. Therefore, libraries no. 3 and no. 4, it was intermediate (from ~35% to ~3%); and without a proper correcting factor, species r would, apparently, be in libraries no. 5 and no. 6, the variability among species was much twice more abundant in the mixture than species s. The effect of lower (from ~24% to ~8.5%). As stated in the Introduction, here the genome size on RSA estimation is also known (Crampton-­Platt et al., RSA is the relative DNA concentration of the species in the mix- 2016; Krehenwinkel et al., 2017; Tang et al., 2014), but it is difficult ture, not their relative biomass. Consequently, all sources of varia- to consider it because measuring the genome size is not an easy task tion between the fresh biological material and the extracted DNA TA B L E 1 Summary information of single-­species libraries GARRIDO-­SANZ et al. Number of reads after Origin Number of raw quality control step Number of candidate Run Library Species Order Family Cultured/Wild (country) reads (paired-­end) (paired-­end) mito-­reads (paired-­end) 1 1 Papilio machaon Lepidoptera Papilionidae Wild Spain 434,520 432,712 14,774 1 2 Drosophila virilis Diptera Drosophilidae Cultured Spain 4,714,902 4,652,442 152,568 1 3 Drosophila Diptera Drosophilidae Cultured Spain 2,283,768 2,275,386 50,944 melanogaster 1 4 Drosophila Diptera Drosophilidae Cultured Spain 1,668,424 1,646,524 109,458 mojavensis 1 5 Bactrocera oleae Diptera Tephritidae Wild Spain 580,996 577,720 23,204 1 6 Linepithema humile Hymenoptera Formicidae Wild Spain 1,422,342 1,402,800 73,580 1 7 Bombus terrestris Hymenoptera Apidae Wild Spain 1,994,938 1,987,920 43,872 1 8 Apis mellifera Hymenoptera Apidae Wild Spain 1,262,388 1,236,630 223,202 1 9 Acyrthosiphon pisum Hemiptera Aphididae Cultured USA 684,688 596,972 104,344 2 1 Atta colombica Hymenoptera Formicidae Cultured Denmark 3,272,710 3,267,688 333,064 2 2 Bemisia tabaci Hemiptera Aleyrodidae Cultured Spain 2,513,212 2,501,586 23,096 2 3 Cimex lectularius Hemiptera Cimicidae Wild Spain 3,506,722 3,474,010 122,826 2 4 Drosophila Diptera Drosophilidae Cultured Spain 2,909,608 2,903,546 74,984 melanogaster 2 5 Drosophila Diptera Drosophilidae Cultured Spain 1,797,470 1,786,926 82,794 mojavensis 2 6 Drosophila virilis Diptera Drosophilidae Cultured Spain 1,336,884 1,332,888 49,146 2 7 Drosophila suzukii Diptera Drosophilidae Cultured Spain 2,510,384 2,501,862 79,498 2 8 Linepithema humile Hymenoptera Formicidae Wild Spain 2,164,404 2,157,154 95,936 2 9 Plutella xylostella Lepidoptera Plutellidae Wild Spain 4,250,124 4,244,328 85,882 2 10 Solenopsis invicta Hymenoptera Formicidae Wild Argentina 3,661,374 3,648,438 146,812 2 11 Vollenhovia emeryi Hymenoptera Formicidae Wild Japan 3,487,834 3,480,802 66,664 2 12 Wasmannia Hymenoptera Formicidae Wild Spain 3,335,212 3,293,182 179,000 auropunctata | 3 4 | GARRIDO-­SANZ et al. (e.g., DNA-­to-­biomass ratio) are ignored (Matesanz et al., 2019; Tang regions from the Mito1794 reference (hereafter, FilteredMito1794). et al., 2015). Low-­complexity regions were identified using dustmasker (-­level 45) All libraries were prepared with the TruSeq DNA PCR-­Free LT Kit (Morgulis et al., 2006) and replaced with Ns using an in-­house py- of Illumina following the manufacturer's instructions (Ref. 15037063) thon script. and sequenced using Illumina MiSeq with the 2 × 150 bp chemistry To avoid confusion, we recapitulate below the name and mean- in three different runs. Both the libraries and the sequencing runs ing of the three different databases of mitochondrial genomes that were also used in a previous study (Garrido-­Sanz et al., 2020). In the we used for the mapping of reads: parent study, the obtained reads were mapped against a reference database of whole-­genomes instead of a reference database of mi- • Mito1794: the original mitogenomes of 1794 species. togenomes as we do here. • FilteredMito1794: as Mito1794, but with the low complexity re- We pre-­processed the raw reads through a quality control step gions removed. using fastqc v0.11.7 (Andrews, 2015) and trimmomatic v0.36 (Bolger • Mito331: a subset of Mito1794 with only one mitogenome per et al., 2014) to trim the reads to a 150 bp length and to remove those family; this reference was only used to obtain the candidate mito-­ shorter than 140 bp. Only paired reads were kept. reads from the total of reads of each sample. As we did not know to which extent the filtering of raw reads 2.2 | Reference genomes and mitogenomes was useful, we conducted four different kinds of mapping of reads to reference mitogenomes in the single-­species We downloaded all mitogenomes of insects available at the NCBI libraries: RefSeq database on 1 August 2019, plus the complete genomes of the 17 species selected for the study from the same database. • Raw reads against Mito1794 reference database Ten species had the complete genomes but not the mitogenomes • Raw reads against FilteredMito1794 reference database on RefSeq, and we downloaded their mitogenomes from GenBank • Candidate mito-­reads against Mito1794 reference database (accessed on 2 August 2019). We used high-­quality mitogenomes • Candidate mito-­reads against FilteredMito1794 reference of as many species as possible to test the ability of the method to database find the selected species among the many others in the reference database. Species with several mitogenomes were deduplicated, and In all cases, the mapping was conducted using BWA with mem we obtained the mitogenomes of 1794 species of insects, compris- algorithm and default options. Because the mapping of a sample was ing 1174 genera, 331 families, and 27 orders (hereafter, Mito1794; done against each mitogenome individually, we obtained the same Table S2). number of SAM files as references used. We subsequently used SAMtools to remove reads that did not map to any reference (view -­F2308). 2.3 | Mapping of reads into references In MMG studies, a small proportion of shotgun reads map into the 2.4 | Assignment of mapped reads to species mitogenome (e.g., Tang et al., 2014), hence most reads are not useful and slow down the mapping process. Thus, it seems reasonable to In general, a read mapped to several reference mitogenomes (e.g., eliminate the reads that are not mitochondrial before the mapping homologous sequences in several species), so an algorithm was step (Crampton-­Platt et al., 2015, 2016; Zhou et al., 2013). For this needed to assign reads to species. For this purpose, we used the γ–­δ purpose, we created a reference database with one mitogenome per algorithm described in Garrido-­Sanz et al. (2020). Briefly, what the family (hereafter, Mito331) where the representative species per γ–­δ algorithm does is to quantify the similarity between the query family was chosen randomly. We then mapped the raw reads against read and the reference (i.e., the mapping ratio) and then decide the Mito331 reference data set using a permissive criterion and kept whether a read is informative or non-­informative. It is informative the putative mitochondrial reads (hereafter, candidate mito-­reads). when it is very similar (mapping ratio above γ) to species i and not The mapping was done using BWA 0.7.15-­r1140 (Li, 2013) with mem very similar (mapping ratio below δ) to the rest of species; in this algorithm and an alignment score of zero (-­T0). samtools 1.10 (Li et al., case, the read is assigned to species i. It can be non-­informative for 2009) was subsequently used to filter the paired-­end reads with no two reasons, either because the read is not similar enough to any mapping reads (view -­F2316 -­b) and recovered the mito-­reads in species (mapping ratio below γ for all species) or because it is too fastq format (bam2fq). similar (mapping ratio above δ) to two or more species; in this case, As low-­complexity regions are prone to misclassify the reads the read is discarded. In all cases γ > δ. (Lu & Salzberg, 2018; Pearman et al., 2019), we prepared a new set The γ–­δ algorithm has never been applied before to MMG of filtered mitogenome references by removing low-­ complexity data, hence the appropriate values of γ and δ are unknown, so the GARRIDO-­SANZ et al. | 5 single-­species libraries were used to find the best combination of the 2. The re-­ arrangement of Equation (1) allows the estimation of parameters γ and δ. The tested values were all the combinations of NMi. γ = {0.99, 0.98, 0.97} and δ = {0.98, 0.97, 0.96} provided that γ > δ. To find the best values of γ and δ we relied on the criterion that xi ⋅ Gi the number of recovered species had to be one in the single-­species NMi = (2) Mi ⋅ (1 − xi ) libraries. The reads in single-­species libraries were divided into a train- ing set with 75% of the reads randomly selected, and a test set 3. Gi and Mi are known for species with sequenced genomes, with the remaining 25% of the reads. The training set was used for but xi is not. In our experimental setting, xi can be estimated the calibration of the procedure; the test set was used to assess in the single-­species libraries as the ratio between the number the goodness of fit of the model and to calculate the summary of reads that map into the mitogenome (RMi) divided by the statistics. total number of reads of species i (RGi). A situation that can arise in real samples, as opposed to the ar- tificial samples used here, is that the mitogenomes of some species RMi in the sample are not in the reference database. We explored this xi = (3) RGi situation by running again the complete pipeline with all the single-­ species libraries using the best set of input data and parameters. Ideally, no read should be assigned to any species because the mi- We obtained RMi by mapping the reads of species i to its mitoge- togenome of the only species in the library is not in the database. nome when this mitogenome was the only one used as the reference However, the reads might eventually be wrongly assigned to other in the mapping. Regarding RGi, we assumed that all reads of the single-­ species in the database and, thus, generate false positives. The out- species library of species i belong to species i. come of this experiment should reveal the robustness of the γ–­δ al- In the comparison of actual versus estimated RSA using the gorithm in the assignment of reads to species. mixed-­species libraries, we multiplied the actual relative abundance of species i (Table S1) by NMi, and then renormalized the values to sum 1. 2.5 | Quantification of the RSA in mixed-­ Finally, we estimated the importance of knowing or ignoring the species libraries and the need for a species-­specific individual values of Gi and Mi of each species in the mixture in the es- correction factor timated RSA by comparing the results obtained with the correction factor of Equation (2) (NMi ) with another factor that uses the mean In the literature, the relative abundance of one species is sometimes value of G (G = 338 Mbp) and M (M = 16.3 kbp) for all the species in compared among different samples and on other occasions, the rela- the mixture (NMi): tive abundance of several species is compared within one sample (within-­species and across-­species RSA, respectively, following Ji xi ⋅ G NMi = (4) et al., 2020). Here, we present the comparison of actual versus esti- M ⋅ (1 − xi ) mated RSA in the mixed-­species libraries using both approaches. As we observed in the Introduction, from a conceptual point of view, 2.6 | Rarefaction of the input samples the quantitative estimation of within-­species RSA in MMG is easier than the across-­species RSA, because in the latter the mitochondrial We only multiplexed six mixed-­species libraries in a single Illumina DNA copy number can vary widely between species. MiSeq run (Table S1), with the consequence of a high economic With the single-­species libraries we estimated the mitochondrial cost per library. However, from a practical point of view, it would DNA copy number (NM) of each species in the following way: be interesting to use fewer reads per library and still have a good quantitative estimation of RSA. To test this possibility, we randomly 1. Let xi be the ratio of the genomic mitochondrial information rarefacted the mixed-­species samples at various proportions of the divided by the total (haploid) genomic information for species original number of reads {0.5, 0.1, 0.05, 0.01} and run the new data i. The mitochondrial information is the mitogenome length (Mi) sets through the entire pipeline. We repeated each simulation 100 times NMi; the total genomic information is the sum of the times using different subsets and from every simulation we recorded nuclear genome length (Gi) and the mitochondrial the number of recovered species. information. 2.7 | Statistical analyses and hardware Mi ⋅ NMi All statistical analyses were performed with r 3.4.2 (R Core xi = (1) Gi + (Mi ⋅ NMi ) Team, 2016) in rstudio 1.0.143 (RStudio Team, 2015). Plots were 6 | GARRIDO-­SANZ et al. created using the r packages ggplot2 (Wickham, 2016) and ggpubr focal species in all the libraries, whereas with the Mito1794 database (Kassambara, 2018). there appeared several false positives in many libraries (2.5 ± 1.2 We run the complete pipeline on a server with two Intel Xeon species per library using all raw reads or 1.7 ± 0.9 species using only E5-­2620 v3 processors with six cores each and hyperthreading tech- candidate mito-­reads; Table 2a). The masked regions mostly be- nology, so a maximum of 24 threads were available. longed to non-­coding regions of the mitogenome, including the con- trol region (Table S8). The use of only candidate mito-­reads instead of all reads produced a loss of c. 6% of informative reads (Table 2b) 3 | R E S U LT S but reduced ~18 times the execution time (Table 2c). In summary, the elimination of the repetitive regions from the genomes removed all 3.1 | Species identification the false positives and the mining of candidate mito-­reads reduced 18-­fold the execution time of the pipeline with a moderate loss of The 21 single-­ species libraries (Table 1) generated informative reads. Therefore, in the subsequent steps, we used both 2,371,091 ± 1,210,091 (mean ± SD) paired-­end reads. A propor- the FilteredMito1794 database and only the candidate mito-­reads tion of 0.012 ± 0.027 reads were eliminated in the trimming step, (Figure 1). remaining a proportion of 0.987 ± 0.027 reads available for further We evaluated the goodness of fit of the model with the test analysis. set (i.e., the remaining 25% of reads not used in the previous cal- The results that follow correspond to the application of the γ–­δ ibration) using the best set of input data and parameters (i.e., the algorithm for the assignation of reads to species on the training set FilteredMito1794 database, the candidate mito-­reads and the pa- (i.e., 75% of the sequenced data). We also eliminated from the fol- rameters γ = 0.99 and δ = 0.96). The number of identified species per lowing results the reads assigned to species that could legitimately library was 1 in all cases (Table S9) and the proportion of informative be attributed to physical contamination in the laboratory or the se- reads was 0.0046 ± 0.0056. quencing. These contaminants were species sequenced in different The absence of the mitogenome of the focal species in the ref- libraries of the same Illumina run and the fly Ceratitis capitata that erence database did not produce many false positives in the single-­ contaminated the library of Bactrocera oleae (see the discussion for species libraries (Table 3). When there were no congeneric species the reason behind this contamination). of the focal species in the reference database (six out of 17 spe- Because the MMG method must recover only one species in cies), no read was assigned to any species; when there were con- single-­species libraries, we fixed the values of γ = 0.99 and δ = 0.96 generic species in the database, in six cases no reads were assigned in the γ–­δ algorithm as this combination was the only one to pro- to any species and in four cases some reads were assigned to an- vide the expected result (Table S3). All the other tested combina- other species of the same genus (Bactrocera, Drosophila, Plutella, and tions reported false positives, such as Bactrocera biguttula in libraries Solenopsis); only in one species (Drosophila melanogaster) appeared of B. oleae, Drosophila formosana in libraries of D. melanogaster and some reads belonging to species of a different genus (Exorista sorbil- Solenopsis richteri in libraries of S. invicta. Results of all tested γ, δ lans, Diptera:Tachinidae). and input data combinations (raw reads versus candidate mito-­reads The six mixed-­ species libraries (Table S1) generated and Mito1794 versus FilteredMito1794 references) are provided as 3,376,087 ± 424,238 paired-­end reads. A proportion of Supporting Information (Tables S4–­S7). 0.003 ± 0.001 reads were eliminated in the trimming step and a pro- Filtering out the repetitive regions of the mitogenomes (i.e., portion of 0.925 ± 0.006 in the mito-­reads mining step. Therefore, FilteredMito1794) had a dramatic effect on the number of identified only a proportion of 0.075 ± 0.006 of the raw reads were candidate species. With the FilteredMito1794 database, we only detected the mito-­reads retained for further analysis. TA B L E 2 Summary of the results per library (mean ± SD) on the training data set of single-­species libraries for the four combinations of input data assessed in this study (raw reads and candidate mito-­reads mapped to Mito1794 and FilteredMito1794 databases) and using γ = 0.99 and δ = 0.96. (a) Number of recovered species per library, (b) Relative proportion of informative reads per library, and (c) processing time per library (format h:mm:ss; the time necessary to find the candidate mito-­reads is included in the processing time). Reads from contaminant species have not been considered Reference database Metric Input data Mito1794 FilteredMito1794 (a) Number of identified species Raw reads 2.52 ± 1.21 1±0 Mito-­reads 1.71 ± 0.90 1±0 (b) Relative proportion of informative reads Raw reads 0.0049 ± 0.0057 0.0047 ± 0.0056 Mito-­reads 0.0049 ± 0.0057 0.0046 ± 0.0057 (c) Processing time Raw reads 7:19:43 ± 3:46:53 7:21:10 ± 3:53:34 Mito-­reads 0:25:13 ± 0:17:07 0:24:37 ± 0:16:37 GARRIDO-­SANZ et al. | 7 estimation was very poor (r ≤ .67 and p > .05 for all samples; Figure 3a). Thus, it seems clear the need for a species-­specific cor- rection factor that considers a variable ratio of mitochondrial to nuclear DNA (Table 4). When we modified the actual RSA with the NMi correction factor (Equation 2), the correlation between actual and estimated RSA across-­species became significant in all samples (r ≥ .84 and p < .05 for all libraries; Figure 3b). The use of the NMi correction factor (Equation 4) instead of NMi provided an even better quantitative estimation of RSA across-­species (r ≥ .91 and p < .005 for all libraries; Figure 3c). The use of rarefacted samples showed that in the libraries with a more variable species abundance (libraries no. 1 and no. 2), the use of just half of the total available reads reduced the number of the identified species (Figure 4a) and promoted the presence of low-­abundant false positives, like Bactrocera biguttula above the de- tection limit ε = 0.0001 in library no. 2. On the contrary, when the abundance of species was less variable (libraries no. 5 and no. 6), the expected number of species was obtained with half the reads (Figure 4f) or even with 10% of reads (Figure 4e). 3.3 | Computer use The total consumed time by running the entire pipeline with the mixed-­species libraries ranged between 66 and 82 min (Table S11). Most of the processing time was devoted to the mapping of the reads to the references (95%) and only 3% of the time was used by the γ–­δ algorithm (Table S11). 4 | DISCUSSION Mito-­metagenomics (MMG) proved to be able to set apart and quan- tify the relative DNA abundance of insect species in artificial mix- tures, even when the species were congeneric. The estimation was F I G U R E 1 MMG pipeline applied in this study. In brackets, the as good as the one obtained with whole genomes instead of mitog- tools used in each step enomes by Garrido-­Sanz et al. (2020), using the same DNA libraries and bioinformatic methods as here. However, to be able to quantify In the mixed-­species libraries, we recovered all species included the RSA in a sample with several species (across-­species RSA) it was in the libraries except Papilio machaon in library no. 2 (Table S10). necessary to correct the raw reads by the variable amount of mito- As in the single-­species libraries, in the mixed-­species libraries, we chondrial to nuclear DNA (mitochondrial DNA copy number) among found reads of Ceratitis capitata and discarded them as laboratory species. contamination. In addition, in libraries no. 2 and no. 3, one single read was attributed to Bactrocera biguttula (Table S10), a species not handled in the laboratory; therefore, an analytical detection limit of 4.1 | Species identification ε = 0.0001 would be useful for the elimination of all false positives. The MMG approach used here recovered only the focal species from all the single-­species libraries, with no false positives when 3.2 | Estimation of the relative species abundance low-­complexity regions were filtered out from the mitogenomes in mixed-­species libraries (except for genuine contaminants, see below). Without this filtering step, some reads were attributed to non-­focal species; these reads The within-­species RSA was well estimated for all species (r ≥ .97 were sequences with biased composition, probably from repeti- and p < .05 for all species; Figure 2), but the across-­species RSA tive regions (e.g., microsatellites) that mostly matched non-­coding 8 | GARRIDO-­SANZ et al. TA B L E 3 List of species detected on the single-­species libraries when the mitogenome of the focal species is in the reference database (column A) and when it is not (column B). For each detected species we indicate its name and the number of assigned reads (in brackets). The number of congeneric species of the focal species included in the database is provided in column C. Libraries are divided into four groups: Group 1, species without congeneric species in the database and without false positive species; Group 2, species with congeneric species in the database and without false positive species; Group 3, species with congeneric species in the database but with false positive of the same genus; and group 4, species with congeneric species in the database but with false positive of a different genus Species used to (A) Focal species mitogenome (B) Focal species mitogenome (C) Number of congeneric Group Run Library prepare the library present in database not present in database species in the database 1 1 9 Acyrthosiphon pisum Acyrthosiphon pisum (154) None 0 2 1 Atta colombica Atta colombica (19412) None 0 2 3 Cimex lectularius Cimex lectularius (1082) None 0 1 6 Linepithema humile Linepithema humile (218) None 0 2 8 Linepithema humile Linepithema humile (913) None 0 2 11 Vollenhovia emeryi Vollenhovia emeryi (1399) None 0 2 12 Wasmannia Wasmannia auropunctata (17) None 0 auropunctata 2 1 8 Apis mellifera Apis mellifera (3597) None 7 2 2 Bemisia tabaci Bemisia tabaci (349) None 1 1 7 Bombus terrestris Bombus terrestris (1877) None 2 1 4 Drosophila Drosophila mojavensis (4407) None 18 mojavensis 2 5 Drosophila Drosophila mojavensis (2646) None 18 mojavensis 2 7 Drosophila suzukii Drosophila suzukii (4664) None 18 1 1 Papilio machaon Papilio machaon (312) None 13 3 1 5 Bactrocera oleae Bactrocera oleae (425) Bactrocera biguttula (137) 14 1 2 Drosophila virilis Drosophila virilis (8070) Drosophila littoralis (20) 18 2 6 Drosophila virilis Drosophila virilis (2704) Drosophila littoralis (6) 18 2 9 Plutella xylostella Plutella xylostella (2371) Plutella australiana (127) 1 2 10 Solenopsis invicta Solenopsis invicta (131) Solenopsis richteri (236) 2 4 1 3 Drosophila Drosophila melanogaster (979) Drosophila formosana (312) 18 melanogaster Exorista sorbillans (26) Drosophila mauritiana (3) 2 4 Drosophila Drosophila melanogaster (1384) Drosophila formosana (478) 18 melanogaster Exorista sorbillans (53) Drosophila mauritiana (3) regions on the reference mitogenomes (Table S8; Faber & Stepien, captured together in fly traps. A more throughout discussion about 1998; Wolff et al., 2012). Some popular tools do implicitly or explic- this problem is provided in Garrido-­Sanz et al. (2020). The removal itly filter out low complexity regions from the reference genomes: of the genuine contaminant species in artificial libraries as we did Kraken masks low complexity regions when adding references to the here was possible because we knew the identity of the species in the database (Wood & Salzberg, 2014); and BLAST filters both query se- mixture, but it is impossible in real samples. quences and references (Altschul et al., 1990, 1997; Camacho et al., 2009). We also found contaminant species in all the sequenced librar- 4.2 | Quantification of the RSA and the need for a ies (Tables S9–­S10). The origin of these reads of contaminants can species-­specific correction factor be tag jumping during the sequencing reaction (Schnell et al., 2015) or actual contamination in the laboratory. The first reason is prob- With the mixed-­species libraries, we estimated the within-­species ably the cause of finding reads in single-­species libraries belonging RSA with high statistical confidence (Figure 2). Similar good results to other species sequenced in the same run but different libraries. have been reported in previous studies that used mock samples The second reason is behind the presence of reads attributed to (Bista et al., 2018; Ji et al., 2020). On the contrary, the across-­species Ceratitis capitata in libraries where Bactrocera oleae was also pres- RSA estimation within a sample was not statistically significant in ent, because the two dipterans, which are agricultural pests, were any sample (Figure 3a). These results contrast with the study of GARRIDO-­SANZ et al. | 9 F I G U R E 2 Scatter plot of the 0.55 r = 0.97 , p < 0.005 estimated versus the actual RSA for each 0.0075 r = 0.99 , p < 0.005 species of the mixed-­species libraries (i.e., 0.50 within-­species RSA). Each plot shows the 0.0050 Pearson correlation coefficient (r) and the 0.45 corresponding p-­value. The coordinate at 0.0025 the origin of all regression lines was not different to 0 0.40 0.0000 A. pisum A. mellifera 0.00 0.03 0.06 0.09 0.14 0.16 0.18 r = 1 , p < 0.005 0.20 r = 0.97 , p < 0.005 0.3 0.15 Estimated relative concentration 0.2 0.10 B. oleae B. terrestris 0.1 0.05 0.3 0.4 0.5 0.6 0.10 0.15 0.20 0.25 0.035 r = 0.99 , p < 0.005 0.25 r = 0.99 , p < 0.005 0.030 0.20 0.025 0.15 0.020 0.10 0.015 D. melanogaster 0.05 D. mojavensis 0.010 0.03 0.06 0.09 0.12 0.05 0.10 0.006 r = 0.99 , p = 0.015 r = 0.97 , p < 0.005 0.005 0.004 0.004 0.003 0.003 0.002 0.002 0.001 0.001 L. humile 0.000 P. machaon 0.01 0.02 0.03 0.04 0.05 0.00 0.03 0.06 0.09 Actual relative concentration Gueuning et al. (2019) that reported good quantitative estimations one reported here, as the result was almost the same in both cases. of across-­species RSA in artificial mixtures of wild bees. One possible explanation is that Tang et al. (2015) dealt only with Using the mitochondrial DNA copy number correction factor, the wild bees (a group of a few Hymenoptera families), so the interspe- RSA across-­species correlated significantly with the real values in the cific differences might be low compared to our study which included six artificial samples analysed (Equation 2; Figure 3b), even when species from four insect orders. Nevertheless, the effect of the mi- the mean genome and mitogenome sizes were used instead of the tochondrial DNA copy number on the estimation of RSA deserves species-­specific value (Equation 4; Figure 3c). Other studies report- more research effort if DNA-­based techniques are to provide good ing RSA estimation across-­species also used some correction factor quantitative results, both using mito-­metagenomics and amplicon before comparing the expected and observed number of reads, but metabarcoding targeting variable copy number regions. none included the genome size, as we did here. For instance, Gómez-­ The method used here to estimate the correction factor for the Rodríguez et al. (2015) and Tang et al. (2015) considered the mitoge- mitochondrial DNA copy number (i.e., preparation and sequencing nome size of the species and Tang et al. (2015) and Lang et al. (2019) of a single-­species library to a depth of c. one million reads) has a the number of reads from the genome. Tang et al. (2015) is the only cost that is not negligible. Ideally, there should be a method to in- study that provides both the goodness of fit with and without the dependently estimate the mitochondrial DNA copy number of each use of the correction factor, and the effect is very different from the species that did not involve sequencing. Such methods do exist 10 | GARRIDO-­SANZ et al. (a) (b) (c) r = 0.67 , p = 0.067 r = 0.93 , p < 0.005 0.4 r = 0.95 , p < 0.005 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0.0 Lib. 1 0.0 Lib. 1 0.0 Lib. 1 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.5 0.5 r = 0.65 , p = 0.12 0.5 r = 0.97 , p < 0.005 r = 0.99 , p < 0.005 0.4 0.4 0.4 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.0 Lib. 2 0.0 Lib. 2 0.0 Lib. 2 0.0 0.2 0.4 0.6 0.0 0.1 0.2 0.3 0.4 0.0 0.2 0.4 0.6 0.5 r = 0.48 , p = 0.23 r = 0.91 , p < 0.005 r = 0.95 , p < 0.005 0.4 0.4 0.4 Estimated relative concentration 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0.0 Lib. 3 0.0 Lib. 3 0.0 Lib. 3 0.1 0.2 0.3 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.2 0.4 r = 0.52 , p = 0.19 r = 0.97 , p < 0.005 r = 0.98 , p < 0.005 0.4 0.4 0.4 0.2 0.2 0.2 0.0 Lib. 4 0.0 Lib. 4 0.0 Lib. 4 0.1 0.2 0.3 0.4 0.0 0.2 0.4 0.0 0.2 0.4 0.6 r = 0.38 , p = 0.4 r = 0.84 , p = 0.017 r = 0.91 , p < 0.005 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0.0 Lib. 5 0.0 Lib. 5 0.0 Lib. 5 0.08 0.12 0.16 0.20 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.5 r = 0.44 , p = 0.33 0.5 r = 0.92 , p < 0.005 0.5 r = 0.95 , p < 0.005 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0.0 Lib. 6 0.0 Lib. 6 0.0 Lib. 6 0.10 0.15 0.20 0.25 0.1 0.2 0.3 0.4 0.5 0.0 0.2 0.4 0.6 Actual relative concentration GARRIDO-­SANZ et al. | 11 F I G U R E 3 Scatter plot of the estimated versus the actual RSA in each mixed-­species library (i.e., across-­species RSA). At the top, it is indicated the way we conducted the actual RSA. (a) Original expected data. (b) Corrected expected data after applying the NMi correction factor. (c) Corrected expected data after applying the NMi correction factor. Rows from top to bottom correspond to mixed-­species libraries from 1 to 6. Each plot shows the Pearson correlation coefficient (r) and the corresponding p-­value because there is an interest in medicine to measure the mitochon- 4.3.1 | More complex mixtures drial DNA copy number for its relationship with several diseases. In medicine, the mitochondrial DNA copy number is usually esti- Real samples can contain hundreds of species, and that might affect mated using quantitative PCR (qPCR; Thyagarajan et al., 2012); this the ability of the method to detect the less abundant ones. With the method requires two primer pairs, one for a mitochondrial marker sequencing depth achieved here (~3.4 million raw reads per sample; and one for a single-­copy nuclear marker. These primers are known Table S1) we were able to detect three (out of four) species with an for humans, but it would be costly to generate them for every spe- expected RSA below 1‰ (Table S10). The subsequent rarefaction cies in an environmental sample, especially for those species whose experiment showed that with fewer reads more species become un- genome has not yet been sequenced (but see Liu et al., 2018). The detected (Figure 4). In consequence, it seems that at least 3.4·106 qPCR itself is cheap, but the preparative work for each species reads are needed to detect most species with an RSA above 1‰. would be long. Having hundreds of species in the mixture would not hamper the It is important to emphasize that the RSA used here is based on quantitative ability of the method, as most species would be above a the relative proportions of DNA of the species in the mixture, but 1‰ abundance. However, ultra-­rich samples with thousands of spe- what is needed in most ecological applications is the relative propor- cies would require a higher sequencing depth to ensure the detec- tions of biomass (or individual counts) of the different species. The tion of most species. reason behind our choice was to simplify the problem of obtaining the actual RSA from high-­throughput sequenced reads in several steps. There is one bias caused by the variable mitochondrial copy 4.3.2 | Limited amount of DNA available number of the different species and there is another, independent, bias caused by the variable DNA content of the biomass of different The Illumina TruSeq kit used here to prepare the libraries requires species. We addressed here the first bias and obviated the second 1 μg of DNA and this might be a problem with small specimens or in one. From our results, the RSA based on the biomass could be ob- DNA-­poor samples. However, today there are alternative methods tained by multiplying our estimates by the biomass-­to-­DNA ratio of that provide good results with just 1 ng of DNA, like the Illumina each species, if known. There are very little data in the literature Nextera DNA Flex kit (Sato et al., 2019), albeit potential biases about the proportion of DNA to biomass in different species, but it should be tested in future experiments for these kits. can be very variable; for instance, Pornon et al. (2016) reports a very different DNA yield from the same number of pollen grains of three plant species. 4.3.3 | Absence of mitogenomes in the In metabarcoding applications, some authors use empirical reference database correction factors based on mixtures of known relative biomass of several species rather than in mixtures of DNA (e.g., Matesanz Our results showed that the proposed methodology was robust in et al., 2019; Thomas et al., 2016). In these cases, the correction fac- the absence of the mitogenomes of species in the reference data- tor solves for the mitochondrial copy number among species and base. Of course, the species without their mitogenome in the refer- also for the DNA-­to-­biomass ratio and the differential amplification ence database will never be found, but their reads will not generate efficiency caused by PCR. This method is undoubtedly practical many false positives, even for species with close relatives in the but does not differentiate the relative importance of each source reference database (Table 3). The presence of species without their of bias. reference genome in the mixture is likely to occur frequently in real samples. The unassigned reads (or also when the prior composi- tion of the sample is unknown) can be further explored by mapping 4.3 | Mito-­metagenomics in real samples them against other databases, like COI barcodes from BOLD; thus, the identity of more species will be revealed, albeit not their relative The present study is based on artificial mixtures of a low num- abundance. ber of species whose mitogenomes are already assembled. Thus, it is fair to question the value of our proposal in real samples with many more species, with a limited amount of DNA, where the prior 4.3.4 | Incomplete genomes species composition is unknown, or when the reference mitoge- nomes are obtained in the same experiment and are only partially In several MMG studies, the reference mitogenomes are assembled assembled. from the same mixtures in which the RSA is intended to be quantified 12 | GARRIDO-­SANZ et al. species i divided by the total number of reads from that single-­species library; NMi is the estimated number of mitochondrial DNA copies for species i; NMi is as NMi but calculated using the mean TA B L E 4 Summary table of the species in single-­species libraries data used for the obtention of the correction factors NMi and NMi. xi is the number of reads mapping into the mitogenome of (Crampton-­Platt et al., 2016; Zhou et al., 2013). In these cases, the Equation 4) NMi (from mitogenomes are assembled de novo and, normally, they are incom- 95.5 287.8 217.6 145.6 1115.6 46.0 170.4 220.9 187.1 plete. We do not see any impediment in using mitogenomes as- sembled in this way if all of them have a similar length and quality. However, we would advise against the simultaneous use of mitog- Equation 2) enomes with disparate length or quality for quantification purposes, because that would bias the RSA towards the species with better NMi (from mitogenomes (Tang et al., 2015). On the contrary, the use of all avail- 180.5 209.5 30.3 144.6 84.4 66.4 117.4 823.9 340.1 able partial mitogenomes would be fine for identification purposes. xi (from Equation 4.3.5 | Estimation of the mitochondrial DNA copy number (NMi) 0.0046 0.0089 0.0082 0.0022 0.0070 0.0104 0.0105 0.0137 0.0511 3) Perhaps the most difficult problem in real samples is the estima- tion of NMi. The rationale that we propose for the estimation of NMi Number of reads mapped to (Equations 1–­4) seems reasonable, but the devil is in the detail: the estimation of the variables needed to calculate NMi is paved with dif- ficulties for species without a reference genome. First, the estima- tion of the proportion of reads that belong to the mitogenomes (xi, the mitogenome Equation 3) is biased, because we assumed that all reads belong to the same species (RGi, Equation 3); however, there is always DNA 64,448 22,839 16,262 7,218 4,052 1,992 20,424 48,967 3,151 that comes from other sources, like food, gut bacteria, parasites, etc. Consequently, the number of reads of the entire genome RGi is overestimated and, hence, xi underestimated. Second, the size of the mitogenome and the whole genome is generally unknown for most Number of raw reads species; even for the best studied species their whole genome is far from complete (e.g., Paris et al., 2020), and the estimated size (Gi ) is 1,262,388 1,422,342 1,668,424 2,283,768 1,994,938 4,714,902 an underestimation of the real size. (paired-­end) 684,688 434,520 580,996 Nevertheless, despite the above problems, the correction factor NMi helped to reproduce the expected across-­species RSA in our li- braries. Similar results were obtained using the mean values of the mitogenome and whole genome sizes (NMi, Equation 4). The apparent length of whole genomes and mitogenomes of all species considered here Mitochondrial genome lack of effect of the species-­specific genome and mitogenome sizes might be caused to the low variability of the mitogenome size (coef- ficient of variation, CV = 9%) and moderate variability of the whole genome size (CV = 47%; Table 4). On the contrary, the proportion of size (kbp) reads mapping into the mitogenome (xi, Equation 3) was much more 15.2 15.8 16.3 17.0 20.0 17.4 14.9 14.9 16.1 variable among species (CV = 113%). Given the previous considerations, we suggest the use of the Complete genome size correction factor NMi instead of NMi, for species without a reference genome and to estimate the three necessary variables (xi, Mi, Gi ) in the following way. • xi. The proportion of reads belonging to the mitogenome of spe- (Mbp) 278 206 144 472 197 220 249 250 542 cies i could be estimated by shotgun sequencing a single-­species DNA extract. The value of xi would be an underestimation of the Drosophila mojavensis real value but given the high interspecific variability, the obtained Acyrthosiphon pisum Linepithema humile Bombus terrestris xi values should still be useful for correction purposes. melanogaster Bactrocera oleae Drosophila virilis Papilio machaon Apis mellifera • Mi. Ninety per cent of the 1794 mitogenomes used here have a Drosophila length of 14.9 to 17.0 kpb (i.e., a rank of 2.1 kbp or 13% of the Species mean Mi; Table S2). Consequently, we recommend the use of the mean value M for the group of species of interest (Table S12). GARRIDO-­SANZ et al. | 13 F I G U R E 4 Number of identified 8 8 species using different proportions of reads in the mixed-­species libraries. Each 7 7 simulation was performed 100 times with different subsets, except when the 6 6 entire library was used. Letters from a to f 5 5 indicate mixed-­species libraries from 1 to 6. Grey dashed lines indicate the expected 4 (a) 4 (b) number of recovered species in each 1 0.5 0.1 0.05 0.01 1 0.5 0.1 0.05 0.01 library Number of identified species 9 8 8 7 7 6 6 5 5 4 (c) 4 (d) 1 0.5 0.1 0.05 0.01 1 0.5 0.1 0.05 0.01 7 7 6 6 5 (f) 5 (e) 4 1 0.5 0.1 0.05 0.01 1 0.5 0.1 0.05 0.01 Relative proportion of initial sample • Gi. The length of the whole genome is more variable across spe- and Mario Cáceres for his advice on a preliminary version of this cies than the length of the mitogenomes: 90% of the 115 whole manuscript. Financial support was provided by Spanish Government genomes of insects available at RefSeq (Table S13A) have a length grant TIN2017-­8 4553-­C2-­1-­R and Generalitat de Catalunya grant between 0.14 and 0.98 Gbp. However, if the insects are split by AGAUR 2017-­SGR-­1001. orders the variability of Gi is smaller for most insect orders (Table S13B). Consequently, we would advise using the mean value Ḡ for AU T H O R C O N T R I B U T I O N S each group of taxa (e.g., insect orders). Lidia Garrido-­Sanz and Josep Piñol conceived the experiment. Lidia Garrido-­Sanz performed the bioinformatic analyses and the other two authors supervised it. All authors were involved in the analysis, 4.4 | Concluding remarks interpretation of data and writing of the manuscript. The approach presented here to identify insect species and to es- DATA AVA I L A B I L I T Y S TAT E M E N T timate their relative abundance in complex mixtures using MMG Sequenced DNA samples have been made available at the Dryad worked well with artificial samples of known composition for a select repository: https://doi.org/10.5061/dryad.t1g1j​wsz7. The γ–­δ algo- group of species whose mitogenomes are sequenced to an advanced rithm script is available at GitHub: https://github.com/Lidia​GS/g-­d_ degree. The key for the accurate estimation of the across-­species algor​ithm. RSA was a correction factor for the mitochondrial copy number of each species. We are aware that the proposed methodology is not ORCID immediately applicable to most real samples, so its real value should Lidia Garrido-­Sanz https://orcid.org/0000-0003-2622-9674 be tested on more of such samples. Miquel Àngel Senar https://orcid.org/0000-0002-0316-5420 Josep Piñol https://orcid.org/0000-0002-4067-3301 AC K N OW L E D G E M E N T S We are grateful to several entomologists who provided the speci- REFERENCES mens: Xavier Espadaler, Nicolás Pérez, Alfredo Ruíz, Francesc Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Mestres, Aleix Valls, Francisco Beitia, Carlos Hernández-­C astellano, Basic local alignment search tool. Journal of Molecular Biology, Joan Josep Ibañez, Carlos Pradera, Rasmus S. Larsen, Jacobus J. 215(3), 403–­410. https://doi.org/10.1016/s0022​-­2836(05)80360​-­2 Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., Miller, W., Boomsma, Luis Calcaterra and Misato O. Miyakawa. We also thank & Lipman, D. J. (1997). Gapped BLAST and PSI-­BLAST: A new gen- Anna Barceló and Roger Lahoz of the Genomics facilities of the UAB eration of protein database search programs. Nucleic Acids Research, for the preparation of the DNA libraries and sequencing the DNA 25(17), 3389–­3402. https://doi.org/10.1093/nar/25.17.3389 14 | GARRIDO-­SANZ et al. Andrews, S. (2015). FastQC: A quality control tool for high throughput se- Krehenwinkel, H., Wolf, M., Lim, J. Y., Rominger, A. J., Simison, W. B., & quence data. http://www.bioin​forma​tics.babra​ham.ac.uk/proje​c ts/ Gillespie, R. G. (2017). Estimating and mitigating amplification bias fastqc in qualitative and quantitative arthropod metabarcoding. Scientific Bista, I., Carvalho, G. R., Tang, M., Walsh, K., Zhou, X., Hajibabaei, M., Reports, 7, 17668. https://doi.org/10.1038/s4159​8-­017-­17333​-­x Shokralla, S., Seymour, M., Bradley, D., Liu, S., Christmas, M., & Lang, D., Tang, M., Hu, J., & Zhou, X. (2019). Genome-­skimming Creer, S. (2018). Performance of amplicon and shotgun sequenc- provides accurate quantification for pollen mixtures. ing for accurate biomass estimation in invertebrate community Molecular Ecology Resources, 19(6), 1433–­ 1446. https://doi. samples. Molecular Ecology Resources, 18, 1020–­1034. https://doi. org/10.1111/1755-­0998.13061 org/10.1111/1755-­0998.12888 Li, H. (2013). Aligning sequence reads, clone sequences and assembly Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: a flexible contigs with BWA-­MEM. arXiv, 1303.3997v1 [q-­bio.GN]. trimmer for Illumina sequence data. Bioinformatics, 30(15), 2114–­ Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, 2120. https://doi.org/10.1093/bioin​forma​tics/btu170 G., Abecasis, G., & Durbin, R. (2009). The sequence alignment/map Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, format and SAMtools. Bioinformatics, 25(16), 2078–­2079. https:// J., Bealer, K., & Madden, T. L. (2009). BLAST+: Architecture doi.org/10.1093/bioin​forma​tics/btp352 and applications. BMC Bioinformatics, 10, 421. https://doi. Liu, R., Jin, L., Long, K., Tang, Q., Ma, J., Wang, X., Zhu, L. I., Jiang, A., org/10.1186/1471-­2105-­10-­421 Tang, G., Jiang, Y., Li, X., & Li, M. (2018). Analysis of mitochondrial CBoL Plant Working Group (2009). A DNA barcode for land plants. DNA sequence and copy number variation across five high-­altitude Proceedings of the National Academy of Sciences of the United States species and their low-­altitude relatives. Mitochondrial DNA Part B, of America, 106(31), 12794–­12797. https://doi.org/10.1073/ 3(2), 847–­851. https://doi.org/10.1080/23802​359.2018.1501285 pnas.09058​45106 Lu, J., & Salzberg, S. L. (2018). Removing contaminants from databases Crampton-­Platt, A., Timmermans, M. J. T. N., Gimmel, M. L., Kutty, S. N., of draft genomes. PLoS Computational Biology, 14(6), e1006277. Cockerill, T. D., Khen, C. V., & Vogler, A. P. (2015). Soup to tree: The https://doi.org/10.1371/journ​al.pcbi.1006277 phylogeny of beetles inferred by mitochondrial metagenomics of Matesanz, S., Pescador, D. S., Pías, B., Sánchez, A. M., Chacón-­L abella, J., a borean rainforest sample. Molecular Biology and Evolution, 32(9), Illuminati, A., Cruz, M., López-­Angulo, J., Marí-­Mena, N., Vizcaíno, 2302–­2316. https://doi.org/10.1093/molbe​v/msv111 A., & Escudero, A. (2019). Estimating belowground plant abun- Crampton-­Platt, A., Yu, D. W., Zhou, X., & Vogler, A. P. (2016). dance with DNA metabarcoding. Molecular Ecology Resources, 19(5), Mitochondrial metagenomics: letting the genes out of the bottle. 1265–­1277. https://doi.org/10.1111/1755-­0998.13049 Gigascience, 5(1), 15. https://doi.org/10.1186/s1374​2-­016-­0120-­y Morgulis, A., Gertz, E. M., Schäffer, A. A., & Agarwala, R. (2006). A fast Elbrecht, V., & Leese, F. (2015). Can DNA-­based ecosystem assessments and symmetric DUST implementation to mask low-­ complexity quantify species abundance? Testing primer bias and biomass-­ DNA sequences. Journal of Computational Biology, 13(5), 1028–­ sequence relationships with an innovative metabarcoding pro- 1040. https://doi.org/10.1089/cmb.2006.13.1028 tocol. PLoS One, 10(7), e0130324. https://doi.org/10.1371/journ​ Paris, M., Boyer, R., Jaenichen, R., Wolf, J., Karageorgi, M., Green, J., al.pone.0130324 Cagnon, M., Parinello, H., Estoup, A., Gautier, M., Gompel, N., & Faber, J. E., & Stepien, C. A. (1998). Tandemly repeated sequences in the Prud’homme, B. (2020). Near-­chromosome level genome assem- mitochondrial DNA control region and phylogeography of the pike-­ bly of the fruit pest Drosophila suzukii using long-­read sequencing. perches Stizostedion. Molecular Phylogenetics and Evolution, 10(3), Scientific Reports, 10, 11227. https://doi.org/10.1038/s4159​8-­020-­ 310–­322. https://doi.org/10.1006/mpev.1998.0530 67373​-­z Garrido-­Sanz, L., Senar, M. À., & Piñol, J. (2020). Estimation of the rela- Pearman, W. S., Freed, N. E., & Silander, O. K. (2019). The advantages tive abundance of species in artificial mixtures of insects using low-­ and disadvantages of short-­and long-­read metagenomics to infer coverage shotgun metagenomics. Metabarcoding and Metagenomics, bacterial and eukaryotic community composition. bioRxiv. https:// 4, e48281. https://doi.org/10.3897/mbmg.4.48281 doi.org/10.1101/650788 Gómez-­Rodríguez, C., Crampton-­Platt, A., Timmermans, M. J. T. N., Piñol, J., Mir, G., Gomez-­Polo, P., & Agustí, N. (2015). Universal and Baselga, A., & Vogler, A. P. (2015). Validating the power of mito- blocking primer mismatches limit the use of high-­ throughput chondrial metagenomics for community ecology and phylogenet- DNA sequencing for the quantitative metabarcoding of arthro- ics of complex assemblages. Methods in Ecology and Evolution, 6(8), pods. Molecular Ecology Resources, 15(4), 819–­830. https://doi. 883–­894. https://doi.org/10.1111/2041-­210X.12376 org/10.1111/1755-­0998.12355 Gregory, T. R. (2020). Animal genome size database. http://www.genom​ Pornon, A., Escaravage, N., Burrus, M., Holota, H., Khimoun, A., Mariette, esize.com J., Pellizzari, C., Iribar, A., Etienne, R., Taberlet, P., Vidal, M., Gueuning, M., Ganser, D., Blaser, S., Albrecht, M., Knop, E., Praz, C., & Winterton, P., Zinger, L., & Andalo, C. (2016). Using metabarcod- Frey, J. E. (2019). Evaluating next-­generation sequencing (NGS) ing to reveal and quantify plant-­pollinator interactions. Scientific methods for routine monitoring of wild bees: Metabarcoding, mi- Reports, 6, 27282. https://doi.org/10.1038/srep2​7282 togenomics or NGS barcoding. Molecular Ecology Resources, 19(4), R Core Team (2016). R: A language and environment for statistical comput- 847–­862. https://doi.org/10.1111/1755-­0998.13013 ing. R Foundation for Statistical Computing. https://www.R-­proje​ Hebert, P. D. N., Ratnasingham, S., & de Waard, J. R. (2003). Barcoding ct.org/ animal life: cytochrome c oxidase subunit 1 divergences among RStudio Team (2015). RStudio: Integrated development for R. RStudio Inc. closely related species. Proceedings of the Royal Society of London. http://ww.rstud​io.com/ Series B: Biological Sciences, 270, S96–­S99. https://doi.org/10.1098/ Saitoh, S., Aoyama, H., Fujii, S., Sunagawa, H., Nagahama, H., Akutsu, M., rsbl.2003.0025 Shinzato, N., Kaneko, N., & Nakamori, T. (2016). A quantitative pro- Ji, Y., Huotari, T., Roslin, T., Schmidt, N. M., Wang, J., Yu, D. W., & tocol for DNA metabarcoding of springtails (Collembola). Genome, Ovaskainen, O. (2020). SPIKEPIPE: A metagenomic pipeline for 59(9), 705–­723. https://doi.org/10.1139/gen-­2015-­0228 the accurate quantification of eukaryotic species occurrences Sato, M. P., Ogura, Y., Nakamura, K., Nishida, R., Gotoh, Y., Hayashi, and intraspecific abundance change using DNA barcodes or mi- M., Hisatsune, J., Sugai, M., Takehiko, I., & Hayashi, T. (2019). togenomes. Molecular Ecology Resources, 20, 256–­267. https://doi. Comparison of the sequencing bias of currently available library org/10.1111/1755-­0998.13057 preparation kits for Illumina sequencing of bacterial genomes Kassambara, A. (2018). ggpubr: ‘ggplot2’ Based publication ready plots. R and metagenomes. DNA Research, 26(5), 391–­398. https://doi. package version 0.2. https://CRAN.R-­proje​c t.org/packa​ge=ggpubr org/10.1093/dnare​s/dsz017 GARRIDO-­SANZ et al. | 15 Schnell, I. B., Bohmann, K., & Gilbert, M. T. P. (2015). Tag jumps illumi- risk. Cancer Epidemiology, Biomarkers & Prevention, 21(9), 1574–­ nated -­reducing sequence-­ to-­sample misidentifications in me- 1581. https://doi.org/10.1158/1055-­9965.EPI-­12-­0138-­T tabarcoding studies. Molecular Ecology Resources, 15(6), 1289–­ Wickham, H. (2016). ggplot2: Elegant graphics for data analysis. Springer-­ 1303. https://doi.org/10.1111/1755-­0998.12402 Verlag. http://ggplo​t 2.org Schoch, C. L., Seifert, K. A., Huhndorf, S., Robert, V., Spouge, J. L., Wolff, J. N., Shearman, D. C. A., Brooks, R. C., & Ballard, J. W. O. Levesque, C. A., Chen, W., Bolchacova, E., Voigt, K., Crous, P. (2012). Selective enrichment and sequencing of whole mitochon- W., Miller, A. N., Wingfield, M. J., Aime, M. C., An, K.-­D., Bai, F.-­ drial genomes in the presence of nuclear encoded mitochon- Y., Barreto, R. W., Begerow, D., Bergeron, M.-­J., Blackwell, M., drial pseudogenes (Numts). PLoS One, 7(5), e37142. https://doi. … Schindel, D. (2012). Nuclear ribosomal internal transcribed org/10.1371/journ​al.pone.0037142 spacer (ITS) region as a universal DNA barcode marker for Fungi. Wood, D. E., & Salzberg, S. L. (2014). Kraken: ultrafast metagenomic se- Proceedings of the National Academy of Sciences of the United quence classification using exact alignments. Genome Biology, 15, States of America, 109(16), 6241–­6246. https://doi.org/10.1073/ R46. https://doi.org/10.1186/gb-­2014-­15-­3-­r46 pnas.11170​18109 Zhou, X., Li, Y., Liu, S., Yang, Q., Su, X. U., Zhou, L., Tang, M., Fu, Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C., & Willerslev, E. R., Li, J., & Huang, Q. (2013). Ultra-­deep sequencing enables (2012). Towards next-­ generation biodiversity assessment using high-­fidelity recovery of biodiversity for bulk arthropod sam- DNA metabarcoding. Molecular Ecology, 21(8), 2045–­2050. https:// ples without PCR amplification. Gigascience, 2(1), 4. https://doi. doi.org/10.1111/j.1365-­294X.2012.05470.x org/10.1186/2047-­217X-­2-­4 Tang, M., Hardman, C. J., Ji, Y., Meng, G., Liu, S., Tan, M., Yang, S., Moss, E. D., Wang, J., Yang, C., Bruce, C., Nevard, T., Potts, S. G., Zhou, X., & Yu, D. W. (2015). High-­throughput monitoring of wild bee diversity S U P P O R T I N G I N FO R M AT I O N and abundance via mitogenomics. Methods in Ecology and Evolution, Additional supporting information may be found online in the 6, 1034–­1043. https://doi.org/10.1111/2041-­210X.12416 Tang, M., Tan, M., Meng, G., Yang, S., Su, X. U., Liu, S., Song, W., Li, Y., Wu, Supporting Information section. Q., Zhang, A., & Zhou, X. (2014). Multiplex sequencing of pooled mitochondrial genomes-­a crucial step toward biodiversity analysis using mito-­metagenomics. Nucleic Acids Research, 42(22), e166. How to cite this article: Garrido-­Sanz, L., Senar, M. À., & https://doi.org/10.1093/nar/gku917 Piñol, J. (2021). Relative species abundance estimation in Thomas, A. C., Deagle, B. E., Everson, J. P., Harsch, C. H., & Trites, A. artificial mixtures of insects using mito-­metagenomics and a W. (2016). Quantitative DNA metabarcoding: improved estimates of species proportional biomass using correction factors derived correction factor for the mitochondrial DNA copy number. from control material. Molecular Ecology Resources, 16(3), 714–­726. Molecular Ecology Resources, 00, 1–­15. https://doi. https://doi.org/10.1111/1755-­0998.12490 org/10.1111/1755-­0998.13464 Thyagarajan, B., Wang, R., Barcelo, H., Koh, W.-­P., & Yuan, J.-­M. (2012). Mitochondrial copy number is associated with colorectal cancer