Isopoda

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.

We used amino acid translations to align sequences to the BIOCODE reference dataset using MACSE v1.00
[47 (link)]. Quality-filtered sequences were sequentially aligned and added to the reference dataset using the option “enrichAlignment”. This alignment strategy is only reasonable because the studied COI fragment is highly conserved at the amino acid levels. To further optimize computing time, sequences were split into subsets containing 500 sequences that were aligned in parallel thanks to a computer farm and then progressively merged into a single final alignment using the option “alignTwoProfiles”. MACSE can detect and quantify interruptions in open reading frames due to: (1) nucleotide substitutions that result in stop codons and (2) insertion or deletion of nucleotides (non multiples of three) that induce frameshifts. Sequences with stop codons are likely bacterial sequences, pseudogenes or chimeric sequences. On the other hand, frameshifts may be caused by sequencing errors that are frequent with the 454 platform
[48 (link)]. MACSE can also detect and quantify insertions and deletions that do not lead to interruptions in open reading frames. COI is relatively conserved and indels are relatively uncommon. For example, only 0.9% of the sequences in BIOCODE dataset (including platyhelminthes, gastropods and isopods) display a deletion of one codon in their COI sequence, and none of the sequences in the BIOCODE dataset have codon insertions. As a result, we decided to keep all sequences from the 454 dataset which satisfied the following criteria: no stop codons, no frame shifts, no insertions and less than four deletions. For the final dataset we retained all sequences with a single frameshift when they had no stop codon, no insertions and no deletions to account for sequencing errors. Alignment of these sequences with frameshift required insertion or deletion of a nucleotide either at the first, second or third codon position. However, because the correct position could not be known, we chose to remove these codons all together.

We sequenced 28 transcriptomes from 20 invertebrate taxa that lack genomic resources, but are well-suited for answering questions about the function, development, and evolution of eyes and other light-interacting structures (Table 1). For example, we generated transcriptomes from RNA expressed by the eyes and skin of certain cephalopod mollusks (squid and octopus). These animals may have the most complex light-influenced behaviors of any invertebrate [51 ,52 (link)], but it appears that the eyes of cephalopods tend to contain only a single spectral class of photoreceptor ([53 (link)]; though see [54 (link)] as an exception). Additional physiological complexity may be suggested by the results of high throughput sequencing. It is also possible that certain visually-influenced behaviors in cephalopods – such as dynamic camouflage – may be influenced by molecular components that are expressed outside of their eyes. For example, past work suggests that certain cephalopods express LIT genes in their light-producing photophores [55 (link)] and in certain dermal cells [56 (link)].
We also sequenced transcriptomes for a range of arthropods. We chose to study stomatopods (mantis shrimp) because they have an unsurpassed ability to distinguish different aspects of light. Certain species are maximally sensitive to twelve distinct wavelength peaks and some species can identify both linearly and circularly polarized light [57 (link)-60 (link)]. Similarly, we chose to study odonates (damselflies and dragonflies) because they have physiologically complex eyes [61 (link)] and display a diversity of visually-influenced behaviors [62 -64 (link)]. To study the degeneration of eyes in arthropods from subterranean environments, we examined certain species of isopods and crayfish in which closely related species or populations live either above or below ground. Specifically, we sequenced tissues from the eye-bearing, surface-dwelling isopod Caecidotea forbesi and its eyeless, cave-dwelling congeneric C. bricrenata. We also sequenced transcriptomes for different populations of the isopod Asellus aquaticus, which has a surface-dwelling form and multiple cave-dwelling populations with typical cave morphologies like degenerated eyes [65 (link),66 (link)]. Likewise, we generated transcriptome data from a pair of surface (Procambarus alleni) and cave (P. franzi) freshwater crayfish. Crayfish have previously been the focus of molecular evolutionary studies of opsin in cave/surface comparisons [67 (link)]. To study the evolution of sexually dimorphic eyes, we generated a transcriptome for the RNA expressed by developing eyes from the ostracodEuphilomedes carcharodonta, a species in which males have compound eyes, but females do not [68 ,69 (link)]. Other species in this family of ostracods exhibit a similar, but independently evolved eye dimorphism, suggesting that these ostracods may be a promising system for the study of sex-specific convergent phenotypic evolution [70 (link)].
Lastly, we sequenced transcriptomes for Tripedalia cystophora, a cubozoan cnidarian (box jellyfish). Cubozoans are the only cnidarians with camera-type eyes and, for that reason, have been the subject of numerous studies of visual neurobiology [71 (link)-74 (link)], morphology [75 (link),76 (link)], and behavior [77 (link),78 (link)]. Transcriptomic resources will aid these efforts. Further, as cnidarians, cubozoans may help us understand the evolutionary origins of the metazoan phototransduction cascade [79 (link)-81 (link)].

The genome of the Hepatincola symbiont of A. vulgare was sequenced from a female isopod collected from a natural population in Availles-Thouarsais, France (46° 51’ 37” N, 0° 8’ 28” E) in 2014. Hepatincola was known to be present in this population from our previous metabarcoding study [33 (link)]. DNA was extracted from both pairs of midgut glands using phenol-chloroform extraction. In total, 4.5 µg of DNA were used for size selection using AMPure XP beads (Beckman Coulter, USA) at a bead:sample ratio of 0.7x to enrich in long fragments. In total, 3.5 µg of DNA were recovered and used for library preparation using the Oxford Nanopore Ligation Sequencing kit SQK-LSK 108 (Oxford Nanopore Technologies, UK). The library was sequenced on an R9.4 flowcell on the MinION sequencer for 58 h. The run was stopped and restarted several times to optimize pore use. Basecalling was done using Albacore v2.0.1 using a quality threshold of Q7. After discarding low quality (A. vulgare (Accession: GCA_004104545.1) using Minimap2 v2.15 [38 (link)], resulting in 1,083,710 reads. The reads ≥1 kb were assembled using Canu v1.7 [39 (link)], producing a 1.37 Mbp circular contig containing two 16S rRNA genes 99% identical to the 16S rRNA gene sequence of Hepatincola from P. scaber (Accession: AY188585). This initial assembly was first polished with Nanopore reads using Nanopolish v0.11.1 [40 (link)] and subsequently with Illumina reads using Racon v1.4.3 [41 (link)]. Additional Illumina polishing with Polypolish v0.4.3 [42 (link)] did not detect additional errors. Illumina reads mapping onto the Hepatincola genome were extracted from A. vulgare shotgun metagenomic datasets from our previous study [43 (link)], in which DNA from different tissues (including the midgut glands) of the same Hepatincola-infected individual had been included.

Twenty terrestrial isopod species were screened for the presence of Hepatincola. Most tested individuals came from population cages maintained by the UMR CNRS 7267 at the University of Poitiers (France), except for Philoscia muscorum and Porcellionides pruinosus, for which field-collected individuals from Ensoulesse (France) were also included (Supplementary Table S1). For four species, several populations could be tested, resulting in a total of 25 terrestrial isopod populations. One male and one female of each population were tested for the presence of Hepatincola. To this end, the midgut glands were dissected in Ringer solution under a stereomicroscope. As there are two pairs of midgut glands per individual, one pair was used for DNA extraction and diagnostic PCR, while the second pair was fixed for TEM. DNA was extracted using phenol-chloroform extraction and PCR was performed by amplification of the 16S rRNA gene using the Hepatincola-specific forward primer 137F (5’-ACACGTGGGAATTTGGCT-3’) in combination with the “universal” reverse primer 520R (5’- ATT-ACC-GCG-GCT-GCT-GG-3’) [37 (link)].

The genome of the Hepatincola symbiont of P. dilatatus petiti was assembled from Illumina shotgun metagenomic datasets from our previous study [44 (link)]. Specifically, we used two datasets corresponding to the metagenomes of pooled midgut glands from seven P. dilatatus petiti males and females, respectively. After removal of potential host reads by mapping against a custom database containing all isopod sequences available on NCBI and unpublished sequences produced by the laboratory UMR CNRS 7267, the remaining reads were assembled using Megahit v1.0.3 [45 (link)] with custom parameters (--min-count 2 --k-min 21 --k-max 127 --k-steps 1). Contigs greater than 1500 bp were grouped into bins using MetaBAT2 v2.12.1 [46 (link)]. Bin quality was checked with CheckM v1.0.13 [47 (link)] and only bins with a completeness >80% and a contamination rate <5% were retained. RNAmmer v1.2 [48 (link)] was used to predict ribosomal RNAs for each bin, which led to the detection of two partial 16S rRNA genes with >99% identity to the 16S rRNA gene sequence of Hepatincola from P. scaber (Accession: AY188585) in a bin from the midgut glands of female P. dilatatus petiti. The raw reads were mapped onto the contigs of the bin using BOWTIE2 v2.2.9 [49 (link)] with the --very-sensitive option. Mapped reads were de novo assembled using SPAdes v3.13.0 [50 (link)] with the --careful option. This produced 66 contigs with a combined length of 1.27 Mbp. The contigs were ordered using Mauve v2.4.0 [51 (link)] using the complete Hepatincola genome from A. vulgare as a reference. This resulted in a genome scaffold of 1,229,614 bp.