Adult parasite samples of C. rogercresseyi corresponded to two previously characterized strains with contrasting sensitivities to the antiparasitic drug azamethiphos (Núñez-Acuña et al., 2020). The sensitivity to this drug was tested through a specific bioassay validated for this species, representing drug resistance/susceptibility as the EC₅₀ index, effective drug concentration affecting 50% of the population (SEARCH C, 2006; Marín et al., 2018), and the efficacy in field treatments in salmon farms. The selected resistant strain had an EC₅₀ index of 13.72 and azamethiphos treatment efficacy in salmon farms < 60%, while the susceptible population had EC₅₀ of 0.4778 and in-field treatment efficacy > 92%. Bioassays for each population consisted of the exposure of ten parasites (five females and five males) to a gradient of different concentrations of azamethiphos, including 1, 3, 8, 10, 30, 100, 200, and 300 ppb (µg/L water) of a drug. In addition, control unexposed groups of ten parasites were also included. Exposures were applied in triplicates (three groups of ten parasites per drug concentration) in a petri dish for 30 min. This exposure time was selected because it corresponds to the therapeutic time used on salmon farms. All the concentrations were used to calculate the EC₅₀ indexes, and the parasites exposed to 8 ppb and un-exposed controls were used to collect samples fixed in RNA Later solution (Ambion^®, Thermo Fisher Scientific™, Waltham, MA, USA) for molecular analyses.

The same groups of sea lice samples were used to perform DNA and RNA sequencing, including the control (unexposed samples to azamethiphos) and parasites exposed to 8 ppb of azamethiphos. Three pools of males and females from the same conditions were used for each type of sequencing, combining five parasites of each sex. DNA was extracted using the DNeasy^® Blood & Tissue Kit (QIAGEN, Hilden, Germany) following the manufacturer’s protocol. Sample homogenization was performed using ceramic beads in a Mixer Mill MM 200 homogenizer (RETSCH, Haan, Germany) at a frequency of 24 Hz for 5 min. DNA concentrations were measured using a Qubit 4 fluorometer (Invitrogen, Thermo Fisher Scientific™, Waltham, MA, USA) and the Qubit 1X dsDNA BR Assay kit. The quality of the extracted DNAs was observed in the TapeStation 2200 instrument (Agilent Technologies Inc., Santa Clara, CA, USA) by calculating the DIN (DNA Integrity Number). Only samples that were more significant than 9 in DIN were used for library preparations. The initial amount of DNA samples for library preparation was 1 μg, and then the Genomic DNA by Ligation kit (SQK-LSK109) was used for Nanopore library construction (Oxford Nanopore Technologies, ONT, Oxford, United Kingdom) following the manufacturer’s protocol and washing the samples applying the AMPure XP beads (Beckman Coulters, Brea, CA, USA). High-molecular weight on the DNA samples was expected. Thus 500 ng of each library was sequenced in the R.9.4. FLO-MIN106D flow cell on the MinION sequencing platform (ONT). Three sequencing replicates were conducted on each pool in separate runs.

Whole-transcriptome data corresponded to an RNA-sequencing study that was previously published using the same samples of this present study (Núñez-Acuña et al., 2020), SRA accession numbers: SRX9398676, SRX9398677, SRX9398678, SRX9398679, SRX9398646, SRX9398647, SRX9398644, SRX9398645, SRX9398656, SRX9398667, SRX9398674, SRX9398675). Briefly, total RNA from pools of parasites exposed and unexposed to azamethiphos drug was extracted by the Ribopure™ kit (Ambion^®, Thermo Fisher Scientific™, Waltham, MA, USA). Next, RNA integrity was calculated in a TapeStation 2200 instrument (Agilent Technologies Inc., Santa Clara, CA, USA), calculating the RNA Integrity Number (RIN) using the R6K Screentape kit. Only samples with RIN > 8.0 were used for dscDNA library construction, which was prepared from 2 μg of high-quality RNA of each sample using the TruSeq Total RNA kit (Illumina^®, San Diego, CA, USA). Then, RNA-sequencing was carried out in a MiSeq NGS sequencer (Illumina^®, San Diego, CA, USA) in a 2 x 250 paired-end format on three replicates for each experimental condition.

Adapter sequences from long DNA reads obtained by Nanopore sequencing were trimmed out using the Porechop software with default parameters and including a search in both forward/reverse directions and adapters in the middle of reads (https://github.com/rrwick/Porechop). Trimmed reads were imported to the CLC Genomic Workbench software (Qiagen Bioinformatics, Hilden, Germany) for mapping to the C. rogercresseyi genome previously published for this species (Gallardo-Escárate et al., 2021) (GenBank genome assembly accession number: ASM1338718v1). In parallel, RNA-seq reads obtained from the Illumina sequencing platform were imported into the same software and trimmed by quality, keeping only reads with Q30 values > 30. The Illumina reads had greater quality (fewer sequencing errors) than Nanopore reads; thus, polishing of DNA MinION reads was achieved by the MiSeq reads using a combination of the racon program (Vaser et al., 2017) and minimap long-read assembler, as was previously described (Li, 2016). The minimum sequence length was 1, the POA windows size was 500 in the racon polishing algorithm.

The DNA polished reads were mapped to the reference genome using the “long read mapping” tool in CLC using the following parameters: match score = 2; mismatch cost = 4; gap open cost = 4; gap extend cost = 2; long gap open cost = 24; long gap extend cost = 1; score bonus for global alignment = 0. Mapped reads on the genome were used to find genes in each mapping by using the “transcript discovery tool” in CLC with the following parameters: minimum length of ORF = 100; gene merging distance = 500; minimum reads in gene = 10; minimum predicted length = 250, exon merging distance = 100. Gene lists for the susceptible and resistant strains were saved in separate files.

Duplicated genes, including collinear genes (duplicated in different chromosomes – or different regions from the same chromosome) and tandem-repeated genes, were identified through an evaluation of the presence of syntenic blocks within the C. rogercresseyi genome. First, the Blastall software was used to compare the homology of all the predicted genes in the resistant and the susceptible strains, then MCScanX software (Wang et al., 2012) was used by each strain to establish significant alignments on genes considering an E-value < 1E-5 from the BLAST results. Significant collinear blocks were considered if the duplicated region on the genome involved at least five predicted genes and the E-value of the whole region was < 1E-5. Alignments were associated with chromosome positions in MCScanX by exporting a simplified GFF file from CLC Genomics, including the chromosome number, gene name (correlative number as ID), start and end positions of the genes for each strain. Tandem genes were identified if at least two consecutive genes had an E-value of < 1E-20 in the respective chromosome for each strain.

A set of collinear and tandem genes was exported by each sea louse strain, and the sequence was retrieved to identify gene functions and GO terms. Blast-X analyses were conducted on the collinear and tandem genes using a simplified version of the UniProt database (http://uniprot.org), considering all the protein sequences for arthropod species and, in parallel, compared to the SwissProt revised protein databases for all eukaryote species. Sequences with lower E-values in the two databases were used for further analyses. Sequences with significant BLAST hits were used for Gene Ontology annotation and Enrichment of GO terms on biological processes, molecular functions, and cellular components in the TBtools package (Chen et al., 2020).

The genes with significant BLAST hits by each strain were used as a reference to evaluate the expression levels by mapping back the RNA-seq Illumina reads corresponding to each strain. The “RNA-seq analysis” tool in CLC Genomics was used to calculate the TPM values of the duplicated genes of untreated and parasites exposed to azamethiphos. Expression differences between these conditions and between resistant and susceptible groups were evaluated using GLM-linear models in CLC Genomics, considering as significant those genes with variations in |Fold Change| > 4 and P-value < 0.01. Volcano plots and heatmaps based on a hierarchical cluster of TPM values were constructed using TBtools to visualize gene expression differences. Heatmaps were clustered based on Manhattan’s distance and with an average linkage. P-values were corrected through false discovery rate (FDR) calculation following the procedure of Benjamini and Hochberg (1995), applying multiple testing criteria with a minor modification named “automatic independent filtering”, which was previously reported (Love et al., 2014).

Due to the importance of discovering novel potential mechanisms for pharmacological resistance in the sea louse C. rogercresseyi, genes related to cuticle formation and penetration resistance (Balabanidou et al., 2018) were evaluated in the genomes obtained by DNA resequencing. Cuticle protein genes and similar genes associated with cuticle formation were searched in the BLAST results from the duplicated genes in resistant and susceptible strains. Three genes related to this function were found in the resistant strain and two in the susceptible group. The MEME Suite package (Bailey et al., 2015) was used to observe the positions of the coding sequence, UTR regions, and a specific motif of the complete cuticle protein genes. Syntenic alignments involving these genes were involved were highlighted in the SynVisio Synteny Browser software (Bandi, 2020). Expression analyses of cuticle genes were conducted using the same methodology applied in the previous point for duplicated genes.

Copy Number Variants (CNVs) were identified in the genome mappings obtained by Nanopore reads through the “depth-of-coverage” method, applying the CONTRA algorithm (Li et al., 2012) incorporated in the CLC Genomics Workbench. In this analysis, mapped reads on the genome from the resistant strains were defined as “case” and mapped reads from the susceptible strain as “control” to compare the two lice populations. As the target genes for CNVs prediction, we used the overlapping genes presented in both strains with an average coverage of > 20, which corresponded to 10,330 genes. The “graining level” parameter for the prediction of CNVs was “intermediate”, and the statistical cutoff for significant CNV calls was P-value < 0.05, and |fold change| > 2.0. P-values in this section were also corrected through FDR correction using the method mentioned in Section 2.4. This means that if the fold change was > 2.0, it represented an amplification in the resistant strain, and if the fold change was < -2.0, it represented a deletion in this population.

Expression analyses of genes containing significant CNVs were conducted by applying the same method for duplicated genes described in Section 2.4. In addition, an association between the position of significant CNVs with the gene density in the sea louse reference genome, mapping coverage by each strain, and collinear genes obtained in the resistant strain was conducted by constructing a CIRCOS plot in the TBTools package (Chen et al., 2020).

A resequencing strategy was performed, including DNA and RNA sequencing from two previously described C. rogercresseyi strains with contrasting drug sensitivities to azamethiphos (Núñez-Acuña et al., 2020; Núñez-Acuña et al., 2021b). One strain was resistant to this drug, and the other was susceptible. Genomic DNA reads obtained through Nanopore sequencing yielded 581,919 high-quality trimmed reads with an average length of 2,327 bp, representing 1.352 Gbases of DNA reads for the resistant strain. Transcriptome RNA sequencing obtained from the same strain using Illumina sequencing yielded 99,389,810 trimmed reads with an average length of 150.07 bp. Illumina sequencing has higher quality than third-generation sequencing technologies, such as Nanopore, so it was used to polish long reads obtained from the MinION sequencer to improve long DNA reads. The obtained polished window was 22.52% of the original Nanopore reads for the resistant strain. In this polished window, the quality of nanopore sequences was increased from an average PHRED score 17.22 to 39.65 after correction. After the polishing, 470,394 DNA reads (80.83% of the total number of reads) from the resistant strain were successfully mapped to the C. rogercresseyi draft genome as a reference. Nevertheless, the percentage of base pairs mapped to the reference was 93.68% because the mapped reads were longer than the 111,525 unmapped reads (2,693 versus 767 bp, respectively). The number of base pairs covered by this genome mapping was 1.267 Gbases, which represented 2.5X of the genome size, considering that the reference length was 0.506 Gbases. Regarding the full-gene sequences, the average coverage of the genes from this strain was 6.94X.

Regarding the susceptible strain, the genomic DNA obtained from Nanopore sequencing yielded 840,131 reads with an average length of 2,408 bp, representing 2.023 Gbases. Transcriptome sequencing of the susceptible strain using Illumina sequencing yielded 100,992,178 reads with an average length of 149.73 bp. Reads polishing resulted in a window of 21.71% with respect to the original number of base pairs from Nanopore reads, increasing the quality from 18.1 to 36.1 in the PHRED score. Therefore, 82.88% of polished reads were correctly mapped to the C. rogercresseyi genome, corresponding to 92.88% of the total length. The average length of the mapped 696,294 reads was 2,699 bp, and 1,002 bp for the 143,837 not-mapped reads. Considering this number, 1.879 Gbases of Nanopore DNA sequencing were mapped, representing 3.7 X of the reference genome. Furthermore, the average coverage of the predicted genes in this strain was 6.51X.

Putative genes were predicted through the “Transcript Discovery” tool of the CLC Genomics workbench from the mappings obtained from Nanopore reads to the reference genome. The resultant number of genes for the resistant strain was 18,438, with an average length of 15,677 bp, whereas for the susceptible strain, the number of genes was 18,525 with an average length of 15,575 bp. In addition, differential expression analyses were conducted from RNA sequencing analyses for each strain, comparing transcriptomes from parasites exposed to the azamethiphos drugs versus non-exposed lice (control), resulting in 631 differentially expressed genes (|fold change| > 4 and P-value < 0.01) in the resistant strain and 754 in the susceptible strain ( Figure 1 , first track).

Gene duplications were identified in the lice genome by sequence complementarity of predicted genes, by evaluating the microsynteny arrangements throughout all the chromosomes, and then categorized as collinear or tandem genes. From the 18,438 predicted genes in the resistant strain, 602 were collinear corresponding to 3.26%, and were grouped in 42 alignments between two chromosomes consisting of at least five genes per local arrangement (see red links in Figure 1 ). Additionally, tandem repeated genes, defined as genes subsequently duplicated in the same chromosome, corresponded to 274 genes or 1.49% of the total number of genes in this strain. Regarding the susceptible strain, the collinear genes were much lower, corresponding to 84 genes out of 18,525 genes, being 0.45% and grouped in six different alignments of at least five genes (see green links in Figure 1 ). Tandem duplicated genes were 520, corresponding to the 2.81% from the 18,525 total genes.

Gene categories according to their origin were similar in both strains, except for some genes with origins by Whole-Genome Duplications (WGD) or Segmental events, which were higher in the resistant strain (21.16%) than in the susceptible strain (18.3%). These categories also implied a slight reduction in the % of the Dispersed genes in the resistant strain compared to the susceptible lice ( Figure 2A ). The collinear genes are among the Whole-Genome Duplications (WGD) or Segmental events. Collinear and tandem duplicated genes were annotated by Gene Ontology ( Figure 2B ). Enrichment analyses of GO terms in all the duplicated genes on both strains showed that the nucleus is a cellular compartment enriched for collinear and tandem duplicated genes. Regarding molecular functions, the most enriched GO terms were similar for both types of duplicated genes, most of which were related to DNA integration, metabolic process and recombination, and transposition. Regarding this last term, the transposable elements (TEs) were highly abundant among the duplicated genes after Blast results ( Supplementary File 1 ). Molecular functions were also similar for both tandem and collinear duplicated genes, being different terms related to catalytic activity on DNA, nucleases, and polymerase activity.

Transcriptome expression patterns were analyzed in the R and S lice strain, focusing on collinear and tandem-duplicated genes ( Figure 3 ). The 686 collinear genes (including those collinear genes in R and those in the S strains) and the 794 tandem-duplicated genes had expression patterns with differential expression trends in R/S strains and samples exposed to azamethiphos ( Figures 3A, D ). Hierarchical clustering of TPM (transcripts per million of reads) values for the four experimental groups was observed to be associated with the drug exposure and not according to the resistant/susceptible strain. However, the clustering of features observed in the heatmaps showed various groups of genes that were differentially expressed when comparing the resistant and susceptible lice strains. Statistical comparisons were conducted in different pairs of samples to validate these differences ( Figures 3B, C, E, F ). Among all the collinear genes, 75 were differentially expressed (|fold change| > 4; P-value < 0.01) in the resistant population when comparing the exposed and control samples; the other 75 genes were differentially expressed when conducting this comparison in the susceptible population. However, these collinear genes were not the same, with only 19 differentially expressed in exposed vs. control samples in both populations ( Figure 3B ). Other statistical comparisons consisted of the combination of the TPM values for exposed and control samples and then comparing the fold change between resistant and susceptible populations (R vs. S in Figures 3B, E ); and combining TPM values for both populations and then comparing then fold change between exposed and control samples (E vs. C in Figures 3B, E ). The R vs. S comparison resulted in only 15 differentially expressed collinear genes, whereas the E vs. C comparison resulted in 29 differentially expressed collinear genes. Finally, 40 collinear genes were exclusively differentially expressed in the exposed vs. control comparison in the lice from the resistant strain, and 51 collinear genes were exclusively expressed in the susceptible strain ( Figure 3B ). Concerning tandem-duplicated genes, 93 were differentially expressed in the exposed vs. control comparison for the resistant strain and 107 were differentially expressed in the susceptible strain ( Figure 3E ). From these genes, 26 were differentially expressed according to this comparison in both strains. The R vs. S comparison for tandem genes resulted in 21 differentially expressed genes and the E vs. C comparison for 42 genes. Exclusive differentially expressed genes in the resistant population were 51, whereas 68 were for the susceptible lice. There was only one collinear gene and one tandem gene with differential expression values for all possible comparisons, which was an annotated gene for collinear genes (not found in public databases) and a retrotransposable element ORF2 for tandem genes.

One group of genes of particular interest was the cuticle protein genes, and another with functions on the cuticle formation was evaluated to identify duplication in lice genomes ( Figure 4 ). Five collinear duplicated cuticle genes were found in both strains (white links in Figure 4A, D ). Three cuticle-formation genes were collinear in the sea lice-resistant strain: cuticle protein 4, present in chromosome 4 and with duplication in chromosome 9; larval cuticle protein 65, in chromosome 5 and duplicated in chromosome 4; and cuticlin 4, in chromosome 10 with duplication in chromosome 5 ( Figure 4A ). In susceptible lice population, there were two collinear cuticle genes: cuticle protein 9 and cuticle protein 19, both present in chromosome 11 and duplicated in chromosome 4 ( Figure 4D ). The full sequence length of these genes ranged from 8,030 to 32,756 bp, with a coding sequence ranging from 378 to 1,566 bp. Additionally, conserved Motif 1 was repeatedly found throughout all the cuticle protein genes ( Figure 4B ). The expression patterns of these five genes differed depending on the strain. In the resistant strain, cuticle proteins 7, 8, and the larval cuticle protein 65 were upregulated in sea lice exposed to the delousing drug azamethiphos. In contrast, only cuticle protein 8 showed a similar pattern in the susceptible lice after exposure to the same drug ( Figure 4C ). The other two genes were downregulated in the susceptible strain. The cuticlin-4 gene was downregulated in both strains after exposure to azamethiphos, while the cuticle protein 19 showed no significant change in the resistant strain, but was upregulated in the susceptible lice ( Figure 4C ).

Coverage normalization analyses resulted in 10,330 genes with standard coverage for Copy Number Variants (CNVs) identification among both resistant and susceptible strains ( Figures 5A, B ). These 10,330 common genes in both strains were used for coverage comparisons using the susceptible strain as a reference and obtaining 134 genes with differential coverages indicating the presence of CNVs. Establishing a cutoff of |adjusted Fold Change| > 2 and P-value < 0.01, a total of 106 genes containing CNVs with duplication in the resistant population (Gain) and 28 genes containing CNVs with duplication in the susceptible lice (Loss) were found ( Figure 5C ). Most of the CNVs were in chromosome 6 (16 CNVs), 4 and 7 (12 CNVs in each one), which was not directly associated with the gene density on those chromosomes or the presence of collinear genes ( Figure 5D ). Additionally, there were 17 region CNVs, which were those variants involving wider genome regions, including more than a single-gene portion of the genome, and were distributed throughout all the chromosomes in the genome (“Gain” and “Loss” annotations in Figure 5D ). The highest fold change value for a region CNV was 6.39 (P-value = 2.23 E – 308) and corresponded to a region including the gene Rotatin, located in chromosome 6, amplified in the resistant strain (Gain). This gene was located on the first chromosome region, in-between an area of the genome that is also collinear with other genomic regions from chromosomes 3 and 4 ( Figure 5D ).

Genes located in areas with CNVs in the genome were extracted for expression analyses using the Illumina RNA-seq data for these strains ( Figure 6 ). Fold change values were calculated for both the strains comparing control vs. samples exposed to azamethiphos, resulting in some genes with the most significant differences for both populations, such as Rotatin, Allatostatin-A, Gag-Pol polyprotein, and Ribonuclease H, which were the genes with the most differential values ( Figures 6A, C ). Gene ontology analyses from the genes containing CNVs between susceptible and resistant populations showed that the nucleus was the only enriched cellular component, and various polymerase/catalytic activities were the most enriched molecular functions. The most enriched biological process was DNA metabolic process, followed by the cellular macromolecule metabolic process and DNA integration/recombination process ( Figure 6B ).