A draft genome assembly of the solar-powered sea slug Elysia chlorotica

Elysia chlorotica, a sacoglossan sea slug found off the East Coast of the United States, is well-known for its ability to sequester chloroplasts from its algal prey and survive by photosynthesis for up to 12 months in the absence of food supply. Here we present a draft genome assembly of E. chlorotica that was generated using a hybrid assembly strategy with Illumina short reads and PacBio long reads. The genome assembly comprised 9,989 scaffolds, with a total length of 557 Mb and a scaffold N50 of 442 kb. BUSCO assessment indicated that 93.3% of the expected metazoan genes were completely present in the genome assembly. Annotation of the E. chlorotica genome assembly identified 176 Mb (32.6%) of repetitive sequences and a total of 24,980 protein-coding genes. We anticipate that the annotated draft genome assembly of the E. chlorotica sea slug will promote the investigation of sacoglossan genetics, evolution, and particularly, the genetic signatures accounting for the long-term functioning of algal chloroplasts in an animal.


Background & Summary
Many species of sacoglossan sea slugs are able to intracellularly sequester chloroplasts from their algal food, a phenomenon known as kleptoplasty, that is not observed in other clades of animals. In some sacoglossan species, the captured chloroplasts (usually called kleptoplasts) are maintained and capable of photosynthesis for one to several months, earning these molluscs the title of "solar-powered sea slugs" [1][2][3] . Among them, Elysia chlorotica, where the kleptoplasts are obtained from the filamentous alga Vaucheria litorea, is particularly interesting because it can retain functional chloroplasts in the cells of its digestive diverticula and survive without food supply for ten months to one year 2,3 . The mechanism that keeps 'stolen' chloroplasts functioning requires special proteins produced by nuclear genes of the algal host 4 . While there is a great deal of evidence using polymerase chain reaction (PCR) [5][6][7][8][9][10] , western blot 11,12 , RNA-seq 13 , and fluorescent in situ hybridization (FISH) investigations 14 that algal nuclear genes are present in the sea slug, genomic resources are scarce for E. chlorotica, limited to a mitochondrial genome assembly 7 , a few transcriptomes 13,15,16 and a low-coverage genome sequencing dataset of eggs 17 . There are no nuclear genome assemblies, even fragmented ones, publicly available for E. chlorotica so far. From an evolutionary perspective, although Mollusca represents the second largest animal phylum with around 85,000 extant species 18 , a fairly limited number of mollusc genomes have been sequenced yet [19][20][21][22][23][24][25][26][27][28][29][30][31] , with only 23 genomes publicly available on NCBI genome database (https://www.ncbi.nlm.nih.gov/genome/ browse/#!/overview/mollusca; access on December 5, 2018). Particularly, no reference genome has been generated for any sacoglossan mollusc.
In this study, we present the first draft genome assembly for the representative solar-powered sea slug E. chlorotica, which was assembled from Illumina short and PacBio long reads using a hybrid and hierarchical assembly strategy. We anticipate that this well-annotated draft genome assembly and the massive sequencing data generated in this study will serve as substantial resources for future studies of the evolution of sacoglossan molluscs, and particularly, for the investigation of the genetic basis underlying the long-term maintenance of algal chloroplasts in these sea slugs.

Sample collection, library construction and sequencing
Specimens of the sea slug E. chlorotica (NCBI taxonomy ID 188477; Fig. 1) were collected from a salt marsh near Menemsha on the island of Martha's Vineyard, Massachusetts in 2010. From there, the animals were shipped to Tampa, Florida, for maintenance in aquaria containing sterile, artificial seawater (1000 mosm; Instant Ocean, VA, USA) on a 14:10 light-dark cycle at 10°C as described in Pierce et al. (2012) 13 .
In Tampa, total DNA, including slug genomic, mitochondrial and algal chloroplast DNA, was extracted from a whole adult specimen that had been starved for at least 2 months using a Nucleon Phytopure DNA extraction kit (GE Healthcare UK limited, Buckinghamshire, UK) according to the manufacturer's instructions. The same kit was used to extract total DNA, including slug genomic and mitochondrial DNA, from a batch of~1000 larvae that had not hatched from the egg capsules and never been fed. The larvae do not have any chloroplasts. A total of 11 Illumina DNA paired-end (PE) libraries were constructed according to the standard protocol provided by Illumina (San Diego, CA, USA), including six libraries with short-insert sizes (170 bp × 2, 500 bp × 2, and 800 bp × 2) from the adult DNA, and five mate-paired libraries with long-insert sizes (2.5 kb × 2, 5 kb × 2, 10 kb × 1) from the larval DNA. Sequencing was performed for all the 11 libraries on the HiSeq 2000 platform according to the manufacturer's instructions (Illumina, San Diego, CA, USA), using the modes of PE100 for all the short-insert libraries and PE49 + PE90 for each of the five mate-paired libraries. A total of 296.73 Gb of Illumina reads were produced (Data Citation 1 and Data Citation 2), which can cover the estimated haploid genome size of E. chlorotica by k-mer analysis for 516 times (Table 1).
In addition, 9.45 Gb (16X) PacBio long-reads with a mean subread length of 1.2 kb and N50 subread length of 1.7 kb were sequenced for another DNA sample (Data Citation 1 and Data Citation 2), which was extracted from another starved adult specimen, that had been shipped frozen to Okazaki, using the CTAB method 32 and purified with DNeasy Plant mini kit (Qiagen). Three libraries were constructed according to the 6 kb library construction protocol and sequenced in 93 SMRT cells on the PacBio RS platform using the C2 chemistry following the manufacturer's instructions (Pacific Biosciences, Menlo Park, CA, USA).

Estimation of genome size and heterozygosity
Prior to downstream analyses, all the Illumina reads were submitted to strict quality control using SOAPnuke (v1.5.3) 33 . Duplicated reads arising from PCR amplification during library construction, adapter-contaminated reads and low-quality reads were removed using parameters -l 7 -q 0.4 -n 0.02 -d -t 10,0,10,0 for the short-insert (i.e. 170 bp, 500 bp and 800 bp) data and -l 7 -q 0.35 -n 0.05 -d -S for the long-insert (i.e. 2.5 kb, 5 kb and 10 kb) data, yielding a total of 176.32 Gb of clean Illumina reads ( Table 1).
All of the clean reads from the six short-insert libraries, except those derived from algal chloroplasts, slug mitochondria and the previously reported endogenous retrovirus of E. chlorotica 34  considered to be derived from organelles or retrovirus if a pair of reads was mapped to any of the three genomes (Data Citation 3-5) by BWA-MEM (v0.7.16) 35 , and such read pairs were discarded, resulting in a total of 62.8 Gb Illumina data for k-mer analysis. The haploid genome size of E. chlorotica was estimated to be around 575 Mb according to k-mer frequency distributions generated by Jellyfish (v2.2.6) 36 using a series of k values (17, 19, 21, 23, 25, 27, 29 and 31) with the -C setting, which was calculated as the number of effective k-mers (i.e. total k-merserroneous k-mers) divided by the homozygous peak depth ( Table 2). But this estimated haploid genome size might be an underestimate, as certain parts of the E. chlorotica genome (e.g. GC-extreme regions) may have failed to be sequenced due to technical limitations 37 , and/or repetitive sequences may not have been resolved properly by k-mer analysis given that mollusc genomes are generally known to be repeat rich [19][20][21][22][23][24][25][26][27][28][29][30][31] .

Genome assembly
The E. chlorotica genome was assembled by a hybrid and hierarchical assembly strategy as described below: (i) Clean reads from the Illumina short-and long-insert libraries were assembled into contigs using ALLPATHS-LG (v52488) 40 with default parameters except setting HAPLOIDIFY = True, which yielded an initial assembly with a total length of 776 Mb and a contig N50 of 1.7 kb. This initial assembly was~35% longer than the estimated genome size of 575 Mb, indicating that, for some genomic regions, two haploids were assembled separately due to high heterozygosity. Thus, (ii) we used HaploMerger2 (v20151124) 41 to separate the two haploid sub-assemblies from the initial ALLPATHS-LG assembly, and an assembly with a total length of 575 Mb and contig N50 of 1.9 kb was produced. Next, (iii) we assembled a separate genome with the PacBio long-reads alone. Sequencing errors in the PacBio reads were first corrected by the clean Illumina reads from 170 bp and 500 bp short-insert libraries using PacBioToCA (v8.3) 42 with parameter -length 400, and 5.36 Gb (9.32 X) error-corrected PacBio reads were retained (Table 1). Then the error-corrected PacBio reads were assembled using Canu (v1.4) 43 with parameters minReadLength = 400 minOverlapLength = 400 contigFilter = 2 400 1.0 1.0 2, which produced an assembly with a total length of 469 Mb and contig N50 of 4.4 kb. (iv) The PacBio assembly was merged with the above HaploMerger2 assembly with Metassembler (v1.5) 44 , which resulted in an improved assembly with a total length of 535 Mb and contig N50 of 5.1 kb. (v) These resulting contigs were further assembled into scaffolds using the distance information provided by read pairs from the Illumina short-and long-insert libraries with SSPACE (STANDARD-3.0) 45 . Specifically, prior to scaffolding, the read pairs were aligned to the contigs using BWA (v0.6.2), and the insert size of each library was inferred from the statistics of a pre-run of SSPACE based on satisfied pairs in distance and orientation within contigs. Then scaffolding was performed with SSPACE using the estimated insert size of each library with the minimum allowed insert size error setting to be 0.3 for short-insert libraries and 0.5 for long-insert libraries. Subsequently, (vi) intra-scaffold gaps were filled using PBJelly from PBSuite (v15.8.24) 46,47 with the error-corrected PacBio long reads by setting minReads = 3, followed by using GapCloser (v1.10.1) 48 with the Illumina short-insert paired-end reads by setting library insert sizes according to SSPACE estimation as described above. (vii) The gap-filled scaffolds were submitted to HaploMerger2 again to reduce redundant sequences, followed by polishing with all the Illumina shortinsert clean reads by PILON (v1.22) 49 . Finally, (viii) potential contaminants in the assembly including sequences from algal chloroplasts, slug mitochondria and adaptor/vector as identified by the NCBI contamination-screening pipeline were removed by an in-house script. The improvements of assembly generated at each step of the assembly process were presented in Table 3.
The final result was a genome assembly with a total length of 557 Mb, comprising 9,989 scaffolds (Data Citation 6,7). The contig and scaffold N50s of this assembly were 28.5 kb and 442.0 kb, respectively, and unclosed gap regions represented 3% of the assembly (Table 4), exhibiting a continuity comparable to other published molluscan genomes (Data Citation 8-21). In addition, GC content of the E. chlorotica assembly excluding gaps was estimated to be 37.7%.

Protein-coding gene annotation
We applied a combination of homology-based, transcriptome-based and de novo prediction methods to build consensus gene models for the E. chlorotica genome assembly. For homology-based prediction, protein sequences of Aplysia californica, Caenorhabditis elegans, Crassostrea gigas, Drosophila melanogaster, Lottia gigantea and Homo sapiens were first aligned to the E. chlorotica assembly using TBLASTN (blast-2.2.26) 53 with parameters -F F -e 1e-5. Then the genomic sequences of the candidate loci together with 2 kb flanking sequences were extracted and submitted to GeneWise (wise-2.4.1) 54 for exonintron structure determination by aligning the homologous proteins to these extracted genomic sequences with settings of -sum -genesf -gff -tfor/-trev (-tfor for genes on forward strand and -trev for reverse strand). For transcriptome-based prediction, we collected published RNA-seq data from Pierce et al.    Finally, gene models from the above three methods were combined into a non-redundant gene set using a similar strategy as Xiong et al. (2016) 59 . Briefly, the homology-based gene models were first integrated with the transcriptome-based models to form a core gene set, followed by integration with the de novo models. Then de novo models not supported by homology-based and transcriptome-based evidence were also added to the core gene set if BLASTP (blast-2.2.26; parameters -F F -e 1e-5) 53 hits could be found in the UniProtKB/Swiss-Prot database (v2018_05). Finally, genes showing BLASTP (blast-2.2.26; parameters -F F -e 1e-5) hits to transposon proteins in the UniProtKB/Swiss-Prot database (v2018_05) were removed from the combined gene set. We ultimately obtained 24,980 protein-coding genes with up to 22,717 (90.9%) supported by RNA-seq signal ( ≥ 5 RNA reads).
To assign gene names for each predicted protein-coding locus, we first mapped the protein sequences of all the 24,980 genes to the UniProtKB/Swiss-Prot database (v2018_05) using BLASTP (blast-2.2.26) 53 with parameters -F F -e 1e-5. Then the best hit of each gene was retained based on its BLASTP bit score, and the gene name of this best hit was assigned to the query E. chlorotica gene. A similar process was performed against the NCBI nr database (v20180315). In addition, we performed functional annotation with InterProScan (v5.29-68.0) 60 to examine motifs, domains, and other signatures by searching against the databases including ProDom 61 , PRINTS 62 , Pfam 63 , SMART 64 , PANTHER 65 and PROSITE 66 . To determine what pathways the E. chlorotica genes might be involved in, protein sequences of the E. chlorotica genes were searched against the KEGG database (v87) 67 with BLASTP (blast-2.2.26) using parameters -F F -e 1e-5. As a result, 21,452 (85.9%) of the predicted protein-coding genes were successfully annotated by at least one of the four methods (Table 6).

Code availability
The bioinformatic tools used in this work, including versions, settings and parameters, have been described in the Methods section. Default parameters were applied if no parameters were mentioned for a tool.

Data Records
Raw reads from Illumina and PacBio sequencing are deposited in the NCBI Sequence Read Archive (SRA) database with accession number SRP156455 and Bioproject accession PRJNA484060 (Data

Technical Validation
To evaluate the quality of the E. chlorotica assembly, we first aligned the 62.8 Gb Illumina short-insert clean reads which were previously used for k-mer analysis to the assembly using BWA-MEM (v0.7.16) 35 , and observed that 96% could be mapped back to the assembled genome with 85% of the mapped reads being aligned in proper pairs as counted by samtools flagstat (SAMtools v1.7) 68 Fig. 3). Moreover, single-position coverage analysis by samtools bedcov based on the above BWA-MEM alignment with PCR duplicates removed by Picard (v2.10.10) 70 revealed that 98% of the assembly excluding gaps were covered by ≥ 5 reads and resulted in a per-position coverage distribution with its peak at 98X (Fig. 4). Of note, a minor peak at 49X was also observed (Fig. 4), implying that either a small amount of redundant sequences was still present in the final assembly, or the genomic regions corresponding to the minor peak  were highly polymorphic between two haploids so that reads from the unassembled allele could not be aligned to the assembled allele by BWA-MEM due to low sequence identity. We also used Picard CollectInsertSizeMetrics (v2.10.10; setting MINIMUM_PCT = 0.5) to analyze the insert size distribution of each Illumina library based on BWA-MEM (v0.7.16) alignment of read pairs from each library, and observed that the estimated insert sizes of all the libraries matched their expected fragment sizes (Fig. 5). These results indicated that most sequences of the E. chlorotica genome that were captured by the sequencing platforms are present in the current assembly with proper orientation.  The quality of the E. chlorotica assembly was further assessed by REAPR (v1.0.18) 71 , a tool that evaluates the accuracy of a genome assembly using mapped paired-end reads. Specifically, all the shortinsert clean reads and the 10 kb mate-pair reads were aligned to the final assembly by reapr smaltmap for calling error-free bases and scaffolding errors, respectively. It is noteworthy that REAPR recommends using the longest insert data with sufficient fragment coverage for calling scaffolding errors while data from multiple long-insert libraries are available 71 . In our case, the fragment coverage of the 10 kb matepair library calculated by REAPR peaked at 752X (Fig. 6), far beyond the minimum requirement of 15X. Ultimately, REAPR judged 80.45% of the bases in the E. chlorotica assembly as error free (i.e. bases covered by ≥ 5 perfectly and uniquely mapped reads), and identified 123 collapsed repeats and a total of 2,943 fragment coverage distribution (FCD) errors. An FCD error usually represents incorrect scaffolding, a large insertion or deletion in the assembly 71 . Considering the high heterozygosity (3.66%), the high repeat content (32.6%) and especially the high tandem repeat content (10.3%) of this genome, which likely affect the performance of read mapping, we believe that the accuracy of the E. chlorotica assembly is acceptable. As a comparison, more than 7,000 FCD errors are recently reported in an improved genome assembly of the Atlantic cod Gadus morhua, of which the assembly size (643 Mb) and the tandem repeat content (11% of assembly) are actually comparable to the E. chlorotica assembly 72 .
Next, we employed Benchmarking Universal Single-Copy Orthologs (BUSCO, v3.0.2) 73 , a software package that can quantitatively measure genome assembly completeness based on evolutionarily informed expectations of gene content, to evaluate the completeness of the E. chlorotica assembly and 14 other published molluscan genomes using 978 genes that are expected to be present in all metazoans. We found that 927 (94.7%) of the expected genes were present in the E. chlorotica assembly with 913 (93.3%) and 14 (1.5%) identified as complete and fragmented, respectively. Only 51 (5.3%) genes were considered missing in the E. chlorotica assembly. The completeness of the E. chlorotica assembly based on BUSCO assessment was overall comparable to other published molluscan genome assemblies (Table 4).
Finally, we evaluated the completeness of the annotated gene set of E. chlorotica with BUSCO (v3.0.2) and DOGMA (v3.0) 74 , a program that measures the completeness of a given transcriptome or proteome based on a core set of conserved domain arrangements (CDAs). BUSCO analysis based on the metazoan dataset showed that 968 (98.9%) of the expected genes were present in the E. chlorotica gene set with 948 (96.9%) identified as complete. A higher number of expected genes were identified by BUSCO in the annotated gene set than in the E. chlorotica genome assembly, probably because searching genes in a transcriptome or proteome is simpler than in a genome. Meanwhile, DOGMA analysis based on PfamScan Annotations (PfamScan v1.5; Pfam v32.0) 63 and the eukaryotic core set identified 93.3% of the expected CDAs in the annotated gene set. These results demonstrate the completeness of the annotated gene set of the E. chlorotica assembly.