Sex-biased gene expression and evolution of candidate reproductive transcripts in adult stages of salmon lice (Lepeophtheirus salmonis)

Adult female L. salmonis were collected from Bay Management Area 2A (BMA2a), Bay of Fundy, New Brunswick, Canada, and returned to the Huntsman Marine Science Center in St. Andrews, New Brunswick. Egg strands were removed and hatched using a seawater hatching system at 10 °C as described by Igboeli et al. (2013). The resulting copepodid larvae were used to infect 27 Atlantic salmon (Salmo salar) smolts obtained from Cooke Aquaculture hatchery in Elmsville, New Brunswick, housed in 10 L re-circulating aquaria at 12 °C ± 0.5 °C and a salinity of 33 ± 2 ppt containing one fish per randomly selected tank. All fish were anesthetized individually using MS-222 (Sigma-Aldrich, Oakville, Ontario, Canada) at a concentration of 0.1 g/L before being placed in 10 L of fresh seawater containing L. salmonis copepodids (30 per fish) for 5 min. Resulting attached lice were allowed to develop to the pre-adult II stage at which point sex could be determined macroscopically. At 22 days post infection (dpi), 10 fish were anesthetized as before, and all pre-adult II male lice (n = 13) were removed to create a virgin female group (n = 12). Six pre-adult II females (n = 6) were also collected from the host during this sampling, flash-frozen, and stored alongside the males at −80 °C for downstream analyses. Once the lice were collected, the anesthetized fish recovered in a salt water bath and then were returned to their original tanks. Another 10 fish (i.e., those infected with both male and female lice) were allocated to the mated-female group (females, n = 15; males n = 11). The remaining seven fish were used to monitor louse development and male-to-female interactions (mate-guarding) to ensure optimal sampling points throughout the study, which was terminated at 36 dpi when all lice had developed to the adult stage. All individuals from the mated-females group had received a spermatophore; similarly, all individuals from the virgin group had an absence of a spermatophore. The presence of egg strands at various developmental stages was observed for the majority of females. Fish were euthanized using 0.25 g/L MS-222 and lice were collected and flash frozen prior to being stored at −80 °C as before.

Total RNA was extracted from each louse individually using TRIzol^®, according to the manufacturers’ instructions (Thermo Fisher Scientific™, Burlington, Ontario, Canada; Chomczynski 1993). An in-solution DNase treatment, TURBO DNA-free DNA removal kit (Ambion™, Foster City, California, USA) was used following the manufacturer’s instructions to remove residual traces of DNA. All samples were analyzed for quantity and purity using a Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientific™, Burlington, Ontario, Canada), whereas a subset of 24 samples was randomly selected for assessing RNA integrity using ExperionStdSens Analysis Kit™ (Bio-Rad Laboratories, Mississauga, Ontario, Canada) as per manufacturer’s instructions. One sample was removed from the study because of potential RNA degradation. In total, 57 RNA extractions were executed, creating a study design with stage, sex, and adult female fertilization status as explanatory variables.

DNase-treated RNA was synthesized to 2 μg of cDNA in an 40 μL solution using iScript™ Reverse Transcription Supermix for RT-qPCR (Bio-Rad) as per manufacturer’s instructions. An identical protocol was used on pooled RNA without the reverse-transcriptase enzyme (no-RT controls) to confirm an absence of genomic contamination. A pool of cDNA containing 3 μL from each sample was prepared for testing primer specificity and efficiency. Primers were designed to specifically target sex-biased sequences from Eichner et al. (2008), Poley et al. (2016), or GenBank sequences related to sexual dimorphism or reproduction (Table 1). A 1%–2% agarose gel electrophoresis procedure using either a TAE or TBE buffer was completed using PCR product to ensure primers amplified a single product from synthesized cDNA before RT-qPCR analysis. Transcript-specific primers, which did not produce a fluorescent band or had multiple bands, were removed from the study. Serial dilutions for 16 genes of interest and two reference genes, elongation factor 1α (EF1α) and ribosomal protein subunit 20 (RPS20), were generated using a 6-step 5-fold dilution series. Efficiency curves were performed in 96 well plates using an Eppendorf Realplex with a 95 °C polymerase activation step for 30 s followed by 40 cycles of 95 °C for denaturation and 60 °C for combined annealing and extension, each lasting 15 s. Three-hundred and eighty-four well plates were loaded using the Aurora Biomed VERSA 10 for RT-qPCR amplification using the Bio-Rad CFX384™ Real-Time System. SsoAdvanced SYBR Green Supermix (Bio-Rad) was used to perform 11 μL qPCR reactions with 1 μL template and a 0.1 μmol/L concentration of each primer pair following the same thermal regime described above. Melt curve analysis was performed by increasing the temperature in 0.5 °C increments every 5 s from 65 °C to 95 °C to corroborate single product formation from gel electrophoresis. Individual C_q values above 37 cycles were considered to be unreliable positives (i.e., inconsistent amplification across replicates; Ct range of –RTs, etc.) and were determined to be non-expressing for this analysis. Reference gene stability was measured using geNorm and mean normalization of RT-qPCR data was completed using qbase-PLUS (Biogazelle, Gent, Belgium) on a log₂ scale relative to the geometric mean of EF1α and RPS20. Individual primer pair efficiencies were used when generating the calibrated normalized relative quantity (CNRQ) values. Gene expression data were visualized using histograms and quantile–quantile plots to confirm normal distribution within each group. In addition, a Levene’s test was used to determine homogeneity of the variances. A t test was performed to determine statistically significant differences between mated and virgin adult female lice groups, and a two-way analysis of variance (ANOVA) was performed with stage and sex as explanatory variables using a significance cut-off of p < 0.05. In cases where a significant result was obtained, Tukey’s HSD test p < 0.05 was used to determine the groups with different expression profiles.

To assess the evolutionary rate of L. salmonis sex-biased genes, putative full-length coding DNA sequences (FLcDNA) were obtained from a BLASTN search against L. salmonis sequences in either the NCBI Non-Redundant (NR) or Transcript Shotgun Assembly (TSA) databases using 90 contig sequences obtained from a previous microarray study on sex-biased gene expression (Poley et al. 2016). Balanced groups containing 30 FLcDNA sequences each were obtained for male-biased, female-biased, and unbiased genes selected based on largest sequence and an open reading frame (ORF) with ≥50% coverage where the transcript had to show concordant sex-bias expression across three populations of lice (Poley et al. 2016). Transcripts with ≥90% identity were selected for further analysis. An ORF finder (ncbi.nlm.nih.gov/orffinder/) was used to predict ORFs of selected transcripts. A TBLASTN using the predicted L. salmonis protein sequences was performed to obtain orthologs in C. rogercresseyi, the closest related organism to L. salmonis with sufficient genomic characterization to run this analysis. Accession numbers for L. salmonis and C. rogercresseyi orthologs can be found in Supplementary Material 2. In order for a pair of L. salmonis and C. rogercresseyi genes to be considered orthologous, results from a TBLASTN required the following criteria be met: E-value ≤10⁻¹⁰, similarity score ≥75%, coverage score ≥50%, and transcript length ≥300 base pairs (Leong et al. 2010). Protein sequences of orthologous pairs were aligned using EMBOSS WATER with default input settings and MARKx10 for FASTA output. The result was used as a template for manual trimming of the nucleotide sequences to match their respective protein alignment. RevTrans (1.4 server) was used to create codon-aware DNA alignments using default parameters, and the FASTA format output was selected (Wernersson and Pedersen 2003). Additionally RevTrans served as a quality assurance check for the manual nucleotide trimming as nucleotide sequences were required to be an exact match to the unaltered protein alignment supplied to RevTrans that acted as a scaffold for the nucleotide alignment.

Alignments from the RevTrans output were used to analyze non-synonymous and synonymous mutation ratios (d_N/d_S) using the YN00 program of the phylogenetic analysis by maximum likelihood (PAML) software (Yang and Nielsen 2000; Yang 2007). Estimated evolutionary rates (d_N/d_S) were calculated for the 30 male-biased, female-biased, and unbiased genes from pairwise comparisons of L. salmonis and C. rogercresseyi. The YN00 program was run using the default settings with following changes; verbose: one detailed output, noisy: ≥9, and weighting—yes. Codon usage bias data for L. salmonis sequences was calculated using the program codonW (version 1.4.2 available from: codonw.sourceforge.net/) with default settings for frequency of optimal codons (F_op). Both d_N/d_S and codon usage bias data were visualized using histograms and quantile–quantile plots and a Levene’s test for homogeneity of variance. One-way ANOVAs were completed for both the evolutionary rate (frequency of d_N/d_S, d_N, and d_S) and codon usage bias (frequency of optimal codons) with male-biased, female-biased, and unbiased genes as explanatory variables using a significance cut-off of p < 0.05. In cases where a significant result was obtained, Tukey’s HSD test (p < 0.05) was used.

All animals were cared for in accordance with the Guide to the Care and Use of Experimental Animals available at: ccac.ca/Documents/Standards/Guidelines/Experimental_Animals_Vol1.pdf. Animal care and use protocols were submitted for review and approval to the UPEI Animal Care Committee (UPEI Animal Care Protocol No. 10-014).

The L. salmonis transcriptome at the pre-adult II stage has been revealed to contain extensive sex-biased gene expression (ca. 35%–43%) (Poley et al. 2016). However, further characterization is required to understand the stability of sex-biased expression throughout different stages. Sex-biased genes were characterized for pre-adult II and adult L. salmonis from BMA2a in the Bay of Fundy, New Brunswick using RT-qPCR for differential gene expression of sex and stage. A total of 16 candidate transcripts (Table 1) were selected from previous salmon louse transcriptomic studies based on their high degree of sex-biased expression, putative involvement in sexually dimorphic properties and were hypothesized to have an influence on mating responsive gene expression. Putative functions include: embryogenesis, embryonic development, vitellogenesis, sperm storage, chorion hardening, feeding response, and semiochemical/pheromone detection (Eichner et al. 2008; Poley et al. 2016). Many of these genes had similar putative functions to known mating responsive genes from other arthropods with homologous sequences and, therefore, we attempted to determine whether this trend was general or taxon-specific (Mack et al. 2006; Rogers et al. 2008).

A principal component analysis (PCA) was performed to visualize the relationships between sex and stage using log₂ CNRQ data for the 16 candidate genes (Fig. 1). The first and second principal components (PCs) explained 36.4% and 24.2% of the total variation in gene expression profiles, respectively. Both sex and stage were separated along PC1 (x-axis), where adult males and adult females were the most distant groups. The majority of genes assayed here were male-biased or female-biased and, therefore, contributed to the separation of samples along PC1. Genes influenced by stage only, unaffected by either sex or stage, and three male-biased genes odourant-binding protein, protein sarah, and vitellogenin-like (Table 1), controlled differences associated with PC2. Some pre-adult II female samples overlapped with male samples suggesting pre-adult II female expression of some genes is more similar to that of males, especially at the corresponding pre-adult II stage. Male samples also had less overall separation between stages compared with females along PC1, indicating stage differences in male expression was smaller than in females.

Four transcripts influenced only by sex included nucleoplasmin-like protein, vitellogenin-like protein, putative odourant-binding protein, and protein sarah (Fig. 2; Table 1). Nucleoplasmin-like protein displayed female-biased expression, whereas the remaining three genes were expressed highest in males. Sex-biased expression, however, was not always observed to be stable throughout the pre-adult II to adult life stages as displayed by the expression profiles of embryo cathepsin L, trypsin-like 4, vitellogenin 1 (LsVit1) and vitellogenin 2 (LsVit2) (Table 1; Fig. 3). All four transcripts were significantly upregulated at the adult female stage when compared with other sex and stage combinations, with embryo cathepsin L, vitellogenin 1, and vitellogenin 2 being expressed significantly lower in pre-adult II females when compared with pre-adult II males (p < 0.001; LsVit1 had no detectable expression (n.d.)) and adult males (p < 0.001) (Fig. 3). Conversely, Trypsin-4 was expressed highest in adult females followed by pre-adult II females and then males (p = 0.005) demonstrating a conservation of sex-biased expression that is further amplified across stage. Vitellogenins are precursors of egg yolk glycoproteins produced within fat bodies for many arthropods and have vital roles in progeny fitness (Schneider 1996; Eichner et al. 2008; Dalvin et al. 2009; Dalvin et al. 2011). In L. salmonis, vitellogenins (LsVit1 and LsVit2) are expressed throughout the adult female subcuticular tissue and are involved in embryogenesis. Additionally, these proteins provide an essential energy source throughout the planktonic larval life stages (Eichner et al. 2008; Dalvin et al. 2011; Edvardsen et al. 2014). Expression knock-down of a separate novel female specific L. salmonis yolk-associated protein (LsYAP) has been shown to disrupt proper embryonic development, highlighting their importance to progeny development and fitness (Dalvin et al. 2009). These results are consistent with a study by Eichner et al. (2008) who observed a gradual increase in expression of LsVit1 and LsVit2 from the PAII to adult females as well as six time points across adult stage development.

A total of six transcripts had no detectable expression in ≥2/3 of individuals from one sex across both stages. These transcripts (C-type lectin, chorion peroxidase, kunitz 1, and peroxidase) were, therefore, considered to be male-specific, whereas carboxypeptidase and vitellogenin receptor (Table 1; Table 2) were similarly classified as female-specific (Table 2). Although some individuals from the putatively non-expressing sex formed PCR products, all Cq values were high (between 34 and 39) and (or) melt curves were unreliable potentially representing artifacts rather than true expression from the opposite sex. Carboxypeptidase and vitellogenin receptor were classified as female-specific in expression and may be involved in two essential processes for female lice. Vitellogenin receptors are another egg-yolk-associated protein that facilitates the endocytosis of vitellogenins into developing oocytes (Schneider 1996; Tufail and Takeda 2009). As mentioned, successful incorporation of vitellogenins are essential to proper development of embryos (Dalvin et al. 2009) and disruption of this vital pathway may elucidate novel reproduction-associated louse management strategies. Carboxypeptidases are putatively involved in feeding and managing blood meals within L. salmonis similar to other insects and parasitic arthropods (Edwards et al. 1997; Warr et al. 2007; Isoe et al. 2009). Carboxypeptidase B proteins have been investigated in Anopheles gambiae as targets for development of a transmission-blocking vaccine for Plasmodium falciparum (Lavazec et al. 2007). An interesting secondary effect reported by Lavazec et al. (2007) was mosquitos having fed on anti-carboxypeptidase mice exhibited reduced reproductive capacity. Furthermore, a recombinant vaccine against a gut-associated glycoprotein BM86 in Boophilus microplus has been shown to be effective in controlling tick infestations (Willadsen et al. 1995). Female-specific carboxypeptidases and other gut/digestion associated proteins may provide interesting targets for further investigation to determine their viability as vaccination options in L. salmonis.

C-type lectin, chorion peroxidase, kunitz-1, and peroxidase are putative seminal fluid proteins (SFPs) in salmon lice providing further evidence for male-specific expression profiles. SFPs are the products of male accessory glands and ejaculatory ducts which are transferred along with sperm in the male ejaculate to a female recipient (Avila et al. 2011). SFPs have been shown to impact multiple processes such as egg production, female longevity, female behaviour, sperm storage, mating plug formation, and receptivity to re-mating across a wide range of arthropod species (Chapman 2001; Avila et al. 2011). Recently, a relevant new gene family has been identified in L. salmonis essential for proper development of the wall surrounding the spermatophore (Borchel and Nilsen 2018). Understanding the evolution of SFPs and how they elicit responses in females may provide practical implications for pest management strategies. Increased selection, or genetic engineering, of males producing SFPs that have greater association with negative impacts on female fitness such as reduced receptivity to re-mating or male sterilization may provide attractive avenues for alternative management strategies (Sirot et al. 2015). These results provided confirmation of the stability of sex-biased gene expression from the pre-adult II to adult stage of L. salmonis, but also determined transcripts may have differing expression profiles dependent on stage, similar to observations by Eichner et al. (2018) in earlier life stages (chalimus II and pre-adult I). The results provide further characterization of sexually dimorphic characteristics of L. salmonis and provide novel insights into sex-biased expression during the final life stages.

Another goal, while exploring the stability of sex-biased expression across stage, was to investigate the impact copulation/reception of a spermatophore had on female transcripts. Mated and virgin females were analyzed for differential expression for all 16 candidate transcripts, but no significant expression differences were observed between copulatory statuses (t test; p > 0.05 for all transcripts). The transfer of male SFPs are known to induce gene expression differences in Drosophila melanogaster and A. gambiae females (Soller et al. 1997; Peng et al. 2005; Mack et al. 2006; Ram and Wolfner 2007; Rogers et al. 2008; Avila et al. 2011). Reception of a spermatophore in D. melanogaster facilitates differential gene expression either transiently (1–3/4–6 h) or for longer periods (≥12 h) (McGraw et al. 2008). Candidate L. salmonis genes with similar functions to known short- and long-term mating responsive genes such as chorion protein 36/38, odourant binding protein 99a, yolk protein 1/2/3 and vitelline membrane 26Aa/32E/34Ca appear not to be impacted post-reception of a spermatophore in L. salmonis (McGraw et al. 2008). These results suggest, similar to a study by Eichner et al. (2008) that copulatory status does not appear to be a trigger for egg production despite having similar functions to genes that are sensitive to copulatory status in other arthropods.

Potential explanations for the lack of differential expression in mated female L. salmonis include, but are not limited to: (1) select candidate genes with similar functions to known mating-responsive genes within Arthropoda are potentially not mating-responsive in L. salmonis (Soller et al. 1997; McGraw et al. 2008; Rogers et al. 2008). (2) Mating-responsive gene expression is known to be transient within Arthropoda (McGraw et al. 2008; Rogers et al. 2008). Because of the limitation of lice numbers post infection, in addition to the difficulty monitoring copulatory status for each mating-pair, inclusion of multiple sampling time points post-mating wasn’t feasible. It is therefore possible the single sampling time point was either too early or late to capture differential expression. (3) Despite having mated and virgin lice isolated in separate tanks, it is possible the recirculation system was unable to remove semiochemicals produced by males from within the mated group tanks. If semiochemical cues are important to louse receptivity to mating then the virgin lice may have falsely perceived the presence of a male and, therefore, potentially influenced their reproductive mechanisms. (4) Lastly, the use of whole homogenates over specific tissues, more specifically the genital segment, may have resulted in a dilution effect of smaller scale expression differences as the greatest number of differentially expressed genes post-mating understandably appear in the reproductive tissues (Mack et al. 2006; Edvardsen et al. 2014).

Continued investigation into mating-responsive gene expression may provide insights into alternative management strategies that are becoming increasingly necessary to be effective and sustainable when managing salmon lice infestations. Reproduction targeted or disrupting strategies are not novel ideas and have been studied for decades (Cardé and Minks 1995). Synthetic semiochemicals provide an avenue to prevent or reduce mating success as a means of control instead of drug treatments susceptible to the acquisition of resistance (Smart et al. 2014). Furthermore, hormone disruption has been proposed as a viable alternative strategy to limit Anopheles spp. reproductive success and, thus, the spread of Plasmodium faliciparum (Childs et al. 2016). A recent study used the CRISPR-Cas9 system to distort the sex-ratio of A. gambiae by utilizing a specific ribosomal sequence targeted Cas9 approach during spermatogenesis to produce an extreme male-bias in the progeny (Galizi et al. 2016). Further efforts into the identification of mating-responsive gene expression, especially through means of whole transcriptome analyses of mated vs. virgin females will be beneficial to markers of insemination success, progeny development, and provide similar targets as described above for alternative management strategies.

The ratios of non-synonymous and synonymous substitutions (d_N/d_S) were calculated for 30 male-biased, female-biased, and unbiased transcripts identified in a previous transcriptomic study (Poley et al. 2016). Male-biased L. salmonis transcripts had significantly higher d_N/d_S values when compared with unbiased transcripts, whereas values for female-biased genes were between those of the two groups (Table 3; Fig. 4). This suggests that L. salmonis sex-biased genes, especially those that are male-biased, are under reduced selective constraint or greater positive selection (Yang and Nielsen 2000; Zhang et al. 2004; Ellegren and Parsch 2007). Analyses using d_N/d_S require a closely related species for sequence comparisons and C. rogercresseyi was selected as it is the closest related species to L. salmonis with sufficient genomic and transcriptomic information. It is important to note, however, that these species have been diverging for ~100–150 million years (Eyun 2017), which has allowed for greater accumulation of synonymous mutations over time, thereby reducing the overall sensitivity of the d_N/d_S analysis. When examining the d_N and d_S variables separately, it was observed that d_S values for male-biased, female-biased, and unbiased genes were not significantly different and had similar saturation of synonymous mutations (d_S). This saturation caused resulting d_N/d_S values to be underestimated because of the stable and consistently high d_S. Similar to a study by Zhang et al. (2004) who also describe a saturation of d_S even within Drosophila spp., we can consider only d_N as it is the driving force on the overall d_N/d_S differences. Male-biased transcripts had significantly greater d_N values when compared with both female-biased and unbiased transcripts (Table 3; Fig. 4). This result supports the previous analysis by Poley et al. (2016), who showed a greater proportion of orphan genes for those with male-biased expression and is consistent with results found in Drosophila spp. (Zhang et al. 2004).

Codon usage bias was also investigated for these gene categories to validate the d_N/d_S results. Codon usage bias can be described as the weak selection pressure that drives the accumulation of a specific codon over other synonymous codons (Akashi 1994; Hambuch and Parsch 2005; Behura and Severson 2013). Often these optimal codons reflect an organism’s tRNA pool and provide the benefits of increased translation efficiency and accuracy (Duret 2000). Codon usage bias is often measured by the frequency of optimal codons (F_op) where highly expressed and (or) heavily conserved genes have the greatest accumulation of optimal codons (Powell and Moriyama 1997). It was therefore expected that this weak selection pressure would have a lesser effect on L. salmonis genes with higher d_N/d_S values. Male-biased L. salmonis transcripts indeed had the lowest frequency of optimal codons in this analysis (Table 3), which is again consistent with results from D. melanogaster (Hambuch and Parsch 2005) and d_N/d_S values presented above.

Overall, these analyses indicate that male-biased L. salmonis transcripts likely follow the commonly observed trend of rapid evolution when compared with unbiased genes (Ellegren and Parsch 2007). Whether this is caused by a reduction in selective constraint or increased sexual, natural, and (or) antagonistic selection pressures is not yet fully understood (Zhang et al. 2004; Ellegren and Parsch 2007; Parsch and Ellegren 2013). The annotations of male-biased transcripts in L. salmonis are enriched for putative reproductive functions (Poley et al. 2016; Table 1) that are likely influenced by strong selective pressures, which is consistent with observed patterns of rapid evolution for reproduction associated genes (Swanson and Vacquier 2002; Haerty et al. 2007; Meisel 2011; Grath and Parsch 2012). Recently an interesting sex-associated single nucleotide polymorphism (SNP) marker was identified in salmon lice suggesting they operate using a ZW sex determining system (Carmichael et al. 2013). Global dosage compensation is a major factor in sex-biased gene expression as the homogametic sex of a given species has twice the number of available gene copies as the opposite sex for those genes present on sex chromosomes, for example D. melanogaster compensate by increasing the transcription of the single male X chromosome two-fold (Baker et al. 1994). Compensation mechanisms appear to be less effective (or selective for only a proportion of genes on the sex chromosome) in ZW systems than XY (Itoh et al. 2007; Mank 2009; Wolf and Bryk 2011). Determining the nature of L. salmonis dosage compensation will help elucidate the necessity and (or) consequences of sex-linked gene expression.

Our analyses have provided a useful starting point to better understand evolutionary aspects of male-biased, female-biased, and unbiased genes; however, they do lack some degree of specificity as transcripts were selected based on sex-biased expression and transcript length, not function. Future analyses targeting sex-biased expression with specific roles and expression breadths (i.e., reproductive genes in testes and ovaries) or sex-linked genes will provide greater insight into these mechanisms. Additionally, optimization and automation of evolutionary rate analyses for sex-specific whole transcriptome data would be useful to globally confirm elevated sex-biased gene evolution in L. salmonis. To our knowledge this is the first time these techniques have been applied to salmon lice and will serve as valuable tools for future research. These techniques can be applied to any gene groups for comparative purposes; for example, those involved in the host-parasite interaction, which are likely under similar elevated selection pressures when compared with somatic genes.