Metazoan mitochondrial gene sequence reference datasets for taxonomic assignment of environmental samples

Mitochondrial-encoded genes are increasingly targeted in studies using high-throughput sequencing approaches for characterizing metazoan communities from environmental samples (e.g., plankton, meiofauna, filtered water). Yet, unlike nuclear ribosomal RNA markers, there is to date no high-quality reference dataset available for taxonomic assignments. Here, we retrieved all metazoan mitochondrial gene sequences from GenBank, and then quality filtered and formatted the datasets for taxonomic assignments using taxonomic assignment tools. The reference datasets—‘Midori references’—are available for download at www.reference-midori.info. Two versions are provided: (I) Midori-UNIQUE that contains all unique haplotypes associated with each species and (II) Midori-LONGEST that contains a single sequence, the longest, for each species. Overall, the mitochondrial Cytochrome oxidase subunit I gene was the most sequence-rich gene. However, sequences of the mitochondrial large ribosomal subunit RNA and Cytochrome b apoenzyme genes were observed for a large number of species in some phyla. The Midori reference is compatible with some taxonomic assignment software. Therefore, automated high-throughput sequence taxonomic assignments can be particularly effective using these datasets.


Background & Summary
Massively parallel sequencing technologies have revolutionized our ability to survey and monitor biological diversity. Samples containing multiple species are collected directly from the environment (e.g., from water, sediments, traps, feces) and variants of one or several sets of genes are inventoried using PCR-based (e.g., metagenetics) or PCR-free (e.g., mitogenomics, metatranscriptomics, metagenomics) approaches.
Two types of gene sequences have been widely used as phylogenetic and taxonomic markers in metazoans: nuclear-encoded ribosomal RNA genes (18S and 28S ribosomal RNA genes 1,2 ) and mitochondrial-encoded genes 3,4 . Nuclear-encoded ribosomal RNA fragments, especially hypervariable regions of the 18S rRNA gene, were prime targets in early metagenetic analyses because broad-range primers well conserved across the eukaryotic domain were available 2,5,6 . As a result, considerable efforts have already been made to build quality filtered and formatted reference sequence datasets of nuclearencoded ribosomal RNA genes for taxonomic assignments 7,8 . Mitochondrial genes, which provide higher taxonomic resolution for most metazoan groups 9 , have been increasingly used following the design of highly degenerate primer sets 4,10-12 and the development of bioinformatics tools to facilitate the assembly of mitogenomes from environmental samples 13 . However, high quality reference datasets that are compatible with taxonomic assignment software are not yet available for metazoan mitochondrial genes. Therefore, at the moment, most of the metazoan metagenetic studies target low-resolution nuclear ribosomal RNA gene as a marker (e.g., SILVA 7 ). Some exceptions, which target high-resolution mitochondrial genes, used Blastn searches against sequences from GenBank for taxonomic assignments 10,12 without explicit taxonomic quality control of the database. This means that highthroughput sequence taxonomic assignment with quality controlled mitochondrial gene reference dataset is currently not feasible. Here, we constructed quality-controlled reference datasets 'Midori' for thirteen protein-coding and two ribosomal RNA genes sequences encoded in the mitochondrial genome.
After downloading the nt datasets from GenBank, we implemented the following informatics procedures to construct the datasets ( Fig. 1): (I) extraction of mitochondria-related gene sequences, (II) extraction of metazoan gene sequences, (III) insertion of taxonomic ranking information onto each sequence, (IV) removal of sequences that did not have species-level taxonomic information, (V) splitting the sequences into thirteen protein and two ribosomal RNA gene sequences, (VI) preparation of a taxonomy rank file.
We have prepared two versions of the 15 datasets: (I) Midori-UNIQUE, which contains for each species a representative sequence of each unique haplotype and (II) Midori-LONGEST, which contains for each species the single longest sequence. Each dataset is composed of two ribosomal RNA (    . In both datasets, COI had the largest number of sequences overall (Table 1). However, Cytb had a higher number of sequences in Chordata, while the lrRNA gene had the highest coverage among Cnidaria, Hemichordata and Placozoa (Fig. 2, Table 2). We provide two formats of Midori-UNIQUE and Midori-LONGEST (www.reference-midori.info) that are compatible with taxonomic assignment software such as RDP Classifier 15 , SPINGO 16 and BLAST+ ref. 17. We believe these datasets will improve the accuracy of taxonomic assignments of the massive numbers of sequences obtained from high throughput sequencing experiments.

Mitochondrial DNA sequence curation
The nt fasta file was downloaded from the National Centre for Biotechnology Information (NCBI) server (ftp://ftp.ncbi.nih.gov/blast/db/FASTA) on 18 September 2015 (Fig. 1). Mitochondria-related gene sequences (including nuclear-encoded mitochondrial genes) were extracted from the fasta file using a custom Perl script (02_ext_seq_mito.pl) available online (www.reference-midori.info/download.php#) along with following scripts described below to build the reference datasets. Next, GenBank flat files of all the mitochondria-related gene sequences were downloaded using NCBI Edirect (efetch -db nucleotide -id gene_id -format gb), and metazoan flat files were extracted using a custom perl script (06_ext_fasta_seq.pl). Next, CDS and rRNA features were extracted from the metazoan flat files using a custom Perl script (09_ext_cds_rna.pl), each combination (feature, gene and product) was counted using MySQL, and the CDS and rRNA feature table (11_Database_mtRecords_final.xlsx) was created. The table with feature combinations was manually examined to assign each gene. Sequences that could not be assigned unambiguously to one of the thirteen protein-coding genes or one of the two ribosomal RNA genes were discarded. Accession number, feature, position, gene, product and gene abbreviations of those assigned sequences were extracted using MySQL. The mitochondrial fasta file that was prepared as described above was partitioned in 15 individual fasta files with sequences of each mitochondrial-encoded gene.

Taxonomic information curation
We extracted the GI number from the mitochondrial fasta file using a custom perl script (15_ext_gi.pl). Then, taxonomy ID and taxonomy ranking were extracted using the gb_taxonomy_tool (https://github. com/spond/gb_taxonomy_tools). These two output files were combined into a taxonomy file using a custom perl script (19_rdp_train_3.pl). RDP classifier 15 utilizes only eight taxonomic rankings (Root, Kingdom, Phylum, Class, Order, Family, Genus and Species); therefore, we extracted only those ranks. At this stage, we performed the following quality controls: (I) removal of sequences that did not have species name in the species rank; (II) removal of sequences containing the following text in the taxonomy ranks: 'cf.', 'aff.', 'sp.', 'environment', 'undescribed', 'uncultured', 'complex', 'unclassified', 'nom.', 'nud.' and 'unidentif' (because these terms indicate uncertainty of species identity); (III) removal of sequences with the following identifiers 'sp0936BC', 'MG98.09', 'sp0942A' and 'EEG-2007' (since these are obviously not Latin names). This quality filtration was performed using a combination of Perl, MySQL and Unix commands.

Creation of taxonomy rank file for RDP classifier
First, we counted the number of occurrences of each taxonomy rank. Next, we extracted the eight taxonomic rankings. Then, the taxonomy rank file was formatted in two steps using two custom Perl scripts (trainset_db_taxid.pl and trainset_db_taxid_parent_2.pl). On some occasions, a conflict of taxonomic names was observed, such as the same genus name for animals in different higher taxonomic groups. Such cases, caused by duplicated taxonomic names above the species level (which are prohibited within the Metazoa but occur through error) caused the taxonomic assignment software to report error messages and abort analysis. In those cases, we made some modifications to the taxonomic name, such as addition of a distinguishing number to one of them (e.g., from Automolus to Automolus01).

Code availability
Scripts used for the present reference datasets preparation are freely available from the site www.reference-midori.info/download.php#.

Data Records
All reference datasets are freely available from the Midori reference web site (www.reference-midori.info) and also Dryad Digital Repository (Data Citation 1). Midori-UNIQUE and Midori-LONGEST (see usage notes for more information) are available in two formats, one compatible with the RDP Classifier 15 and the other compatible with SPINGO 16 . Both formats are compatible with BLAST+ 17 ; however, because the RDP format contains taxonomic rank information, we recommend using an RDP formatted file as an input to build the local BLAST+ 17 database (see usage notes for more information on how to build the local database). The numbers of sequences included in the datasets are listed in Table 1.

Technical Validation
Three kinds of quality filtrations were performed in the present study. First, we included only sequences that have binomial names (genus and species names). Here, we assume that species-level identifications of metazoan are more likely to be performed by well-trained taxonomic specialists, although this step does not ensure the absence of mislabelled sequences, it increases taxonomic accuracy, particularly at higher taxonomic levels. For example, a specimen identified at the species level is more likely to have a correct genus name. Second, we also performed systematic sequence length restrictions by removing extremely long or short gene sequences in the original nucleotide datasets (restriction limits are dependent on the genes, see Methods for details). We observed that sequences at the extremities of the length distribution of each gene were more likely the result of mis-annotations. We also observed that the RDP Classifier consistently provided erroneous sequence taxonomic assignments if very long or very short sequences were included in the reference datasets. Third, we attempted to detect and remove taxonomically mislabelled sequences in the datasets. To do so, we performed a high similarity (99%) clustering of the reference datasets (Table 3) within the typical range of intra-specific variations for mitochondrial-encoded genes. We assumed that, given the low threshold, sequences in each cluster should possess same taxonomic labels at least at the higher phylum-, class-and order-level. The clustering was performed on Midori-UNIQUE using UCLUST 18 and clusters were flagged if they contained sequences of multiple phyla, classes or orders. To identify which sequence was mislabeled in each flagged cluster, we performed a similarity search using the BLAST server 19 (blastn with 'low complexity region filter' and 'mask for lookup table only' function disabled). The distance tree functionality on the BLAST server was used to explore phylogenetic relationships with 100 close matches to each query. Sequences confidently identified as mislabelled were deleted from the datasets. Whenever we could not confidently determine which sequence of the cluster was mislabelled we retained all sequences. Overall, we found that the number of such cases was very low ( Table 3) which indicates that the upstream quality filtration was effective (e.g., the inclusion of sequences with binomial names only). We did not attempt to detect mislabelled sequences at lower taxonomic levels, such as family, genus and species in the present study as they may result from a range of evolutionary factors (e.g., incipient species, hybridization, introgression) but also from issues related to classification and systematics (e.g., synonyms, taxonomic confusions, taxonomic revisions).
In some cases, we observed missing taxonomic information, such as class, order or family name (e.g., accession numbers JN392469 and JF760210). In those cases, estimation of statistical support for the missing taxonomic level is not feasible. Therefore, we recommend including all standard levels of taxonomic names in the GenBank taxonomy.