Introduction

Antarctic krill (Euphausia superba, hereafter called krill), is a keystone species in the Antarctic marine ecosystem. Krill have a circumantarctic distribution in the Southern Ocean and are one of the world’s most abundant animals with an estimated total biomass of 379 million tonnes (Atkinson et al. 2009; Siegel and Watkins 2016). As filter feeders they hold a critical role in the marine Antarctic food web by directly linking primary producers (phytoplankton) to various predators, including baleen whales, penguins, seals, fish and seabirds (Hopkins et al. 1993). Krill also present potentially the largest underexploited source of wild animal protein for human consumption (Tou et al. 2007) and their behaviour of aggregating in large, dense swarms often containing billions of individuals (Nicol 2006) make them an attractive target for the currently expanding krill fishery (Nicol et al. 2012).

The massive abundance of krill has driven significant research effort on this species, but many aspects of krill biology remain elusive. One outstanding question in krill population biology is why sex ratios in wild krill swarms are often biased towards one or the other sex (Quetin and Ross 1984; Watkins et al. 1986, Virtue et al. 1996). For example, female proportion of twelve large krill swarms ranged from 7 to 82% on an Australian research voyage off East Antarctica carried out in 2016 (Kawaguchi, personal communication). Distorted sex ratios directly affect effective population size and have important implications for population dynamics, genetic diversity and reproductive behaviour (Frankham 1995; Mathews 2002). Similarly, virtually nothing is known about sexual development and differentiation of krill on a molecular level.

Krill larval and sexual development on a morphological level has been studied for decades and is relatively well understood (Jia et al. 2014; Kawaguchi 2016). Internal development of the gonads may begin as early as the late larval or early juvenile stage, although at this stage gonads have not differentiated into identifiable male or female versions (Cuzin-Roudy 1987; Kawaguchi 2016). Primary sexual differentiation of the gonads into incipient ovaries and testes may occur in juvenile krill of > 23 mm size (Cuzin-Roudy 1987), although an older study claims to be able to determine sex of individuals down to 10 mm in size (Bargmann 1945). Regardless of the precise lower limit of morphological identification, sex determination through dissection of such small krill is extremely difficult and results are often ambiguous. Unambiguous classification of sex becomes possible with the gradual appearance of external sexual characteristics from the late sub-adult stage onwards. Sexually mature adult males display a claw-like structure called petasma on the first two pairs of pleopods, which are presumably used to transfer spermatophores onto females, whereas mature females have a thelycum, a structure to which the spermatophores are attached (Kawaguchi 2016; Siegel 2016). After the reproductive season in the austral summer, both internal and external sexual structures regress over winter, sometimes to the extent that no external sexual characteristics are observable (Kawaguchi et al. 2007; Thomas and Ikeda 1987), followed by re-maturation during the next reproductive season.

The primary mechanism of sex determination in krill is not known. Sex determination mechanisms could theoretically range from genetic to environmental to parasitic or any combination of these (Duron et al. 2008; Kato et al. 2011; Carmichael et al. 2013; Bachtrog et al. 2014). Furthermore, sex determination mechanisms can evolve rapidly, and may even differ between populations of the same species (Pen et al. 2010). Similar to other Euphausiacea species, krill probably do not have heteromorphic sex chromosomes (Thiriot-Quiévreux and Cuzin-Roudy 1995; Thiriot-Quiévreux et al. 1998), but a sex determination region may still exist within the krill genome. However, the krill genome is very large with an estimated size of 47 gigabases (Jeffery 2011; Deagle et al. 2015) and contains many repetitive elements (Deagle et al. 2015), which are the primary reasons that genome sequencing and assembly has not been carried out (Jarman and Deagle 2016). Initially, we pursued multiple avenues to develop a simple DNA-based test with the aim to determine the sex of krill at any developmental stage, as well as to potentially determine sex ratios in mixed bulk samples of krill containing both sexes. This included trying to identify sex-specific markers in an available restriction site-associated DNA sequencing (RAD-seq) dataset (Deagle et al. 2015) and examining transcriptome data for sex-linked SNPs (for a summary of this work see Online Resource 1). However, due to the extraordinary complexity of the krill genome we were unable to find, or completely discount, the presence of DNA markers so we ultimately shifted our efforts to gene expression differences described in the current paper.

In general, sex differentiation on a molecular level is triggered by a primary signal that differs between males and females (Bachtrog et al. 2014). Irrespective of the underlying mechanism (genetic, environmental or a combination), the primary signal usually triggers differential expression of one or a few key genes. This results in a whole cascade or network of differential gene expression, ultimately leading to different sexual phenotypes (Bachtrog et al. 2014). This differential gene expression network could be observed in krill even if the primary signal remains unknown. Furthermore, by studying krill samples from early developmental stages (furcilial larvae) through to sexually mature and regressed stages, our understanding of initiation, progress and regression of sexual development at a molecular level will be improved. Additionally, for any developmental stage that displays sexual differentiation on a molecular level, a differential gene expression test could be developed to identify the sex of krill.

The first comprehensive krill transcriptome became available in 2015 (Meyer et al. 2015), and in 2017 two additional online krill transcriptome databases were added (the “KrillDB” available at krilldb.bio.unipd.it (Sales et al. 2017) and the “SuperbaSE” available at www.krill.le.ac.uk (Hunt et al. 2017)). These resources could assist in the identification of sex-specific gene expression differences by providing a catalogue of genes present in the genome.

The aim of this study was to develop a morphology-independent genetic test to determine sex of krill as a tool for future studies and to gain insights into sexual development of krill on a molecular level. Specifically, we addressed the following questions:

  1. 1.

    Can sex in krill be determined from gene expression differences, even in the absence of sex-specific external morphological features?

  2. 2.

    At what point in krill development does molecular sex differentiation begin, and is morphological regression associated with a decline in molecular sex differentiation?

Materials and methods

Sex-specific transcriptome analysis

To obtain sex-specific transcriptomes, RNA was extracted from pooled gonad tissues of five sexually mature males and females, respectively, originating from the krill aquarium at the Australian Antarctic Division (AAD) in Tasmania, Australia (Kawaguchi et al. 2010). For RNA extractions, RNeasy mini kits (QIAGEN GmbH, Hilden, Germany) were used and total RNA (8 µg each) was sent to Geneworks, South Australia (www.geneworks.com.au), for Illumina TruSeq 75 bp paired-end sequencing. Both the male and the female samples were sequenced in two technical replica (four groups in total) on two separate lanes.

Transcriptome sequences were aligned to the National Centre for Biotechnology Information’s (NCBI) nucleotide collection database to assess contamination levels, which were found to be low. They were then aligned to the 2015 version of the krill transcriptome (Meyer et al. 2015) using RSEM, a software tool to accurately quantify transcript abundance in the absence of a reference genome (Li and Dewey 2011).

We used the two technical replica of the gonad transcriptomic data to estimate expression differences between testes and ovaries using the function exactTest in R package edgeR (Robinson et al. 2010). While the lack of biological replication prevented us from performing meaningful statistical analyses, the comparison between testis and ovary transcriptomes nevertheless provided a starting point to identify differentially expressed genes in whole krill. The most differentially expressed transcripts (smallest Bonferroni-corrected p values, compare to Online Resource 2) were manually scanned for genes known to be involved in sex determination or differentiation from the literature on other Crustacea for quantitative PCR (qPCR) verification of candidate genes in whole krill.

Krill samples used to verify sex-specific marker gene expression

To verify sex-biased expression of selected transcripts, 16 sexually mature krill “Test” samples from the krill aquarium at the AAD were used (Table 1). The sex of each individual was identified by determining the presence of either a thelycum (females) or a petasma (males) (Kawaguchi 2016). Once differentially expressed genes were identified using these Test krill (see below), 86 additional krill samples of various developmental stages were analysed (Table 1). This included krill from the AAD aquarium as well as freshly collected krill from two voyages covering vast areas of East Antarctica on the RSV Aurora Australis (collected using permits AMLR 08-12-2337 and AMLR 12-17-4037). Trawl coordinates of wild caught krill samples as well as original trawl coordinates of captive krill are shown in Fig. 1 and listed in Online Resource 3. In 40 of these additional 86 samples, the sex could be determined based on external morphology. These samples included sexually mature adults, regressed adults and sub-adult krill (Table 1). The regressed krill were from the AAD aquarium and were either incipient regressed krill (“ReI”) collected under “June” light conditions, or advanced regressed krill (“ReA”) collected under “August” light conditions (simulating early and advanced austral winter conditions, respectively). The remaining 46 samples were juvenile and larval krill (furcilia stages III, IV and VI, determined following Jia et al. (2014)). These developmental stages do not display external sexual dimorphisms and therefore the sex of these samples could not be determined based on external morphology.

Table 1 Details of whole krill samples used to verify gene expression differences using quantitative PCR
Fig. 1
figure 1

Origin of krill samples used in this study. Grey diamonds denote samples that were collected in the wild and immediately stored in RNAlater, whereas grey dots represent original collection coordinates of krill samples that were transported to and kept alive in the krill aquarium in Kingston, Tasmania, where they were later sampled for this study. Shown in black are Australian research stations. The maps were drawn using the R package “maps”, function “maps”

Nucleic acid extraction and qPCR setup

After morphological sex determination of live krill (where possible), Test samples were stored in AllProtect Tissue Reagent (QIAGEN), whereas all other samples were stored in an RNA preserving solution similar to RNAlater (Wit et al. 2012). Samples were then kept at − 20 °C until RNA extraction. For RNA extraction, Test samples were halved lengthwise, whereas all other samples were processed whole. RNA was extracted using Maxwell 16 LEV simplyRNA Tissue Kits (Promega Corporation, Madison, WI, USA) following manufacturer’s instructions with the exception that the volume of homogenization solution used for blending the samples was increased depending on tissue size (details in Online Resource 3). For all samples, 200 µL of homogenate was used for RNA extraction with 10 µL of DNase. RNA was quantified using Qubit RNA BR assay kits with a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and RNA quality was assessed using Agilent RNA 6000 Nano kits with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Note that the electropherograms obtained from the Bioanalyzer showed a large (triple) rRNA peak of 18S size, whereas the 28S rRNA peak usually found in mammals was very weak or absent (Online Resource 4). RNA was stored at − 80 °C until converted to complementary DNA (cDNA) using QuantiTect Reverse Transcription Kits (QIAGEN). For cDNA conversion, 1 µg of RNA in a total of 20 µL reaction volume was used for all samples except larval samples with low RNA yield. Controls without reverse transcriptase were included and no genomic DNA contamination of samples were detected. Converted cDNA was diluted to obtain a final DNA concentration of approximately 10 ng µL−1 (details in Online Resource 3).

Genomic DNA (gDNA) was extracted from 16 Test samples to verify the presence of sex-biased genes (see results) in the genome of both sexes. DNA was extracted with 200 µL of the same homogenate used for the RNA extractions using Maxwell 16 Tissue DNA Purification Kits (Promega Corporation, Madison, WI, USA). DNA was quantified using Qubit DNA BR assay kits with a Qubit 2.0 Fluorometer and diluted to 10 ng µL−1 for PCR amplifications. PCR reaction mixes and conditions were identical to cDNA conditions described below, with the exception that 1 µL of gDNA template (10 ng µL−1) was used in a total reaction volume of 10 µL.

Initially, four male and four female-biased transcripts were selected from the gonad-specific transcriptomic data to determine their sex bias in whole krill using qPCR. As the krill genome contains vast amounts of duplicated sequences (Deagle et al. 2015), five primer pairs were designed for each of the eight amplicons using primer3 (Koressaar and Remm 2007; Untergasser et al. 2012), and each primer pair was then aligned to the krill transcriptome using local BLAST (Madden 2013) to identify a suitable primer pair that exclusively amplified the desired amplicon, or near-identical homologues (primer sequences were separated by 30 “N” bases, word size was set to 7 and e-value to 0.1). Sex bias of the selected amplicon was initially tested using cDNA of the eight male and eight female Test krill samples. When none of the initial male candidates showed sex-biased gene expression in whole krill samples, 10 additional potentially male-biased fragments were selected and primer pairs were designed in the same way.

Primer pairs for eight reference gene candidates were kindly provided by Alberto Biscontin and Cristiano de Pittà from the University of Padova (Biscontin et al. 2016). The primer pairs were tested in the 16 Test samples for stable expression across sexes (data not shown), and the following three reference genes were chosen: 28S rRNA; myotubularin-related protein 4 (MTMR4), ubiquitin carboxyl-terminal hydrolase 46 (USP46) (Table 2; Online Resource 5). None of these three genes displayed statistically significant sex-biased expression (data not shown).

Table 2 Primer sequences used for qPCR analysis of sex biased and reference transcripts. Gene IDs are from the two krill transcriptomes available online (IDs starting with ESG refer to the 2017 krillDB at krilldb.bio.unipd.it; IDs starting with KEST refer to the 2017 SuperbaSE at www.krill.le.ac.uk). Primer efficiencies and R2 calculated from standard curves are also provided

For cDNA fragment amplifications, the PCR reaction mix contained 0.1 µL each of 10 µM forward and reverse primer; 0.1 µL of bovine serum albumin (BSA, 20 mg ml−1), 0.5 µL EvaGreen dye, 20 × in water (Biotium, Hayward, CA, USA); 5 µL 2 × Phusion High-Fidelity PCR Master Mix with HF Buffer (New England Biolabs, Ipswich, MA, USA) and 2 µL of cDNA template (10 ng µL−1) in a total reaction volume of 10 µL. PCRs were done on a LightCycler 480 II (Roche Diagnostics Australia Pty Limited, North Ryde, NSW, Australia) and thermal cycling conditions were 98 °C for 2 min followed by 40 cycles of 98 °C for 5 s, 58 °C for 20 s and 72 °C for 20 s, followed by a melt curve from 50 to 99 °C at a ramp rate of 2.2 °C s−1 with five acquisitions per degree.

Quantitative PCR data analysis

Crossing point (Cp) values were calculated on the LightCycler 480 II using the “Abs Quant/2nd Derivative Max” method with standard settings. Primer efficiencies and R2 values were calculated with standard curves for each sex biased and reference transcript (Table 2). Data were analyzed using qbase + version 3.1 (Biogazelle, Zwijnaarde, Belgium). For male-biased transcripts, interrun calibrations included two female and one male sample, whereas for the female-biased transcripts, only the two female samples could be used, as in the male sample these genes were not sufficiently expressed. To find the most stable reference genes, GeNorm implemented in qbase + was used. Subsequently, MTMR4 and USP46 were used as reference genes to calibrate gene expression of target genes, while 28S was excluded from all analyses due to trawl-specific expression bias. After calibration using the two reference genes, target gene expression was normalized against the mean expression of each target gene across all samples. Therefore, values above 1 correspond to samples with above average gene expression, whereas values between 0 and 1 correspond to samples with below average gene expression. Expression values were exported without logarithm transformation and subsequently analyzed using R (R Core Team 2017). Two PCR technical replica were included for two female and one male sample and found to be very consistent. Mean values of technical replica were reported.

As many of the male samples did not express the female-biased transcripts, the expression bias between sexes could not be calculated in a statistically meaningful way. However, to estimate the extent of the expression bias, Cp values for samples that did not express a transcript were set to the maximum of 35, and log10 transformed expression values were plotted with mean and standard errors separately for each sex or grouped for sex and origin.

For male-biased transcripts, the expression bias between sexes was calculated for all morphologically sexed samples together using a linear model with “sex” as explanatory variable of log10 transformed relative gene expressions, and the function lsmeans (package lsmeans (Lenth 2016)). To address potential differences between different sample types, we used a similar linear model with “sex” and “origin” interacting as explanatory variables, with pairwise comparisons of sex nested in origin and p values adjusted for multiple testing using the method “holm” (Holm 1979).

To classify males and females on a molecular level, gene expression data of the three female-biased genes of all krill with known sex (sub-adult, adult and regressed krill) was used in a Discriminant Function Analysis (R package MASS (Venables and Ripley 2002), function lsa). Juvenile and larval samples were then added to this analysis to predict the sex of the individual samples. This analysis was then repeated using the two male-biased genes.

Results

Transcriptome analyses to identify differentially expressed genes

Comparing the two technical replica of transcriptomic data from the ovary and testis revealed that of the 26,268 transcripts aligned to the 2015 krill transcriptome (Meyer et al. 2015), about three quarters were differentially expressed between the two tissues (10,492 (40%) testis biased, 8986 (34%) ovary biased; p < 0.05 after Bonferroni correction for multiple testing, see Online Resource 6 and Online Resource 2). Raw sequences from both tissues (ovaries: 29,447,654 and 18,233,515 reads; testes: 8,120,993 and 10,465,586 reads) have been deposited in the Australian Antarctic Data Centre (https://dx.doi.org/doi:10.26179/5cd3c8fec9ad8).

Female-biased gene expression

The most differentially expressed transcripts from the gonad transcriptomic data were manually scanned for genes involved in sex determination or differentiation known from the literature. Three transcripts showed female-specific expression in krill samples that could be sexed from external morphological features (Fig. 2) and were annotated to vitellogenin receptor (VGR), nuclear autoantigenic sperm protein (NASP) and polehole-like protein (PHL) (transcript sequences are in Online Resource 5, primer sequences are listed in Table 2). In the gonad transcriptomic data, VGR was the most strongly differentially expressed female-biased transcript, whereas differential expression of NASP and PHL was slightly less pronounced (see Online Resource 2). Interestingly, NASP was present in over 50 highly female-biased variants in the gonad transcriptome, thus overall the expression difference of this gene was remarkable. Primer sequences for transcripts with sex-biased expression are listed in Table 2, including the equivalent gene IDs of the two online available krill transcriptomes.

Fig. 2
figure 2

Relative expression of five sex-biased genes in krill samples where sex could be determined from external morphology (n = 56, i.e. 16 Test samples plus additional 40 samples, see Table 1), separately for females (pink) and males (blue) (mean ± standard error). Gene expression was calibrated using two reference genes, normalised against average expression for each gene, and log10 transformed, i.e. positive values show above average expression, whereas negative values show below average expression. Female-biased genes were not expressed in some male samples and expression values were artificially set to a minimum detectable value (marked with †). In these cases, no meaningful statistical analyses could be performed. For male-biased genes, significant effect of sex is marked with asterisks (***p < 0.001)

These three genes were consistently strongly expressed in females across all krill samples that could be sexed from external morphological features, including sub-adult, sexually mature adult and two stages of regressed adult krill (Table 1; Fig. 3a). In contrast, gene expression of these genes in males was at or near the detection limits, and many samples had to be manually set to the lowest possible values for presentation in Figs. 2 and 3a. As the majority of male samples did not express these genes, we could not quantify the gene expression difference between sexes, however, it is clearly very large. PCR amplification of these three transcripts in gDNA from the 16 Test samples revealed that these three genes are present in the genomes of both male and female krill. This suggests that the observed female-specific gene expression is not due to the presence of these genes on heteromorphic sex chromosomes.

Fig. 3
figure 3

a Relative expression of three female-biased genes in male (blue) and female (pink) krill, collected from different geographic locations and at different developmental stages (see Table 1; values show mean ± standard error for krill of each sex in each location/stage similar to Fig. 2). Sex of juvenile krill (TO2) could be determined based on a discriminant function analysis (DFA) (Fig. 4a). Furcilia stages III, IV and VI (shown together in grey) did not express these genes and sex could not be determined. † Some samples did not express this gene. ‡ None of the samples in this group expressed this gene. b Relative expression of two male-biased genes in the same groups of male and female krill (excluding furcilia stages). Sex of juveniles was determined based on the DFA of the female-biased genes. Significant effects of sex within origin are marked with asterisks (*p < 0.05; **p < 0.01; ***p < 0.001) indicating these genes are not universally elevated in male krill

Discriminant function analysis (DFA) of the three genes showed strong separation of the sexes (Fig. 4a). When juvenile samples with unknown sex were added to the DFA to predict their sex, samples clearly grouped either with one or the other sex with a posterior probability of 1.00 (Fig. 4a). The juvenile females, as classified by this assay all expressed the female-specific genes at similar levels to the developmentally more advanced females, whereas expression in juvenile males was similarly low or absent as in the more advanced males (Fig. 3a). Furcilia larvae of developmental stages III, IV and VI (10 samples each) clustered with male samples in the DFA (data not shown), i.e. none of the larvae expressed these female-biased genes (Fig. 3a). Thus, sex could be determined unambiguously from early juvenile stages onwards, including sub-adult and sexually regressed krill, based on the expression of these three female-specific genes together with two reference genes (Fig. 5).

Fig. 4
figure 4

Discriminant function analyses of krill with known sex (pink: females, blue: males) using the expression of three female-specific genes (a) and two male-biased genes (b). Using these analyses, the sex of juvenile krill was predicted (shown in grey)

Fig. 5
figure 5

Schematic of methods available to determine the sex of krill. ✗: no sex determination possible. Small ✓: sex determination sometimes possible, but results may be ambiguous. Medium ✓: sex determination often possible, with a level of uncertainty. Large ✓: confident and unambiguous sex determination possible. Drawings of telsons indicating furcilia larval stages were adapted from Fraser (1937); drawings of male pleopods showing the development of the male petasma were adapted from Bargmann (1937)

Male-biased gene expression

A total of 14 male-biased candidate transcripts were investigated based on the gonad transcriptomic data. The only two transcripts that showed moderate male-biased expression in krill samples that could be sexed from external morphological features were annotated to an ankyrin repeat protein (ANK) and aquaporin 12A (AQP12A) (Fig. 2; transcript sequences are in Online Resource 5; primer sequences are listed in Table 2). In the pooled gonad transcriptomic data, the expression of both genes was strongly male-biased, although they were not among the most differentially expressed genes (see Online Resource 2). When differential gene expression between sexes was calculated for different sample origins, ANK was differentially expressed in adult krill samples both from the wild and from the aquarium as well as in regressed krill (Fig. 3b). However, there were substantial interactions between the effect of sex and of sample origin, with the result that expression levels in males from one origin were similar to that of females from another, while within each origin sex differences were nevertheless substantial (Fig. 3b). ANK expression levels did not differ between sexes of sub-adult krill that had been sexed morphologically. AQP12A was differentially expressed in krill originating from the krill aquarium, including adult and regressed krill. Sexed krill collected in the wild did not show significant expression differences for AQP12A, including sub-adult and adult krill.

A DFA combining the expression data of ANK and AQP12A for all 56 sexed individuals could not clearly differentiate between sexes (Fig. 4b) and the sex of juvenile samples could not be predicted when the samples were added to this DFA (Fig. 4b). Based on juvenile sex classification calculated from the three female-specific genes, we found that juvenile expression for both male-biased genes did not differ between sexes (Fig. 3b). Thus, any developmental stages before full sexual maturity could not be differentiated based on the expression of these two genes, and confidence for sex classification of sexually mature and regressed krill must be considered low (Fig. 5).

Discussion

We identified substantial gene expression differences between gonads of male and female krill (ovaries vs testes transcriptomes). This result cannot be interpreted strictly as sex-specific gene expression differences, as gene expression between any different tissues usually varies substantially (Parisi et al. 2004). To get a more accurate estimate of global gene expression differences between male and female krill, transcriptomes would need to be produced from RNA extracted from whole krill samples, including biological replication to allow for statistical testing. Nevertheless, our data provided us with a starting point to identify sex-biased gene expression in whole krill and to develop a genetic sex determination test.

We identified three genes that were exclusively expressed in female whole krill throughout development from early juvenile stages onwards. This allows for a molecular test that can unambiguously determine the sex of krill, even if no sex-specific external morphological features are apparent. While sex determination in the absence of external morphological features may be possible to some degree through dissection of the gonads under the microscope (Bargmann 1945; Cuzin-Roudy 1987), it requires extensive microscopy and dissection experience, particularly for small krill, and may not always result in unambiguous classification of sex. Our molecular test, on the other hand, uses RNA extracted from whole krill and is a generic method that delivers clear classification. The test allows studies for the first time to include the effect of sex in juvenile or strongly regressed krill. One example application would be to examine whether juvenile and sexually regressed krill swarms also display skewed sex ratios, similar to sexually mature adult krill swarms (Quetin and Ross 1984; Virtue et al. 1996), which will improve our understanding of krill population dynamics and life history events.

The three female-specific genes we identified for distinguishing krill sex all hold an important role in female sexual development. VGR transports vitellogenin, a yolk precursor, into the developing oocytes, and is thus fundamentally involved in oocyce growth and maturation (Tiu et al. 2008; Subramoniam 2010; Roth and Khalaila 2012). Ovarian-specific expression of VGR has been found in several crustacean species, e.g. the giant tiger shrimp Penaeus monodon (Tiu et al. 2008), two river prawn species, Macrobrachium rosenbergii and M. nipponense (Roth and Khalaila 2012; Bai et al. 2016) and the red swamp crayfish Procambarus clarkii (Shen et al. 2014).

NASP is a histone chaperone involved in histone transport pathways, required for cell proliferation and widely distributed throughout eukaryotes (Richardson et al. 2006; Nabeel-Shah et al. 2014). NASP contains a tetratrico peptide repeat (TPR), which typically mediates protein–protein interactions. Known functions of TRP proteins include protein translocation, transcriptional regulation, pre-mRNA splicing and RNA-binding protein translation repression (D'Andrea and Regan 2003). In several crustacean species (P. monodon, Caligus rogercresseyi, Lepeophtheirus salmonis), NASP expression was found to be female-biased (Karoonuthaisiri et al. 2009; Farlora et al. 2014; Poley et al. 2016), suggesting a role of NASP in female crustacean development. In the nematode Caenorhabditis elegans, NASP forms part of a protein complex which represses male-specific gene transcription in hermaphrodites to promote female development, and may therefore be fundamentally involved in sex differentiation (Grote and Conradt 2006). If NASP holds a similar role in krill, it may be expressed in a female-specific way from the onset of molecular sexual differentiation.

PHL is a homolog of RAF, which is an important component of the Ras/Raf/MAPK pathway. Through a cascade of kinase activations this pathway can ultimately regulate transcription factors (Peyssonnaux and Eychène 2001). The MAPK pathway can be activated by a wide variety of cellular functions, and in vertebrates each kinase exists in several copies, which can additionally express multiple protein isoforms through alternative splicing (Barnier et al. 1995). This results in a very diverse and complex signalling cascade that may allow for tissue and developmental specific modes of gene regulation (Barnier et al. 1995; Peyssonnaux and Eychène 2001). Among many other functions, the MAPK pathway is involved in regulating male sex determination genes in mammals (Bogani et al. 2009; Windley and Wilhelm 2015) and can also play a role in oocyte maturation (Zhang and Liu 2002; Fan and Sun 2004). In drosophila, PHL is involved in the assembly of the eggshell and the specification of the terminal region of the embryo (Jimenez et al. 2002). In the giant tiger shrimp P. monodon, PHL is expressed in the ovaries of juvenile and mature females, but not in testes (Klinbunga et al. 2009), and also appears to be up-regulated in ovaries of the Pacific white shrimp Litopenaeus vannamei (Peng et al. 2015) and the swimming crab Portunus trituberculatus (Yu et al. 2015). Together with our findings of female-specific expression of PHL in krill, this suggests a female-specific role of this gene, such as a possible involvement in germ cell and oocyte development.

All three female-specific genes were found to be highly expressed in females to similar levels from early juvenile stages onwards. While this follows expectations for NASP, which may be required in female sexual differentiation from very early stages onwards, the expression of VRG has previously been found to correlate with advancing sexual maturity (Roth and Khalaila 2012; Bai et al. 2016). Though less is known about PHL expression in crustaceans, it may also increase expression with sexual maturation from juvenile to mature females (Klinbunga et al. 2009). The absence of a correlation of female gene expression with sexual maturity in our study may indicate that ovarian development in juvenile and sub-adult krill was already relatively advanced despite absent or incipient external morphological sexual differentiation. Alternatively, high expression of all three female-specific genes may be required from early sexual differentiation onwards in female krill.

Some care must be taken in interpreting results from regressed krill individuals used in this study. Regression has been found to be particularly strong when food availability was scarce over winter (Ikeda and Dixon 1982), although Thomas and Ikeda (1987) observed similar regression and re-maturation of female sexual organs independent of feeding conditions. In our study, light conditions in the krill aquarium were equivalent to winter conditions found at 66°S, however, feeding conditions were identical to summer conditions (Kawaguchi et al. 2010). This may explain why we only observed weak sexual regression, with many male and female individuals displaying a mature petasma and thelycum, respectively, although no sexually reproducing individuals were observed (data not shown). Thus, either the continuous expression of all three female-specific genes is required throughout winter regression in female krill, or our results may indicate that sexual regression under well-fed aquarium conditions is limited on a molecular level. Further experiments in the krill aquarium combining winter light conditions with reduced food availability could potentially disentangle the two effects. Alternatively, samples collected in the wild over the winter months would give valuable insights into sexual regression on a molecular level. However, winter ice conditions and limited scientific explorations during the winter months make it difficult to obtain such samples, although the krill fishery might be an alternative source for krill sampling in winter (Kawaguchi and Nicol 2007).

None of the three furcilial larval stages tested (furcilia III, IV and VI) expressed the female-specific genes in a total of thirty larvae samples. This indicates that sexual differentiation is initiated after the last larval stage (after furcilia VI, i.e. at some point during early juvenile development). There is the potential that other genes display sexually divergent expression slightly earlier than what we observed. However, the three female-specific genes that we investigated may all be fundamentally involved in female development, thus we suggest that sexual differentiation has probably not been initiated before these three genes are expressed. These findings are in agreement with observations of Cuzin-Roudy (1987), who did not observe sexually differentiated gonads in any furcilia stages nor in juvenile krill smaller than 23 mm that were examined, and contradict findings of Bargmann (1945), who appeared to have determined the sex of krill as small as 10 mm. At what point during early juvenile development sexual differentiation occurs exactly, and whether the gene expression initiation of all three female-specific genes is synchronised, requires further investigations. Ideally, this would involve an aquarium cohort of krill of known age that is regularly sampled after the last larval stage.

We did not identify any genes strongly expressed in whole-male samples but not in females. This may be related to our approach: we used transcriptomic data from ovaries and testes as a starting point to discover sex-biased genes. As we used pooled gonad material of only a few individuals and had therefore no biological replica, genuine male-specific expression patterns may have been masked. Additionally, many genes essential for male sex differentiation, such as androgenic gland hormone, may be exclusively expressed in the androgenic gland (Sagi et al. 1997), and these genes could have been missed in this study. While we investigated the expression of “insulin-like androgenic gland hormone binding protein” (Li et al. 2015) and could not detect any sex bias (data not shown), the androgenic gland hormone itself is not annotated in either of the online krill transcriptome databases, potentially due to relatively large sequence differences to other crustacean species. Another explanation for our difficulties in identifying male-specific genes may lie in the properties of the krill genome: many genes occur in multiple copies (Deagle et al. 2015) and gene duplication can give rise to novel gene function (Zhang 2003). If a male-biased gene exists in multiple copies in the genome, but only certain variants hold a function in male sexual differentiation, detection of sex-biased expression of that particular variant could be challenging. Alternatively, sex-specific gene expression may be more prevalent in females than in males. Further investigations using different approaches are thus needed to identify genes expressed exclusively in males.

ANK and AQP12A, the two genes that were found to show moderate male expression bias in adult krill, are probably not directly involved in male sex differentiation, as neither juvenile nor sub-adult krill showed sex-specific expression. Consequently, neither of these genes is suitable as a molecular determinant of sex in krill lacking external sexual characteristics. Additionally, their function in male differentiation or reproduction is difficult to interpret, as both genes can be involved in a variety of cellular processes. Ankyrin repeats are among the most common existing protein motifs in nature, mediate protein–protein interactions and have a wide variety of functions (Li et al. 2006). An Ankyrin gene has been found to potentially hold a role in sex determination in the Chinese mitten crab Eriocheir sinensis (Cui et al. 2015), although the expression was found to be slightly female-biased, in contrast to our results. It is currently not clear what male-specific role this gene could play in krill. Aquaporins are channel proteins that transport water or other small solubles through cell membranes (Perez Di Giorgio et al. 2014) and can hold a wide variety of functions during reproduction (Huang et al. 2006). Male-biased expression of AQP12A has been found in Salmon lice L. salmonis (Poley et al. 2016), and another Aquaporin gene has also been found to be strongly male-biased in the Salmon louse C. rogercresseyi (Farlora et al. 2014), suggesting a similar male-biased function of AQP12A in our study.

In summary, we identified three genes that are expressed exclusively in females from juvenile stages onwards. This allows for an unambiguous molecular test to determine the sex of krill samples that are lacking sex-specific external morphological features, with the exception of larval samples. Based on our results we suggest that molecular sex differentiation is not initiated until after the last larval stage at some point during early juvenile development. These findings will now allow for the first time to determine sex ratios in krill swarms lacking external sexual morphological features (juveniles, regressed krill), which will contribute to our understanding of krill population dynamics, genetic diversity and reproductive behaviour.