Bivalves

The completeness of the genome assembly was estimated with BUSCO v.3 [89 (link)], using a set of 843 conserved metazoan single-copy orthologs as a reference, and the resulting data about the present, fragmented, duplicated, and missing gene models were compared with previous genome assembly efforts carried out in M. galloprovincialis [18 (link), 19 (link)] (Additional file 1: Data Note 1.3.3). The residual presence of artefactual duplications was assessed with the Kmer Analysis Toolkit [90 (link)]. Consensus gene models were obtained by combining transcript alignments generated with PASA v 2.0.2 [91 (link)], bivalve protein alignments created with SPALN v2.2.2 [92 (link)], and ab initio gene predictions obtained with GeneID [93 (link)], GeneMark-ES [94 (link)], GlimmerHMM [94 (link)], and Augustus [95 (link)]. Evidences derived from these methods were assigned different weights and combined into consensus CDS predictions with EvidenceModeler-1.1.1. Gene models were subjected to an additional round of quality control to refine the annotation of UTRs and alternatively spliced exons (Additional file 1: Data Note 2.1 and 2.2). Gene models were functionally annotated with InterPro [96 (link)], KEGG [97 (link)], Blast2GO [98 (link)], SignalP [99 (link)], and NCBI CDsearch [100 (link)] (Additional file 1: Data Note 2.3). The gene models supported by PASA, but lacking a CDS, were considered as non-coding genes and included in a separate annotation track (Additional file 1: Data Note 2.5).
The completeness and integrity of gene models, as well as the genome assembly size and the number and density of gene models, were compared with several other recently sequenced molluscan genomes (Additional file 1: Data Note 3). Each gene model was assigned a support level (high, mild, or low) based on evidence obtained from Lola gills and digestive gland transcriptomes, as well as from several publicly available M. galloprovincialis RNA-seq datasets (Additional file 1: Data Note 4).

GenBank data was parsed using a combination of command-line and custom Perl scripts using BioPerl modules [22 (link)]. Tabular data was formatted using Python and plotted in R [23 ]. We use the terminology from Nilsson et al., (2005) and refer to taxa identified to the species rank as ‘fully identified’ and all other taxa as ‘insufficiently identified’ [24 (link)]. We also focused on NCBI nucleotide data deposited from 2003, the year COI barcoding was first introduced to the community, to present (2017) [25 (link)].
The names and taxonomic identifications for all Eukaryotes annotated to the species rank were retrieved from the NCBI taxonomy database using the Entrez query "Eukaryota[ORGN]+AND+species[RANK]" with an ebot script [Accessed November 3, 2017] [26 ]. Taxa were filtered according to the contents of the species field so that only fully identified taxa with a complete Latin binomial (genus and species) were retained. Entries that contained the abbreviations sp., nr., aff., or cf. were discarded. The remaining species names were formatted for use in the next query [species list]. For each year from 2003–2017 [year], records in the NCBI nucleotide database containing COI sequences were retrieved using the Entrez query "("CO1"[GENE] OR "COI"[GENE] OR "COX1"[GENE] OR "COXI"[GENE]) AND "Eukaryota"[ORGN] AND [year][PDAT]) AND [species list]” [2003–2016, accessed November 2017; 2017, accessed April 2018]. GenBank records were parsed, retaining information on year of record deposition and number of fully identified records. For fully identified records, sequence length as well as country and/or latitude-longitude fields were parsed.
We also assessed the number of high quality COI sequences that meet the standards developed between the INSDC and the Consortium for the Barcode of Life by looking for the BARCODE keyword in the GenBank record [11 (link)]. For each year from 2003–2017 [year], records in the NCBI nucleotide database containing COI BARCODE sequences were retrieved using the Entrez query "("CO1"[GENE] OR "COI"[GENE] OR "COX1"[GENE] OR "COXI"[GENE]) AND "Eukaryota"[ORGN] AND [year][PDAT] AND “BARCODE”[KYWD]) AND [species list]”. Fully identified and geotagged records were parsed as described above.
For our application example on freshwater biomonitoring, we retrieved a high-level list of relevant groups from Elbrecht and Leese (2017) to facilitate comparisons across studies [27 (link)]. Target freshwater taxa included: Annelida classes Clitellata and Polychaeta; Insecta (Arthropoda) orders Coleoptera, Diptera, Ephemeroptera, Megaloptera, Odonata, Plecoptera, and Trichoptera; Malacostraca (Arthropoda) orders Amphipoda and Isopoda; Mollusca classes Bivalvia and Gastropoda; and Platyhelminthes class Turbellaria. Within these groups there are likely to be non-freshwater taxa included, however, this method allowed us to quickly gauge the representation of freshwater taxa contained therein. These are also the same groupings often used to summarize results from COI freshwater biomonitoring assessments. A detailed look at specific freshwater taxa at finer taxonomic levels is beyond the scope of this paper and will be published elsewhere. For each freshwater target group we queried the NCBI taxonomy database for records identified to the species rank as described above. These taxon ids were concatenated and used to query the NCBI nucleotide database as described above. We assessed the representation of freshwater indicator taxa in the NCBI nucleotide database and level of annotation as described above.
For our application example on IUCN endangered animal species, we retrieved a list of endangered species names from http://www.iucnredlist.org from all available years (1996, 2000, 2002–2004, 2006–2017) filtering the results for native Animalia species [Accessed Dec. 12, 2017]. We excluded insufficiently identified species containing the terms ‘affinis’, ‘sp.’, or ‘sp. nov.’, leaving us with a list of 4,289 endangered animal species as well as 2,089 synonyms. We submitted this combined list of species names to the ‘NCBI Taxonomy name/id Status Report Page’ (https://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi) and retrieved a list of 2,613 taxon ids. For each taxon id, we queried the NCBI taxonomy and nucleotide databases as described above.
To assess the number of COI records unique to the BOLD database compared with the NCBI nucleotide database, we also retrieved records from the BOLD Application Programming Interface (API) as well as from the data releases. Since the BOLD database contains records from several DNA barcode markers such as ITS rDNA for fungi and COI mtDNA for animals, it was necessary to target just the COI records. COI sequences were retrieved from the BOLD API (http://www.boldsystems.org/index.php/API_Public/sequence?) using the terms ‘marker = COI-3P|COI-5P&taxon = ‘ for each Eukaryote phylum except for Arthropoda which was queried separately for each class, and Insecta which was queried separately for each order to enable the download of complete files [Accessed Apr. 26, 2018]. Lists of Eukaryote phyla, Arthropoda classes, and Insecta orders were retrieved from the BOLD taxonomy browser (http://www.boldsystems.org/index.php/TaxBrowser_Home). COI records were also retrieved from the BOLD data releases (http://www.boldsystems.org/index.php/datarelease). All available releases of animal COI records up to and including Release 6.50v1 were individually downloaded and parsed. Note that the records retrieved from the data releases may not be as current as those retrieved through the BOLD API.

Twenty-four well-assembled lophotrochozoan genomes were selected for phylogenetic analysis, include one annelid (Helobdella robusta) as outgroup, 21 bivalves (Archivesica marissinica, Argopecten concentricus, Argopecten irradians, Conchocele bisecta, Crassostrea gigas, Crassostrea virginica, Cyclina sinensis, Gigantidas platifrons, Lutraria rhynchaena, Mactra quadrangularis, Mercenaria mercenaria, Mizuhopecten yessoensis, Modiolus philippinarum, Mytilus coruscus, Pecten maximus, Pinctada fucata, Pinctada imbricata, Ruditapes philippinarum, Saccostrea glomerata, Scapharca broughtonii, Sinonovacula constricta), 5 gastropods (Aplysia californica, Chrysomallon squamiferum, Lottia gigantea, Haliotis rufescens, Pomacea canaliculata), and 2 cephalopods (Octopus bimaculoides and Octopus vulgaris) [22 (link), 26 (link), 52 (link), 113 (link)–132 ]. SonicParanoid v1.3.0 was used to define gene family clusters among different species [133 (link)]. The amino acid sequences of one-to-one single-copy orthologous genes were used to reconstruct their phylogenetic topology. The protein sequences were aligned using MAFFT v7.407 under default settings [134 (link)], and then were concatenated for phylogenetic analysis using a maximum-likelihood method implemented in IQ-TREE v 2.0.6 with the “-m MFP” parameter was applied to each protein partition [135 (link)]. To estimate divergence times, the rooted maximum-likelihood tree, along with a concatenated fourfold degenerate site sequence extracted from single-copy CDS (coding sequence), was used as the input of MCMCtree software implemented in PAML v4.8 [136 (link)]. For calibration, nine nodes were constrained by either fossil records obtained from website of TimeTree.

To evaluate the temporal dynamics of expanded gene families during the evolution of C. bisecta, the nucleotide substitution rates of bivalves were calculated by the branch distance divided by the estimated divergence time using MCMCtree. With default settings of MAFFT and “-automated1” option of trimAl v1.4 [139 (link)], all paralogs of the target gene family were aligned to determine the time required for gene duplication. The Nei-Gojobori pairwise codeml method was used to determine the dN values for all aligned pairs. Divergence times of gene pairs were estimated using the equation T = K/2r [140 (link)], where T is the insertion time, and r is the nucleotide substitution rate. The relationships between different gene pairs are determined following the DupGen_Finder (https://github.com/qiao-xin/DupGen_finder) pipeline, using Nematostella vectensis as a reference genome.

For expansion and contraction analysis, in addition to genome of C. bisecta, we selected 7 representatives with good BUSCO performance out of the 15 collected bivalve genomes. The selected genomes include that of the only two published endosymbiotic bivalves (A. marissinica, G. platifrons) [22 (link), 26 (link)], and 5 asymbiotic bivalves (Crassostrea gigas, L. rhynchaena, M. yessoensis, Modiolus philippinarum, P. fucata) which were separated in different bivalve clades and not known to host chemosynthetic bacteria. Before analysis, HMMSCAN (HMMER v3.1) was applied to identify Pfam domains in protein-coding gene sequences among the selected bivalve. The Pfam domains of the respective species were counted to construct a data frame, while multiple copies of a same domain in the same gene were counted as one.
Gene family analyses in the symbiotic bivalves were conducted using one-tailed Fisher’s exact tests for either expansion or contraction. In detail, for gene expansion/contraction at the protein domain level, we first calculated the counts of each Pfam domain in each genome of the 8 analyzed species, and the Pfam domain counts in each of the symbiotic bivalves (A. marissinica, C. bisecta, G. platifrons) was compared against the background average domain counts of the five asymbiotic bivalve genomes (Crassostrea gigas, L. rhynchaena, M. yessoensis, Modiolus philippinarum, P. fucata), which method was employed by Sun et al. for comparative genomic analysis [26 (link)]. Furthermore, we conducted the same analysis with the gene counts of each KEGG ortholog on each of the three symbiotic bivalves. After that, Pfam domain or KEGG ortholog with a P value less than 0.05 is considered statistically expanded or contracted in the three symbiotic bivalves. Finally, the evolutionary patterns of A. marissinica, C. bisecta, and G. platifrons were compared according to the expansion/contraction results by Fisher’s exact tests.
For phylogenetic analysis of each gene family, we employed Muscle v3.8.31 for multiple sequence alignment [137 (link)], and the phylogenetic trees were constructed with FastTreeMP v2.1.10 [138 (link)]. Specially, for reported expansion events of subfamilies of hemoglobin, which were not included in the KEGG database, we performed additional alignment with the sequences mentioned in Ip et al. [22 (link)] using Diamond, and phylogenetic analysis was conducted as the same.