Background & Summary

Sea lice are marine copepods that negatively impact the salmon aquaculture worldwide. Two of the most studied sea lice species are Caligus rogercresseyi and Lepeophtheirus salmonis1,2. Annually the salmon farming industry accounts $480 million in losses associated with sea lice, representing 10% of production costs3,4,5. The parasitism on farmed fish causes skin damage, immunosuppression, and co-infection of opportunistic pathogenic bacteria6,7,8. Like all ectoparasites, lice spend a large part of their life cycle on a fish host, displaying specific mechanisms for evading the host’s immune response9,10,11.

The life cycle of lice species is complex and consists of several instars divided by moults. For instance, C. rogercresseyi comprises two larval stages (nauplius I, nauplius II and copepodite), four juvenile stages (chalimus I - IV) and one adult stage (female or male)12. During the copepodite stage, the process of host identification occurs, preparing the lice for infestation and settlement8. The successful infestation process on the host allows the parasite access to nutrients for reproduction and adult development13,14. Previous studies have shown that lice have developed physical mechanisms of host recognition. Among these, lice can identify the temperature of the water, salinity changes, and detect the swimming of fish15. Host identification via detection of semiochemicals has also been reported16. In C. rogercresseyi, the presence of advanced chemoreceptors that are capable of identifying specific molecules of different host species has recently been described17,18,19. Herein, the gene family of ionotropic receptors (IRs) are pivotal molecular components for the salmon-louse interaction20,21.

Molecular understanding of C. rogercresseyi is pivotal to develop sustainable salmon aquaculture. However, genomic resources in this species are limited and poorly characterized at functional levels. In 2012, Yasuike et al. (2012) reported a compilation of genomic information on different sea lice genera, including C. rogercresseyi. It was not until 2014 that Gallardo-Escárate et al.22 reported the transcriptome of different life stages during the ontogenetic development as well as differences between male and female adults. This transcriptomic resource served as a basis to identify genes involved in molting, cuticle formation, myogenesis, metabolism, immune response, nervous system development and reproduction. Notably, this gene set has served as a basis for the design of new vaccines2.

The increasing availability of transcriptome data has revealed the importance of non-coding RNAs as key regulators of the mRNA transcription23. To date, microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) have been studied in several arthropod species with special emphasis on parasitic vectors24. Long non-coding RNAs are sequences greater than 200 nucleotides, transcribed in a similar way as coding RNAs25. It has been suggested that the number of lncRNAs has increased during evolution, where organisms with more complex mechanisms have acquired more lncRNAs to control diversifying biological processes. MicroRNAs are transcripts around 22 base pairs in length that play an important role in post-transcriptional gene regulation26. There are studies that show that miRNAs are not only regulators of biological processes, but can also participate in parasite-host interaction processes9,27. In insects affected by viruses, it has been observed that viruses are capable of releasing miRNAs that can regulate the expression of their host genes in order to successfully establish the infection28. For C. rogercresseyi, several miRNAs expressed during the different stages of development have been characterized29,30. Within the profile of miRNAs characterized in C. rogercresseyi, the miRNA annotated as Bantam is highly expressed in the infective stage of copepodid. This suggests that Bantam has a key role in the success of the infection. Taken together, these resources reported for the sea louse C. rogercresseyi represent a valuable tool to develop sustainable control strategies in the salmon industry. What is lacking is an annotated genome that will facilitate an integrated examination of molecular interactions and provide insight in evolutionary and epigenetics processes that underlie critical life history characteristics. In this study, we report the chromosome-scale whole genome sequence of C. rogercresseyi through application of Pacific Biosciences’ single molecule sequencing technique (SMRT) and Phase Genomics’ proximity ligation (Hi-C) based genome scaffolding.

Methods

Sample collection, NGS libraries, and sequencing

Adult female specimens of Caligus rogercresseyi were collected from Atlantic salmon (Salmo salar) at the Caligus Reference Laboratory (CRL), University of Concepción, Chile (Fig. 1). With the aims to reduce the heterozygosity or the number of individuals per pool, female lice were selected for whole-genome sequencing. The samples were frozen in liquid nitrogen to preserve DNA quality, and ten females were used for genomic DNA isolation. High quality DNA was isolated using the Qiagen DNA purification kit (QIAGEN, Germantown, MD, USA) following the manufacturer’s instructions. It is important to note that sea lice are marine copepods exposed to marine environmental conditions and consequently to commensal microorganisms. To reduce bacterial DNA contamination, lice were treated with 20 mg/ml ampicillin (Sigma-Aldrich, USA), 20 mg/ml Kanamycin (US biological, USA), 1x Penicillin-Streptomycin (GIBCO, USA), 100 ug/ml Primocin (Invivogen, USA) for 72 hr prior to the DNA extraction protocol31. Furthermore, lice from different developmental stages were separately collected, fixed in RNA Later solution (Ambion, USA), and stored at −80 °C until RNA extractions.

Fig. 1
figure 1

The sea louse Caligus rogercresseyi. Adult female (right) and adult male (left). Magnification 10x.

Full size image

Genomic DNA libraries were constructed according to the manufacturer’s protocols for genome assembly (Table 1). SMRT sequencing yielded 38.32 Gb long reads from 8 SMRT cells (Table 1S). The subreads N50 and average lengths were 11,093 and 6,824 bp, respectively. Hi-C libraries were constructed from whole animals using Phase Genomics’ Animal Hi-C kit and sequenced on an Illumina’s Hiseq4000 platform to yield 238 million of reads. Short-read sequencing libraries were prepared using an insert size of 150 bp obtained from 1 μg of genomic DNA, after fragmentation, end-paired, and ligation to adaptors, respectively. The ligated fragments were fractionated on agarose gels and purified by PCR amplification to produce sequencing libraries.

Table 1 Sequencing data generated for sea louse C. rogercresseyi genome assembly and annotation.
Full size table

For transcriptome sequencing, RNA libraries were constructed from nauplius I, nauplius II, copepodid, Chalimus I-II, Chalimus III-IV, males and females, and sequenced by Illumina technology according to the manufacturer’s protocols (Table 1). Briefly, total RNA was extracted from 10 parasites from each stage using the Trizol reagent method (Invitrogen, USA). The quality and integrity of extracted RNAs was measured in a TapeStation 2200 instrument (Agilent, USA), using the R6K Reagent Kit based on manufacturer’s instructions. RNA samples >9 in RIN numbers were selected for library preparation. For whole transcriptome sequencing, 2 μg of total RNA was used for dscDNA libraries with TruSeq Total RNA kit (Illumina, USA). RNA libraries quantification was conducted by qPCR using the NEBNext Library Quant Kit for Illumina (New England Biolabs, USA). The sequencing was performed using the MiSeq platform (Illumina, USA) using a 2 × 250 bp paired-end reads scheme (single flow cell per developmental stage). In addition to generating conventional RNA-seq for 6 developmental stages, small-RNA libraries were also constructed using TruSeq Small RNA Kit (Illumina, USA) for each stage, with libraries run in 41 single-end cycles. Small-RNA libraries were simultaneously sequenced using barcodes according to the manufacturer’s protocols. In total 3 flow cells were used to sequence the 6 developmental stages. A total of 52.01 and 28.18 Gb was yielded for transcriptome and miRNome characterization, respectively (Table 1).

De novo assembly of C. rogercresseyi genome

With eight single-molecular real-time cells in the PacBio Sequel platform, we generated 38.32 Gb of high-quality DNA genome information. These long subreads were assembled with the Canu V1.5 package32 using default parameters, yielding a draft genome for the sea louse equivalent to 727 Mb with contig N50 of 43,366 bp and 35.55 GC%. The draft genome was assembled with CANU in 25,608 contigs (Table 2). The size genome assembly made by CANU was comparable with previous genome size reported for closely related species33,34. However, the manual curation of a subset of contigs revealed bacterial DNA contamination. As we previously mentioned, antibiotic treatment was applied to reduce the natural lice microbiota. However, it appeared that some fraction of the bacterial burden still remained despite the antimicrobial compound used. To reduce the bacterial DNA contamination, all contigs assembled by Canu were firstly filtered against NCBI prokaryotic reference sequence database and then against the reference C. rogercresseyi transcriptome (Table 3). For the first filter, BLASTx was applied with an expectation value of 10.0, word size = 3, filter low complexity, protein matrix and gap costs = BLOSUM62, Existence, 11-1, meanwhile that for the second filter a mapping approach was implemented with the following settings using CLC Genomics Workbench V12 (Qiagen, USA): match score = 1, mismatch cost = 2, cost of insertions and deletions = Linear gap cost, insertion cost = 3, deletion cost = 3, length fraction = 0.5, similarity fraction = 0.8, global alignment = No, non-specific match handling = map randomly. Taking advantage of the two filters, we removed all the contigs with a significant match to bacterial DNA, reducing the number of contigs produced by the draft genome for C. rogercresseyi made by Canu from 25,608 to 17,711 contigs. Here, the new dataset yielded a draft genome assembly of 519.19 Mb with an N50 of 38,179 bp (Table 4). Notably, the DNA contamination produced by the natural microbiota found in C. rogercresseyi was ∼30%. This fact shows the importance of the microbiota in louse biology, revealing putative associations with the pathogenesis of this ectoparasite.

Table 2 Genome assembling using PacBio SMRT sequencing in C. rogercresseyi.
Full size table
Table 3 Statistics of transcriptome de novo assembly for the sea louse C. rogercresseyi.
Full size table
Table 4 Statistics of genome assembly and Hi-C analysis for the sea louse C. rogercresseyi.
Full size table

Chromosome assembly of C. rogercresseyi using chromatin interaction mapping analysis

In vivo Hi-C is a technique that maps physical DNA-DNA proximity across the entire genome35,36. The method was introduced as a genome-wide version of its predecessor, 3 C (Chromosome Conformation Capture)37, and has been used as a powerful tool in chromosome-scale genome assembly of many animals in recent years. In this study, Hi-C experiments and data analysis on adult females were used for the chromosome assembly of the sea louse C. rogercresseyi. Here, two Hi-C libraries were prepared and sequenced by Phase Genomics (Seattle, WA, USA), resulting in ∼100x coverage and ∼238 million 150-bp paired-end reads (Table 4). The Hi-C analysis evidenced that 46.70% of high-quality reads analysed showed intercontig signals or Cis-close position (<10kbp on the same contig), and an additional 5.32% of sequence reads revealed a Cis-far conformation (>10Kbp on the same contig). To order and orient the 17,711 contigs Hi-C reads were aligned using Bowtie238 and scaffolding performed using Proximo (Phase Genomics, Seattle, WA, USA). We then applied Juicebox39 for visual inspection and manual correction. We also manually removed 7,897 scaffolds that were microbe-sized and disconnected from the rest of the assembly. We obtained the first chromosome-level high-quality C. rogercresseyi assembly with an N50 scaffold of 29.78 Mb, providing a useful genomic resource for research in sea louse biology and also, to develop novel control strategies applied to the salmon aquaculture (Table 5). In order to visualize the scaffold’s length construction, the in vivo Hi-C data were used to generate 21 pseudo-chromosomes assembled with PacBio consensus long DNA reads (Fig. 2). The largest scaffold was assembled from 1,235 contigs, a size of 36.77 Mb. Meanwhile, the smallest scaffold was 7.98 Mb of length and consisted of 396 original contigs (Fig. 3). Notably, the number of contigs in scaffolds were 16,931 (100% of all contigs in chromosome clusters, 95.6% of all contigs) and 505.27 Mb of genome size (100% of all length in chromosome cluster, 97.33% of all sequence length). The completeness of genome assembly was assessed by the single-copy ortholog set (BUSCO, V3.0.2) against Eukaryota, Metazoa, and Arthropoda40. The results indicated a complete BUSCO of 78.9% [S:75.3%, D:3.6%] and a fragmented BUSCO of 13.5% [M:13.6%, n:303].

Table 5 De novo assembly of C. rogercresseyi genome using chromatin interaction mapping.
Full size table
Fig. 2
figure 2

The sea louse Caligus rogercresseyi genome contig contact matrix using Hi-C data. The blue squares represent the draft scaffold. The color bar illuminates the Hi-C contact density in the plot.

Full size image
Fig. 3
figure 3

The sea louse Caligus rogercresseyi genome. The circos plot shows the genomic features for the 21 pseudo-chromosomes. A) GC content, B) Gene density, C) lncRNA density, D) miRNA density and E) Repetitive elements. The chromosome size is shown in Mb scale.

Full size image

Repetitive element and non-coding gene annotation in the C. rogercresseyi genome

Repetitive elements and non-coding genes in the sea louse genome were annotated by homologous comparison and ab initio prediction. RepeatMasker41 was used for homologous comparison by searching against the Repbase database and RepeatModeler42. According to these analyses, about 269.83 Mb Mb of repeat sequences were annotated, which accounted for 51.9% of the whole genome. Herein, DNA transposons, LINE, and LTR transposable elements were identified (Table 6). Useful genome information for population genetic studies is the identification of simple sequence repeats (SSRs) or microsatellites. The mining of SSRs revealed that the C. rogercresseyi genome has 441,494 SSR sequences, where 65.76% represent dinucleotide motifs (Table 7). The total of SSR sequences accounted for 1.39% of the whole genome, and the genome distribution was correlated with the chromosome size (Fig. 1S). Furthermore, SSRs type dinucleotides, and specifically, the motifs AC/GT were the most abundant, representing the 65.76% of the total microsatellite sequences (Fig. 2S). Trinucleotides and tetranucleotides were found in 32.46% of the SSRs sequences (Fig. 3 and Table 7).

Table 6 Classification and distribution of repeats based on RepeatModeler from C. rogercresseyi genome.
Full size table
Table 7 Simple Sequence Repeats (SSR) of C. rogercresseyi genome using SSR Finder analysis.
Full size table

Protein-coding genes prediction and functional annotation in the C. rogercresseyi genome

For the identification of protein-coding genes, two approaches were employed for the sea louse genome, including homologous comparison and ab initio prediction. For homologous comparison, the protein sequences from Caenorhabditis elegans (GCA_000002985.3), Drosophila melanogaster (GCA_000001215.4), and Daphnia pulex (GL732539.1) genomes were extracted using the respectively published genomes. and aligned against the sea louse genome using TBLASTN (e-value < 1e-5). Gene sequence structure of each candidate genes was predicted using GeneWise43. For ad initio prediction, five tools were used to predict protein-coding genes using the Genome Sequence Annotation Server “GenSAS” (https://www.gensas.org)44. Specifically, Augustus, Braker, GeneMarkES, SNAP, and GlimmerM were used with default parameters. Finally, a non-redundant reference gene set was generated using EvidenceModeler (EVM) and PASA2 tools45. Taken together 25,510 protein-coding genes were identified. (Fig. 3 and Table 8). Additionally, 437 tRNAs were predicted using tRNAscan-SE, and 39 rRNA genes were annotated using RNAmmer via GenSAS. For non-coding RNAs with putative regulatory roles, 5,774 miRNAs and 6,308 long-ncRNAs were identified and annotated within the C. rogercresseyi genome using transcriptome sequencing data (Fig. 4 and Table 9). For functional annotation, the predicted proteins within the sea louse genome were searched by homology against four databases of InterPro46, GO47, KEGG KO48, and Swissprot49. Overall, 88.05%, 68.85%, 64.02%, and 91.02% of genes matched entries in these databases, respectively. A total of 23,686 genes (93%) were successfully annotated by gene function and conserved protein motifs (Table 10).

Table 8 Prediction of protein-coding genes in the sea louse Caligus rogercresseyi genome.
Full size table
Fig. 4
figure 4

Stage-specific transcriptome analysis in the sea louse Caligus rogercresseyi. (A) Transcriptome patterns of coding genes during the lifecycle. The heatmap was based on Transcripts Per Million (TPM) calculation and hierarchical clustering on Manhattan distances with average linkage. White colors mean upregulated coding genes, blue colors downregulated genes. (B) Venn diagram showing shared and unique genes expressed among the six developmental stages. (C) GO enrichment of stage-specific genes (P-value ≤ 10–16;|fold-change| > 5) annotated for key biological processes differentially expressed. The radar plot represents the comparison between two developmental stages according the C. rogercresseyi lifecycle.

Full size image
Table 9 Summary of non-coding RNA annotation in the sea louse Caligus rogercresseyi.
Full size table
Table 10 Statistics for genome annotation of the sea louse Caligus rogercresseyi.
Full size table

Technical Validation

RNA integrity

Before constructing RNA-seq libraries, the concentration and quality of total RNA were evaluated using Agilent 2100 Bioanalyser (Agilent, USA). Three metrics, including total amount, RNA integrity, and rRNA ratio, were used to estimate the content, quality, and degradation level of RNA samples. In this study, only total RNAs with a total amount of ≥10 μg, RNA integrity number ≥8, and rRNA ratio ≥1.5 were finally subjected to construct the sequencing library.

Quality filtering of Illumina sequencing raw reads

The initial raw sequencing reads were evaluated in terms of the average quality score at each position, GC content distribution, quality distribution, base composition, and other metrics. Furthermore, the sequencing reads with low quality were also filtered out before the genome assembly and annotation of gene structure.

Table 11 Software and URLs.
Full size table

Data Records

DNA and RNA sequencing runs were deposited to NCBI Sequence Read Archive (SRA)50,51,52. The assembled genome has been deposited at NCBI assembly with the accession number ASM1338718v153. Additional files containing repeated sequences, gene structure, and functional prediction were deposited in the Figshare database54.