An annotated draft genome for Radix auricularia (Gastropoda, Mollusca).

Introduction

Gastropods are one of the broadest distributed eukaryotic taxa, being present across ecosystems worldwide. They occupy a maximally diverse set of habitats ranging from the deep sea to the highest mountains and from deserts to the Arctic, and have evolved to a range of specific adaptions (Romero et al. 2015, 2016). However, as for molluscs in general, whose species richness is second only to the arthropods (Dunn and Ryan 2015), gastropods are highly underrepresented among publicly available genomes (fig. 1). To date, only eleven mollusc genome sequences—of which six are from gastropods—exist with varying qualities concerning contiguity and completeness (table 1). Any additional genome sequence has therefore the potential to substantially increase the knowledge about molluscs in particular and animal genomics in general.

F

The number of described species (Dunn and Ryan 2015; GIGA Community of Scientists 2014) and the fraction of sequenced genomes (http://www.ncbi.nlm.nih.gov/genome/browse/ on September 1, 2016). Animal phyla were obtained from (Dunn et al. 2014). Phyla with genomic record are displayed. Note the logarithmic scaling.

Table 1

Available Mollusc Genomes

Species	Assembly Length/Estimated Genome Size = % Assembled	#sequences/N50 (*contigs)	Coverage/Technology	Gap [%]	BUSCOs Present	Number of Annotated Proteins
Octopus bimaculoides a	2.4 Gb/2.7 Gb = 89%	151,674/475 kb	92/Illumina	15.1	73.8	23,994
Dreissena polymorpha b	906 kb/1.7 Gb c = 0.06%	* 1,057/855 bp	3/Roche 454	0	0	—
Corbicula fluminea d	663 kb/?	* 778/849 bp	3/Roche 454	0	0	—
Crassostrea gigas e	558 Mb/890 Mb f = 62.7%	7,659/402 kb	100/Illumina	11.8	82	45,406
Mytilus galloprovincialis g	1.6 Gb/1.9 Gb h = 86%	* 2,315,965/1067 bp	17/Illumina	0	1.6	—
Lottia gigantea i	360 Mb/421 Mb j = 85%	4,469/1870 kb	8.87/Sanger	16.9	97.0	23,822
Patella vulgata k	579 Mb/1,460 Mb = 39.7%	295,348/3160 bp	25.6/Illumina	0.00062	16.6	—
Conus tribblei l	2,160 Mb/2,757 Mb = 78%	1,126,156/2681 bp	28.5/Illumina	0	44	—
Aplysia californica m	927 Mb/1,760 Mb j resp. 1,956 Mb n = 53 resp. 47%	4,332/918 kb	66/Illumina	20.4	94.1	27,591
Biomphalaria glabrata o	916 Mb/929 Mb c = 99%	331,401/48 kb	27.5/Roche 454	1.9	89.1	36,675
Lymnaea stagnalis p	833 Mb/1,193 Mb q = 70%	* 328,378/5.8 kb		0	88	—

Note.—An overview from column 2 can be found in supplementary figure 4, Supplementary Material online. Column 5: Fraction of N’s in the assembly. Column 6: BUSCOs: (Benchmarking Universal Single-Copy Orthologs) NMetazoa =843; Present=complete+fragmented.

References: Genome sizes are from the genome publications, if not cited separately.

(Albertin et al. 2015).

(Peñarrubia, Sanz, et al. 2015).

(Gregory 2003).

(Peñarrubia, Araguas, et al. 2015).

(Zhang et al. 2012).

(González-Tizón et al. 2000).

(Nguyen et al. 2014).

(Rodríguez-Juíz et al. 1996).

(Simakov et al. 2013).

(Hinegardner 1974).

(Kenny et al. 2015).

(Barghi et al. 2016).

(Moroz et al.) GCF_000002075.1.

(Lasek & Dower 2013).

(Matty Knight, Coen M. Adema, Nithya Raghavan, Eric S. Loker) GCF_000457365.1.

(unpublished—Ashworth Laboratories 2016) GCA_900036025.1.

(Vinogradov 1998).

The pulmonate freshwater snail genus Radix has a holarctic distribution (Glöer, Meier-Brook 1998; Cordellier et al. 2012) and plays an important role in investigating climate change effects in freshwater ecosystems (Sommer et al. 2012). The number of European species as well as the precise evolutionary relationships within the genus is controversial. This is mainly due to weak morphological differentiation and enormous environmental plasticity across species (Pfenninger et al. 2006). Members of the genus are simultaneously hermaphroditic (Jarne and Delay 1990; Yu et al. 2016) and both outcrossing and self-fertilization occur (Jarne and Delay 1990; Jarne and Charlesworth 1993; Wiehn et al. 2002). The genus Radix is studied in many different fields, including parasitology (Huňová et al. 2012), evolutionary development (Tills et al. 2011), developmental plasticity (Rundle et al. 2011), ecotoxicology (Hallgren et al. 2012), climate change (Pfenninger et al. 2011), local adaptation (Quintela et al. 2014; Johansson et al. 2016), hybridization (Patel et al. 2015) and biodiversity (Albrecht et al. 2012). Despite this broad range of interests, genomic resources, are scarce and limited to transcriptomes (Feldmeyer et al. 2011, 2015; Tills et al. 2015) and mitochondrial genomes (Feldmeyer et al. 2010).

Here, we present the annotated draft genome sequence for Radix auricularia L. (fig. 2). This serves as an important foundation for future genomic and applied research in this scientifically important genus.

F

Photograph of Radix auricularia. Picture by Markus Pfenninger.

Materials and Methods

Sample Collection and Sequencing

Snails were collected from a pond in the Taunus, Germany, identified with COI barcoding (Pfenninger et al. 2006) and kept under laboratory conditions for at least five generations of inbreeding by full-sib mating. Three specimens of R. auricularia (fig. 2) were used for DNA extraction. Pooled DNA was used for preparation of three paired end and three mate pair (2, 5, and 10 kbp insert size) libraries, that were sequenced on an Illumina HiSeq 2000 and 2500 at Beijing Genomics Institute, Hong Kong (supplementary note 1 and table 1, Supplementary Material online). Reads were cleaned of adapter sequences using Trimmomatic 0.33 (Bolger et al. 2014; supplementary note 7, Supplementary Material online) and screened for contaminations with FastqScreen 0.5.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/; last accessed February 22, 2017; supplementary note 8 and fig. 7, Supplementary Material online). Raw reads have been deposited under NCBI BioProject PRJNA350764.

Genome Size Estimation

Genome size was estimated by flow cytometry based on a modified protocol of (Otto 1990; supplementary note 9, Supplementary Material online). Additionally, we estimated the genome size from our sequence data by dividing the total sum of nucleotides used for the assembly by the peak coverage from mapping back the assembly reads with the bwa mem algorithm from BWA 0.5.10 (Li 2013; supplementary note 4, Supplementary Material online). Re-mappings were also used to estimate the repeat content of the genome (supplementary note 5, Supplementary Material online).

Assembly Strategy

Reads were assembled using the Platanus 1.2.1 pipeline (Kajitani et al. 2014) with k-mer sizes ranging from 63 to 88 and a step size of 2. All other assembly parameters were kept at the default value. The output of the Platanus pipeline was filtered for sequences≥500 bp. Afterwards, scaffolding was performed using SSPACE 3.0 (Boetzer et al. 2011) with “contig extension” turned on. To further increase the contiguity of the draft genome, we applied a third scaffolding step, making use of the cDNA sequence data. Transcriptome contig sequences of R. auricularia and three closely related species (supplementary note 10, Supplementary Material online) were mapped sequentially according to phylogeny (Feldmeyer et al. 2015) using BLAT 35 (Kent 2002), with -extendThroughN enabled apart from default settings, onto the scaffolds; the gapped alignments were then used for joining of sequences with L_RNA_scaffolder (Xue et al. 2013). Finally, all sequences with at least 1,000 bp were used as input for GapFiller 1.10 (Boetzer et al. 2012) to close extant gaps in the draft genome. Details of the assembly can be found in supplementary note 11, Supplementary Material online.

Annotation Strategy

Metazoan core orthologous genes were searched in the R. auricularia assembly and all other available mollusc genomes using BUSCO 1.2b (Simão et al. 2015).

The whole annotation process was performed using the MAKER2 2.31.8 pipeline and affiliated programs (Cantarel et al. 2008; Holt and Yandell 2011). Initially, we built a custom repeat library from the assembly using RepeatModeler 1.0.4 (Simit and Hubley 2015) and read data using dnaPipeTE 1.2 (Goubert et al. 2015) with 30 upstream trials on varying coverage depths and then 50 parallel runs on the best-fitting coverage of 0.025 (supplementary note 12, Supplementary Material online). The draft genome and transcriptome of R. auricularia (supplementary note 10, Supplementary Material online) in addition to the BUSCO 1.2b (Simão et al. 2015) annotations of core metazoan genes on the draft genome were used as input for the initial training at the Augustus webserver (Stanke et al. 2004; http://bioinf.uni-greifswald.de/webaugustus/; last accessed February 22, 2017). As additional input for MAKER2, we created two hidden Markov models on the gene structure of R. auricularia. One was generated by GeneMark 4.32 (Lomsadze et al. 2005) and another by SNAP 2006-07-28 (Korf 2004), using the output of CEGMA v2.5 (Parra et al. 2007; summarized results in supplementary table 9, Supplementary Material online). We ran three consecutive iterations of MAKER2 with the draft genome sequence, the transcriptomes (supplementary note 10, Supplementary Material online), models from Augustus, SNAP and GeneMark, the repeat library and the Swiss-Prot database (accessed at May 23, 2016). Between the iterations, the Augustus 3.2.2 (Stanke et al. 2004) and SNAP models were retrained according to the best-practice MAKER2 workflow (supplementary note 13, Supplementary Material online). Finally, all protein sequences from MAKER2 output were assigned putative names by BLASTP searches (Camacho et al. 2009) against the Swiss-Prot database. In addition we used the targeted ortholog search tool, HaMStR v. 13.2.6 (Ebersberger et al. 2009; http://www.sourceforge.net/projects/hamstr/; last accessed February 22, 2017) to screen for 1,031 evolutionarily conserved genes that predate the split of animals and fungi. HaMStR was called with the options -strict, -checkcoorthologsref, and -hitlimit=5. The profile hidden Markov models that served as input for the search are included in the HaMStR distribution.

We created orthologous groups from protein sequences of all six annotated molluscs and 16 additional nonmollusc spiralian species with OrthoFinder 0.7.1 (Emms and Kelly 2015). All proteins were functionally annotated using InterProScan 5 (Zdobnov and Apweiler 2001; Quevillon et al. 2005). The enrichment analyses were performed in TopGO (Alexa and Rahnenfuhrer 2016), a bioconductor package for R (R Development Core Team 2008). We tested for significant enrichment of GO terms in proteins private to Radix and proteins found in all molluscs but Radix. We applied a Fischer’s exact test, FDR correction and filtered by q-values smaller than 0.05. Additional information can be found in supplementary note 6, Supplementary Material online.

Results and Discussion

Genome Assembly

A total of 1,000,372,010 raw reads (supplementary note 1and table 1, Supplementary Material online) were generated and assembled into 4,823 scaffolds (table 2; supplementary table 2, Supplementary Material online). The mitochondrial genome (13,744 bp) was fully reconstructed, evidenced by comparison to the previously published sequence (Feldmeyer et al. 2015). Re-mapping the preprocessed reads revealed that 97.6% could be unambiguously placed, resulting in a per position coverage distribution with its peak at 72× (supplementary note 2, Supplementary Material online; fig. 3). Additionally estimated insert sizes from mate pair libraries match their expected size (fig. 3).

Table 2

Characteristics of the Radix auricularia Genome and Draft Assembly

Parameter	Value
Haploid chromosome number	17 (Garbar and Korniushin 2003)
Estimated genome length	1.51 Gb (Vinogradov 1998)
Flow cytometry	1.58 Gb ± 21.5 Mb (this study)
Sequencing coverage	1.60 Gb
Total assembly length	0.91 Gb single copy or high complexity regions
#scaffolds	4,823
N50	578,730 bp
Gaps	6.4% N
Coverage	72x
Estimation of gene completeness	94.6% of BUSCO genes present
Gene prediction	17,338 genes
Gene space (UTR, Exons, Introns etc.)	200.6 Mb = 21.9% of assembly
Gene length (median)	8.0 kb
Gene fragmentation	147,195 exons
Exon space	25.3 Mb = 2.8% of assembly (1.6% of total genome)
Exon length (median)	125 bp
Protein length (median)	332 AA

Note.—BUSCOs: (Benchmarking Universal Single-Copy Orthologs) NMetazoa =843; Present=complete+fragmented.

F

Re-mapping statistics. For details, see supplementary note 2, Supplementary Material online. (A) Coverage distribution per position. The peak is located at a coverage of 72×. The x-axis is given in log-scale. (B) Insert size distributions for the three mate pair libraries with insert sizes of 2, 5, and 10 kb. The high fraction of mate pairs with insert sizes close to 0 is due to the repetitive nature of the Radix genome (supplementary table 3, Supplementary Material online). In particular, repeat stretches that are not properly resolved in the genome assembly interfere with a proper placement of reads.

The cumulative length of all scaffolds sums up to 910 Mb, which is about 665 Mb below the genome length estimates resulting from flow cytometric analyses (1,575 Mb; supplementary note 3, Supplementary Material online) and from a read-mapping analysis (1,603 Mb; supplementary note 4, Supplementary Material online). Both genome size estimates are consistent. This indicates an approximately uniform coverage of the nuclear genome in shotgun libraries without substantial bias introduced during library generation. This difference in length is most likely caused by a high repeat content in the Radix genome. Within scaffolds, 40.4% of the sequence content was annotated as repeats mostly at the ends of contigs (fig. 4). This, in combination with a pronounced increase of read coverage at contig ends (fig. 4) is typical for collapsed repeat stretches. The overall repeat content of the genome was estimated to be approximately 70% (supplementary note 5). The majority of repeats were either classified as Transposable Elements or as “unknown” (supplementary fig. 3, Supplementary Material online). The difference between genome size and assembly length of this R. auricularia draft assembly resembles that of other published mollusc genomes (supplementary fig. 4, Supplementary Material online). However, when considering contiguity reflected in the N50 value it ranks among the top mollusc genomes (tables 1 and 2). To evaluate completeness of the assembly’s gene space we used BUSCO (Simão et al. 2015) in combination with the provided metazoan set and recovered 94.6% of the subsumed genes. This suggests no conspicuous lack of gene information.

F

Collapsed repeats. (A) Coverage of continuous unambiguous sequence parts of the scaffolds. Outliers from boxplots are not shown. The red line represents the most frequent coverage of 72 × (fig. 3 A). (B) Positions annotated as repeats along continuous unambiguous sequence parts of the scaffolds.

Genome Annotation

The annotation resulted in 17,338 protein coding genes (table 2) of which 70.4% show a significant sequence similarity to entries in the Swiss-Prot database (e-value<10− 10, accessed on May 11, 2016). The number of identified genes is at the lower end compared with other annotated mollusc genomes (Min: Lottia gigantea 23,822; Max: Crassostrea gigas 45,406; table 1; supplementary table 4, Supplementary Material online). Thus, predicted Radix proteins were screened for completeness regarding evolutionary conserved genes using HaMStR (Ebersberger et al. 2009). The analysis resulted in a recovery of 93.7% and is in line with the results from BUSCO. Extrapolating completeness estimates of both tools suggests that the annotation covers the majority of genes being present in the draft genome sequence. We then checked how the differences in protein numbers could be explained. The fraction of orthogroups (cluster of orthologous genes; see “Materials and Methods” section) containing only one sequence per species was highest in Radix (supplementary fig. 5, Supplementary Material online). Moreover, there was a negative correlation (R 2 =0.77; P =0.02) between the number of annotated proteins per species and fraction of orthogroups containing only one sequence per species (supplementary fig. 6, Supplementary Material online). One explanation for this observation could be lineage specific gene duplications in the other mollusc lineages. Additionally, artificial gene fissions in the course of genome annotation may be less common in Radix. This might be attributed to our use of comprehensive transcriptomic data of Radix for guiding gene prediction.

Next to the evolutionarily old genes represented in the BUSCO and HaMStR gene sets, Radix contains 1,481 genes for which we could find no orthologs in the other mollusc and additional nonmollusc spiralian gene sets (supplementary table 5 and note 6, Supplementary Material online). We tested for over-representation of functional categories in genes private to Radix, as well as in genes present in all molluscs but Radix. We identified 17 Gene Ontology (GO) terms to be significantly enriched amongst the 1,481 proteins private to Radix compared with all other mollusc and nonmollusc spiralian protein sets available. Enriched terms include “nucleoside transmembrane transport”, “carbohydrate metabolic process,” and “chitin catabolic process” (supplementary table 6, Supplementary Material online). Among the categories found in all annotated molluscs but Radix (supplementary table 7 and note 6, Supplementary Material online), the “G-protein coupled receptor signaling pathway” is the most prominent one. G-protein receptors are involved in reactions to “hormones, neurotransmitters, and environmental stimulants” (Rosenbaum et al. 2009). The loss of these genes could have led to reduced sensitivity to such stimuli in Radix. Whether the reduced number of G-coupled receptor pathway components is biologically meaningful, or a result of technical and analytical limitations, cannot be determined from the present data. Membrane proteins, for example, are generally more diverse than water soluble proteins in the tree of life (Sojo et al. 2016), so we hypothesize that its proteins could be highly modified and were thus not identified as such in Radix.

Conclusion

Here we present a draft genome of the snail Radix auriculara. The genome is comparable in size to other mollusc genomes and also rich in repeats. This new genomic resource will allow conducting future studies on genome evolution, population genomics, and gene evolution within this genus and higher gastropod and mollusc taxa.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.