Ecological genomics in the Northern krill uncovers loci for local adaptation across ocean basins

Domain: earth_science
Published: 2024-08-01
Editorial outcome: accept

Abstract
Krill are vital as food for many marine animals but also impacted by global warming. To learn how they and other zooplankton may adapt to a warmer world we studied local adaptation in the widespread Northern krill (Meganyctiphanes norvegica). We assemble and characterize its large genome and compare genome-scale variation among 74 specimens from the colder Atlantic Ocean and warmer Mediterranean Sea. The 19 Gb genome likely evolved through proliferation of retrotransposons, now targeted for inactivation by extensive DNA methylation, and contains many duplicated genes associated with molting and vision. Analysis of 760 million SNPs indicates extensive homogenizing gene-flow among populations. Nevertheless, we detect signatures of adaptive divergence across hundreds of genes, implicated in photoreception, circadian regulation, reproduction and thermal tolerance, indicating polygenic adaptation to light and temperature. The top gene candidate for ecological adaptation was nrf-6, a lipid transporter with a Mediterranean variant that may contribute to early spring reproduction. Such variation could become increasingly important for fitness in Atlantic stocks. Our study underscores the widespread but uneven distribution of adaptive variation, necessitating characterization of genetic variation among natural zooplankton populations to understand their adaptive potential, predict risks and support ocean conservation in the face of climate change.

Introduction
Climate change is affecting all life on Earth and forcing species to move or adapt1. Ocean zooplankton are crucial to maintaining food webs and fisheries but face many challenges including increased temperatures and acidification2–4. Many planktonic species are shifting toward higher latitudes4,5, and continued warming is expected to impact marine communities and ecosystem services6. The long-term responses to these changes are unclear, but evolutionary adaptation may be important to sustain populations, particularly when physiological and geographical limits have otherwise been reached7. This warrants the need to better understand adaptation in key zooplankton that strongly influence marine ecosystems. Krill (Euphausiacea; 86 spp.), or euphausiids, are macrozooplankton crustaceans inhabiting all world oceans. Some species include trillions of individuals and are among the most abundant animals on Earth8–10. As grazers of smaller plankton and food for fish and mammals, krill are critical links between primary production and higher trophic levels3. However, polar krill of both hemispheres have recently declined11,12, while boreal species such as the Northern krill spread into new areas13, impacting native biodiversity14,15.

The Northern krill is the largest and most abundant North Atlantic krill species, possibly structured into 3–4 basin-scale gene pools9,16. While many krill species are stenothermal and have narrow latitudinal ranges, the Northern krill has unusually broad thermal tolerance and range17,18. It occurs across a 2–15 °C temperature gradient (Fig. 1a) and breeds within 5–15 °C9, much wider thermal envelopes than for example the Antarctic krill Euphausia superba, which is constrained within −2.0 °C to +4.0 °C and reproductively challenged already at +1.5 °C19. Northern krill from different climates vary in metabolism, nutrition, maturation and timing of reproduction that track local seasonal cycles17,20, ranging from spawning in late winter–early spring in the Mediterranean Sea to summer in the Atlantic Ocean. These phenotypic variations could have genetic bases, making M. norvegica an attractive model for environmental adaptation.

Zooplankton generally have large populations with considerable genetic variation, suggesting they have high potential to adapt to changing environments through natural selection21. Many also exhibit extensive larval dispersal and gene flow, which could either interfere with local adaptation or help introduce adaptive variants21. Insights into the genetic basis of adaptation are still scarce for zooplankton22. In particular, we lack insight into adaptation in euphausiids, due to their large and repetitive genomes ranging between 11 and 48 Gb (4–16× the human genome)23, which until recently have hindered genetic analysis. The publication of the Antarctic krill genome marked a major step forward but revealed extensive genetic homogeneity and limited adaptive variation in this circumpolar and panmictic species24. A comparative analysis of expressed genes in 20 krill species identified mechanisms of adaptation but could only reliably interrogate about 1/3 of all genes and focused on cold-adaptation in polar species25.

Here we used long-read sequencing to assemble the huge genome of the widespread Northern krill. We then performed genome-scale population resequencing with the aim to characterize genetic variation, demographic history, and genetic adaptations to both cold and warm environments at the genome-scale. Insight into adaptive variation can help monitoring processes such as migration and adaptive change, identify hot and cold spots of important functional diversity and support forecasting resilience and risk under continued climate change26.

Results
A highly repeated genome
The Northern krill is diploid and has 19 large metacentric chromosome-pairs that are homomorphic between sexes27. We assembled and scaffolded the nuclear and mitochondrial genomes with Nanopore long-reads, Chromium linked-reads and RNA data from mostly a single specimen (Supplementary Table 1, Supplementary Data 1, 2, Supplementary Figs. 1–3). We tracked contiguity, completeness and accuracy by mapping coding transcripts (Supplementary Figs. 4, 5). The finished genome assembly spanned 19.2 Gb (n = 216,568 scaffolds/contigs) with ~90% of BUSCO genes present and ~0.5 base-level errors/kb (Fig. 1b, Supplementary Table 2). The genome is GC-poor (%GC = 29.9) and repeat-rich. Low-complexity sequences and simple repeats spanned 15% of the genome (Supplementary Data 3). Using a custom transposable element (TE) library, we found that 74% of the genome is repetitive (Fig. 1c; Supplementary Data 3). While we are unable to classify all TEs, retrotransposons (LINEs+LTRs) outweighed DNA-transposons 4:1, similar to the American lobster or black tiger shrimp28,29. The repeatome in the Northern krill exhibits a unimodal pattern of TE divergence (Fig. 1c), whereas the Antarctic krill repeatome is characterized by DNA-transposons and shows signs of several bursts of transposable element proliferation24. These observations hint at divergent genome composition and lineage-specific evolution of the gigantic genomes of euphausiids.

Expansion of cuticular and opsin gene families
We used RNA and comparative data to annotate 25,301 protein-coding genes along with 2283 TEs (mostly expressed retrotransposons), and detected another 14,643 potential yet unannotated genes or TEs (Supplementary Tables 3, 4, Supplementary Data 4, Supplementary Fig. 6). Gene bodies (introns + exons) span 50,276 bp on average and occupy 10% of the genome, while coding sequences cover only 0.22% (Fig. 1b). Orthologous genes between the Northern krill and crustaceans with smaller genomes are 3–8× longer in krill, but have similar amounts of coding sequence (Fig. 1d, Supplementary Fig. 7, Supplementary Table 5, Supplementary Data 5). Compared to the Antarctic krill24, genes are ~2.5× longer in the Northern krill, suggesting proliferation of retrotransposable elements has produced long and repeated introns (Supplementary Fig. 6). We estimated high synonymous divergence (dS = 0.46) between the two species (Supplementary Data 5). Using a decapod molecular clock30, this divergence suggests they split from a common ancestor ~130 MYA, underscoring separate evolution over long time-scales.

We built a crustacean species tree and analyzed gene family evolution. We found 104 rapidly expanding gene families in the Northern krill (p < 0.05; Fig. 1e; Supplementary Fig. 8a, b; Supplementary Data 6), including those related to chitin, cuticular metabolism, regulation of the molting cycle, which are important processes for growth and reproduction in crustaceans. This is notable as renewal of the exoskeleton is unusually frequent and plastic in euphausiids9,31, and similar expansions were independently detected in the Antarctic krill genome24. Moreover, we detected expansions of the opsin gene repertoire, which encodes the light-sensitive receptors in ommatidia. Fourteen opsins have previously been identified from RNA in E. superba32, which are thought to enable vision under the divergent light conditions experienced throughout its life cycle and vertical migrations33, while 16 opsins have recently been inferred from M. norvegica RNA34. We queried our M. norvegica gene-set and the E. superba RNAs32 against a curated crustacean opsin dataset. We detected 19 genes in the former species and 15 putative genes in the latter (Supplementary Figs. 9, 10), including new visual middle wavelength-sensitive (MWS) opsins and non-visual arthropsins. All E. superba opsins have homologs in the M. norvegica genome and all but one previously identified M. norvegica transcript can be anchored unambiguously 1:1 to our gene models (Supplementary Fig. 10). Our findings expand the known opsins in both species and suggest that opsin and molting-gene duplications could be common to all euphausiids.

Ancient whole-genome duplication (WGD) may underlie both genome and gene family expansions but is not commonly reported for crustaceans. To test for WGD, we interrogated divergences between gene-paralogs and searched for Hox-gene duplications, but found no supporting evidence (Supplementary Note 1). Instead, the huge Northern krill genome has likely evolved through TE proliferation and numerous small-scale duplications.

An active DNA methylation system
Epigenetic regulation of the genome may contribute to genome evolution, function and plastic responses to environmental change35. CpG-methylation can silence the expression of harmful transposable elements (TEs), paradoxically allowing them to persist and contribute to genome expansion36. DNA methylation (DNAm) of both TEs and protein-coding genes is ancestral in Arthropoda, but has frequently been lost37. Using DNA from muscle, we characterize the DNAm toolkit and genomic patterns in a euphausiid. The genome had low CpG O/E content (CpG = 0.53), indicating DNAm. We found all genes encoding the canonical methyltransferases and genes responsible for repair or demethylation (dnmt1–3; alkB2; tet2; Supplementary Figs. 12–16; Supplementary Data 7), hallmarks of functional DNAm37. We scanned the Nanopore-reads for signals of CpG-methylation at 75 million CpG-sites. Overall, DNAm rates were higher in genes than intergenic regions (75% vs. 62% of reads being methylated; Fig. 1f; Supplementary Fig. 17) and positively associated with splice isoform variation (n = 2.5 ± 0.05 isoforms/gene with >95% methylation rates vs. 1.8 ± 0.03 isoforms/gene with <5% methylation rates; Supplementary Fig. 17). In contrast, we observed a lack of methylation in the mitochondrial chromosome (4%). DNAm rates are higher across repeats vs. nonrepeated DNA (69% vs. 50%), and appear to target young retrotransposons with similar LTRs, that may have been recently active (Supplementary Fig. 17f). Gene-body methylation is similar to observations in marbled crayfish38, while repeat-oriented methylation is reminiscent of the myriapod Strigamia maritima37. The krill methylome spans both genes and repeats, suggesting dual roles in gene regulation and silencing TEs.

Genome-scale variation is shaped by linked selection and pervasive gene flow
To uncover patterns of genetic variation, we collected 74 Northern krill specimens from eight geographical regions separated by up to 5800 km, covering a wide range of environmental conditions (Fig. 1a, Supplementary Note 2). We sequenced whole genomes to ~3×/specimen (20–30×/population), mapped reads and called 760 million quality-filtered SNPs across 8.4 billion accessible bases (determined with similar filters; Supplementary Fig. 18). We estimated genome-wide nucleotide diversity (π; the average pairwise difference between sequences) and the population mutation rate (θw) to be 1.31% and 1.62% per-base, respectively, similar to those inferred from RNA25 and low compared to marine broadcast spawners such as oysters and sea-squirts39. Levels of variation were similar among populations (Supplementary Fig. 19A). Assuming mutation-drift equilibrium and a mutation rate from snapping shrimp30, we inferred the long-term effective population size (NE) to be 1.53 million, far below the expected census population size (NC) of trillions9. However, Tajima’s D was negative (−0.53), indicating an excess of low-frequency variants compared to expectation under equilibrium, consistent with population expansion shaping genetic variation. We applied the Sequentially Markovian Coalescent to model demographic history from heterozygous sites in the reference specimen40,41 (~5000 longest scaffolds), suggesting populations expanded halfway through the last glacial period42 (Fig. 2b). The rate of decay of linkage disequilibrium between SNPs in turn suggested NE may recently have reached 4–5 million (Supplementary Fig. 19B). We observed ~30% reduction in variation over genes (π = 1.15%), and more so at coding and nonsynonymous sites (Fig. 2a). This effect extended up to 50–100 kb around genes (Fig. 2a), suggesting widespread impact of linked selection on genetic variation.

Population structure was limited but recapitulated geography (the fixation index FST was ≈ 0.06 on average; Fig. 2c) and previously detected mitochondrial gene pools9. FST (relative divergence) increased marginally with geographic distance (Supplementary Figs. 19, 20, Supplementary Data 8). The Mediterranean Sea sample was the most divergent, but the averaged DXY (absolute divergence) was only 1.04× higher compared to distances among Atlantic populations (1.71% vs. 1.64%). Net synonymous divergence (Da) between Mediterranean and Atlantic stocks was very low (<0.1%), reaffirming they are not reproductively isolated43. Moreover, most non-singleton variants were polymorphic in most populations (Supplementary Fig. 21), further indicating extensive gene flow.

Signatures of ancient adaptive divergence across hundreds of genes
To reveal genomic signatures of adaptation, we partitioned the dataset into two major contrasts: i) Atlantic vs. Mediterranean samples (at/me); and ii) North-Eastern vs. South-Western North Atlantic samples (ea/we; Fig. 3a). For each contrast, we computed pairwise divergence in allele frequencies. While most variation segregated at low differences (FST = 0.056; FST(at/me) = 0.017; Fig. 3b; Supplementary Fig. 20), we also detected approximately 8× as many highly divergent variants compared to expectations from simulations of neutral drift (Supplementary Fig. 22). Divergent regions (FST > 0.4) spanned ~1% or less of the genome and were about 2× enriched for gene sequences and extended haplotypes, compared to undifferentiated regions (FST < 0.2; Fig. 3b), consistent with gene-centered signatures of selective sweeps44. We compared FST between genes and similarly sized 50 kb flanking regions 50–100 kb away from genes, which may more often evolve neutrally. At high levels of divergence (the top 0.1% most divergent flanking regions), genes outweighed flanking regions by 7× (at/me) or 2× (ea/we), consistent with natural selection driving adaptive divergence across many genes (Fig. 3c; Supplementary Data 8; Supplementary Fig. 23).

We analyzed the geographic distribution of putatively adaptive variation by defining gene-level haplotypes (having at least four diagnostic SNPs with FST > 0.5). At many divergent genes (exon-wide FST > 0.4), both haplotypes were often present in both groups (n = 319/660; n = 34/34), indicating widespread standing variation (Fig. 3d). Southern or Scandinavian populations were more polymorphic than Barents Sea and Svalbard populations (Fig. 3d), which could reflect genetic drift or ongoing selection at the margin of the Arctic species range45. We found that Mediterranean haplotypes were comparably common in the Norwegian sample (mean frequency of 13%, Fig. 3d-iii), while South-Western haplotypes were frequent in the Swedish sample (mean frequency of 17%, Fig. 3d-vi), suggesting these Scandinavian stocks in particular could contain genetic material that is otherwise rare at high latitudes. We estimated the ages of minor alleles on the divergent haplotypes to learn for how long haplotypes may have been segregating in the species. We first estimated a genome-average recombination rate (r = 0.32 cM/Mb, n = 652 scaffolds) and then applied the Genealogical Estimation of Variant Age (GEVA)46. This coalescent method infers the time to the most recent common ancestor using mutation and recombination rates, without requiring a priori assumptions about demographic history. Most variation originated over 1 MYA, predating multiple glacial cycles, and adaptive variation segregating between Atlantic and Mediterranean populations may predate that segregating in the Atlantic Ocean (Fig. 3e).

Candidate genes for ecological adaptation associated with ecophysiological functions
For each contrast, we ranked all genes by exon-wide FST. The top gene between Atlantic and Mediterranean krill is nose resistant to fluoxetine protein 6 (nrf-6) (Supplementary Data 8), encoding a membrane protein facilitating lipid transport. It is located within a high-divergence region and the Mediterranean samples are fixed for an nrf-6 haplotype with low frequency among Atlantic samples (6%). Exon-wide FST is > 0.8 and high XP-nSL scores indicate loss of variation consistent with a selective sweep among Mediterranean samples (Fig. 4a). We also detect accelerated protein evolution on the Mediterranean haplotype (dN/dS: 0.54 vs. 0.31; Supplementary Fig. 24), consistent with positive selection. We predicted the protein topology and detected 1.7× enrichment of missense variants in its extracellular part compared to synonymous variants (9/16 vs. 5/15; Fig. 4a; Supplementary Fig. 24; Supplementary Data 9), which may alter its function. This is noteworthy as nrf-6 is important for yolk transport into eggs in worm47 and ovary development and oogenesis in fly48, and is overexpressed in the ovaries of sexually precocious crabs49. In the Northern krill, Mediterranean stocks depend on a short spring phytoplankton bloom to accumulate lipid stores and trigger vitellogenesis and reproduction in early spring9,50. The Mediterranean nrf-6 variant could contribute to advanced reproductive timing, making it a strong candidate for local adaptation.

We next queried the ranked lists for enriched gene ontologies using fly homologs. Atlantic-Mediterranean divergence was enriched for genes involved in oogenesis, muscle function, phototransduction and eye development, regulation of circadian rhythm and heterochronic development (Supplementary Fig. 25; Supplementary Data 10). These include homologs for ninaB that synthesize visual pigment51, S that regulate sleep/wake cycles and the transcription factor Kr-h1 that acts in the juvenile hormone signaling pathway to govern vitellogenesis and reproduction in arthropods52 (Fig. 4b). Overall, 30 genes in the 0.1%-percentile (FST = 0.55–0.72; n = 264) were associated with vision-related ontologies, 1.5× more than expected by chance (P = 0.019), and 112 genes in the 1%-percentile (FST = 0.36–0.72; n = 1024), including four MWS opsins (Supplementary Data 10). Photoreception candidates often involved paralogs (Supplementary Fig. 26) and belonged to rapidly expanded gene families, which were otherwise underrepresented in the 1%-percentile (2.98× vs. 0.6×; Supplementary Data 10).

Between the two Atlantic basins, the most divergent gene encodes a larval cuticular protein (homolog of Cpr49Ac; Fig. 4c), while this contrast was enriched for signaling ontologies (Supplementary Figs. 25). Among the top 20 genes we identified three heat-shock chaperones/co-chaperones, including CCT8 (Fig. 4c), that fold proteins and promote proteostasis under thermal stress. Heat-shock proteins contribute to protective thermal tolerance in krill and many other species53. Chaperonins (CCT-n) influence cold shock response in crustaceans and other eukaryotes54. Multiple CCT genes were linked to cold-adaptation in krill25 and CCT8 evolution was implicated in freeze-tolerance in Amur sleeper fish55. North-Eastern krill have a sweep-like signature of extended CCT8 haplotypes, suggesting selection for cold tolerance in these stocks.

Methods
Biological materials
Adult specimens of the Northern krill (Meganyctiphanes norvegica) were collected from eight geographical regions across the North Atlantic Ocean and the Mediterranean Sea by different collaborators and co-authors (Supplementary Table 1; Fig. 1a; Supplementary Fig. 19C). Environmental characteristics and sampling at each location is provided in Supplementary Note 2, while Sea Surface Temperatures in Fig. 1a are from Climate Reanalyzer using NOAA OI SST V2 High Resolution Dataset data provided by the NOAA PSL, Boulder, Colorado, USA75,76. Sample details including sampling coordinates and sequencing yields are given in Supplementary Data 1.

Processing of samples. Live specimens were either preserved in 95-99% EtOH or in RNAlater™ (Invitrogen™) to facilitate DNA and RNA recovery. Some specimens were split behind the carapace so that the cephalothorax and abdomen could be stored separately in either buffer. For most samples, the sex was determined by inspecting petasma or thelycum for male or female characteristics using a stereo microscope. This was either done at the time of collection or after preservation in EtOH. Altogether, 36 samples were determined to be males, 32 to be females. Seven samples were not sexed because the relevant tissue had been removed or damaged in processing.

A single female specimen (Sample IDs “K20” or “swe_1”; ~3 cm long) collected in the Gullmarsfjord in Sweden was used to produce all primary long-read and linked-read DNA data for the main de novo genome assembly, as well as most novel short-read and long-read RNA/cDNA resources used for scaffolding and annotation of the genome sequence. In order to reduce the risk of contamination from food-particles, this specimen and others sampled at the same location were kept in aquaria with filtered running deep water at ambient temperatures (10 °C) at the Kristineberg Marine Research Station without access to food for 24–48 h before they were preserved. Sequences from other specimens were used to assist in specific tasks for the genome assembly (Supplementary Table 1; Supplementary Fig. 2), as well as to produce a comprehensive population dataset to map genetic variation.

Overview of library preparation & sequencing
For the genome assembly and annotation, we produced DNA and RNA/cDNA data using Oxford Nanopore Technologies (ONT) MinION (hereafter: MinION) and PromethION long-read sequencing (hereafter: PromethION). Linked-read DNA data was produced using 10x Genomics Chromium sequencing (hereafter: 10x) (Cat#120257/58/61/62). RNA was also sequenced using Illumina RNA-seq. Low-coverage whole-genome population-scale resequencing (WGS) was carried out using Nextera DNA Flex libraries (Illumina Cat# 20018705) or Illumina DNA PCR-Free libraries (Illumina Cat# 20041795). All DNA/RNA extractions and the MinION sequencing were performed at the Dept. of Medical Biochemistry and Microbiology at Uppsala University. All other library preparation and sequencing was carried out by the National Genomics Infrastructure (NGI) at Science for Life Laboratory (SciLifeLab), Sweden. Extractions were quantified for purity, yield and fragment lengths using the Thermo Scientific NanoDrop and Qubit systems, standard 1% agarose gels, and Agilent TapeStation and Femto Pulse systems, respectively. Summary reports with basic sequence statistics were produced using PycoQC for long-reads77 and MultiQC for short-reads78 with the high-throughput FastQC backend (Simon Andrews, Babraham Institute in Cambridge), respectively, according to standard procedure. See Supplementary Data 2 for additional kit catalog IDs.

High-molecular weight DNA extraction for long-read and linked-read libraries. We produced DNA data using PromethION and 10x sequencing of a single female specimen for the genome assembly. We dissected five abdominal segments for all ethanol-preserved tail muscle tissue and extracted each segment individually (K20-1 through K20-5). Each piece was snap frozen in liquid nitrogen and grinded using a pestle and mortar. High-molecular weight (HMW) DNA was extracted using the Genomic-tip 100/G Midi kit (Qiagen) according to the manufacturer’s protocol for tissues. Briefly, the pieces of tissue were lysed in presence of Protease (Qiagen) and RNAse A (Qiagen) to degrade proteins and RNA, respectively. The DNA was bound to the Genomic-tip resin, purified from RNA, protein and other contaminants, eluted and desalted using isopropanol precipitation. The DNA extractions were finally eluted in 70 μL of TE buffer (Qiagen AE buffer). The extraction yielded approximately 30 μg DNA (Qubit), with a peak of DNA fragment lengths at about 130 to 180 kbp, and with similar yields and qualities across the five abdominal segments (Supplementary Fig. 1A, B).

PromethION long-read sequencing. The ONT SQK-LSK109 Ligation Sequencing Kit was used by NGI to prepare DNA libraries using all five extractions (K201–5; 90% of DNA was consumed), which were then sequenced across ten PromethION FLO-PRO002 R9.4 flow cells. The DNA was sheared to 20 or 30 kb fragments with a Diagenode Diagnostics Megaruptor before preparing most remaining libraries to find an optimal balance and tradeoff between read lengths and yields (Supplementary Fig. 1C). One HMW library was made with the Megaruptor set to shear to 75 kb, followed by size selection for fragments >15 kb using the Sage Science Blue Pippin system. Library preparation for nanopore sequencing was done using the LSK-109 kit, with the following changes to the original protocol: i) End-prep/FFPE incubation times at 20 and 65 degrees were prolonged up to 10 or 15 min, and ii) adapter ligation time was prolonged up to 55 min. Sequencing was performed on a PromethION beta machine, using R.9.4 flow cells. In order to increase the overall yield, we used the flow cells. We first washed them using standard wash kit buffers and then performed nuclease flushes to digest old libraries, before applying new DNA. Two flushed flow cells were run with cDNA libraries instead of DNA (see below).

Sequencing was performed using ONT MinKNOW v3.1.18. Initial base calls were carried out using MinKNOW/Albacore but were updated with calls using ONT Guppy v2.3.7 and the more accurate “flip-flop” algorithm (the dna_r9.4.1_450bps_flipflop_prom model)79. The updated base calls produced 593.8 Gb of DNA data (13% produced on flushed cells), representing 13% more data than the original base calls. The Guppy-recall took 10 days of GPU-compute time across four NVIDIA 1080 GPUs. We sequenced 84.4 M reads with overall N50 = 12 kb. About 57% of sequences had mean base qualities (Q-scores) above 10 using flip-flop (Supplementary Data 2), compared to 16% in the original call. Assuming a genome size of 18.6 Gbp23, the PromethION data corresponded to 32× coverage across the genome, with ~20× coverage of reads 10 kb or longer (Supplementary Data 2; Supplementary Fig. 1D).

10X linked-read sequencing. Four separate 10x GEM libraries were made from the K20-5 HMW DNA extraction by NGI following the manufacturer’s recommended protocol and had mean library fragment lengths of 530–540 bp. Paired-end reads (2×150 bp) were sequenced on a single Illumina NovaSeq 6000 S4 lane, generating 897 Gb of data across 2.97 B read-pairs (~46× coverage; 12.6% marked as duplicates; %GC = 34.7). The standard LongRanger (v2.2.2; 10x Genomics) “basic” command was used to extract the 16 bp barcodes from the forward/R1 reads for each library individually. About 3.5 M of the 4.8 M barcodes (4M-with-alts-february-2016.txt; 10x Genomics) were found to be associated with at least one read in each library and ~96% of all sequences were barcoded. We rewrote the “BX:Z” style barcodes of each of the four libraries using a unique suffix (“_1” through “_4”) in order to keep them separate in downstream analyses and to de-interleave the FASTQ files.

RNA extraction from multiple tissues for long-read and short-read RNA libraries. We produced RNA data using PromethION/minION cDNA-sequencing and Illumina short-read RNA-seq using the K20 specimen for gene and genome annotation. Six libraries were prepared using RNA from small pieces of tissue preserved in RNAlater: i) head/brain; ii) dorsal/caudal muscle tissue on cephalothorax; iii) thoracic exopodite legs (distal segments of the periopods); iv) ventral soft tissue of cephalothorax including gills; v) hepatopancreas (the crustacean digestive gland); and vi) soft-tissue surrounding it. RNA was extracted using the RNeasy Lipid Tissue Mini Kit (Qiagen) according to the standard protocol, using ceramic beads and a BioSpec Products Mini-Beadbeater for tissue homogenization, and eluted in 30 μL RNAse-free water. This yielded RNA extractions with concentrations between 50 and 240 ng/μL and with RIN-values between 9.1 and 10. We used both raw RNA-seq reads and Trinity assembled transcripts80 in downstream analyses.

Short-read RNA sequencing. NGI produced six libraries from the six K20 RNA extractions using Illumina TruSeq Stranded mRNA Library Prep kit with polyA selection (Illumina Cat# 20020594/5, manufacturers’ protocol #1000000040498), which were sequenced as 2×150 bp paired-end reads on one Illumina NovaSeq 6000 S4 lane. This generated 459 M read-pairs that were demultiplexed by NGI (Supplementary Data 2).

Long-read cDNA sequencing. One barcoded and amplified Nanopore cDNA library was produced from each of the six K20 RNA extractions, using the PCR-based ONT Strand-switching protocol with poly-A selection and the LSK109 ligation kit to add barcodes and adapters. We performed cDNA sequencing on PromethION (ONT MinKNOW v3.1.18) using two nuclease flushed flow cells and one new FLO-MIN106D R9.4.1 MinION flow cell (3 M reads) and the standard loading kit and buffers for each instrument. The data was base called using ONT Guppy v4.4.1 and HAC (high accuracy) model (dna_r9.4.1_450bps_hac_prom.cfg). The two PromethION runs generated 8.99 M and 3.91 M reads (10.1+5.2 Gb of cDNA data), with median Q-scores of 10.5 and 10.7 and median lengths of 962 bp and 1060 bp, respectively. The MinION run produced 3.0 M reads (3.5 Gb of cDNA data) with median Q-score of 12.6 and median lengths of 915 bp.

Salt-isopropanol DNA extraction of multiple specimens for population-scale whole-genome resequencing. To produce the population dataset, we extracted DNA from 75 M. norvegica specimens using a salt-isopropanol precipitation protocol. We extracted DNA from abdominal muscle tissue from 1–2 muscle segments of krill samples. Samples were selected based on geographic origin and quality of preservation. We sequenced 6–10 specimens per population. We used a column-free salt-isopropanol precipitation protocol to extract DNA. Briefly, tissues were lysed overnight in ATL Lysis buffer (Qiagen) together with Proteinase K (Qiagen) overnight, followed by RNAse A (Qiagen) treatment, protein separation using Protein Precipitation Solution (Qiagen) and precipitation and washing of the DNA using isopropanol and ethanol, respectively. The washed DNA pellet was finally eluted in 70 μL MilliQ-water. The extractions produced about 7 μg of DNA on average per specimen (Qubit). The protocol is available as a separate document in this publication and at the online platform protocols.io: • https://www.protocols.io/workspaces/wallberg.lab.imbim.uu/resources/isopropanol-hmw-dna-extraction-method

Short-read DNA libraries and resequencing. Barcoded Nextera DNA Flex or Illumina DNA PCR-Free libraries (Illumina) were prepared from 75 DNA extractions (74 specimens for resequencing + the reference specimen) by NGI/SciLifeLab according to the manufacturer’s instructions. Both protocols use a bead-linked transposome to fragment DNA and incorporate adapters in a single step. About 350 ng of DNA was used as input per Nextera library and 300 ng per PCR-Free library. We aimed for ~3× depth of coverage per specimen and sequenced the libraries across six Illumina NovaSeq 6000 S4 lanes, generating 15.1 B 2×150 bp read-pairs with valid library barcodes (18% or 23% of reads marked as duplicates with Nextera or DNA PCR-Free, respectively; %GC = 32.9; 88% of reads with Q-scores above 30). One library (“sva_1” from the Svalbard area) failed to sequence (0.03× coverage; Supplementary Data 1) and was discarded from population genetic analyses (hence 74 samples were carried forward).

A preliminary mitochondrial assembly
A MinION test run of specimen K4 from the Gullmarsfjord produced 1.5 Gbp of long-read across 270,185 reads (N50 = 8229 bp; <0.1× across the genome) ahead of the other sequencing runs. We downloaded mitochondrial M. norvegica template accessions (16S=AY744910; COI=FJ581747; CYTB=AF149775; NADH=AF149775; tRNA_ser=AF149775; tRNA_leu_16S.AF149775) and used BLASTN 2.2.31+ to identify 235 putatively mitochondrial reads shorter than 16 kbp that aligned to any of the templates. These reads were assembled in Canu v1.781 with “genomeSize=17k” that produced a 24.8 kbp contig. We blasted the read pool against this contig and collected additional reads (n=520 in total) and performed a second round of assembly with Canu, generating a new 24.802 bp contig from 83 aligned reads (37× depth of coverage) that the assembler suggested was circular. Manual inspection of the sequence indicated that it was the full mitochondrial chromosome partly repeating itself and we circularized the sequence by hand, producing a 17,360 bp long sequence. We next mapped the long-reads to the mitochondrial sequence with minimap2 v2.17-r94182 (109× mean depth of coverage) and polished it with Nanopolish v0.11.383 using the “variants” and “vcf2fasta” subcommands. We took the high coverage as an indicator of high mitochondrial:nuclear ratio in the extracted muscle tissue. Finally, we mapped RNA-seq data from the published TI787-MN-F specimen to the polished sequence using Bowtie2 v2.3.4.184 and performed one round of polishing with Pilon v1.2285 with the “--fix indels,gaps,local” setting. The resulting mitochondrial sequence was 17,508 bp with %GC = 28. Automated annotation was performed with the online MITOS2 tool86. We used this sequence to screen out potential mitochondrial long-reads from the K20 reference datasets ahead of genome assembly and used a similar approach to assemble the mitochondrial chromosome from this specimen (see mitochondrial re-assembly below).

Pre-processing of long-read DNA data
We used Porechop v0.2.479 to trim the LSK109/NSK007 top and bottom adapters from the PromethION long-reads. Reads with internal adapters can be chimeras and either be split or removed. We first ran the program without splitting potentially chimeric reads:
porechop -t 20 --no_split --adapter_threshold 60 --end_threshold 70 --end_size 200 --in fastq.gz --out fastq.trimmed.gz
We then repeated the scan without “--no_split” to note which reads had middle adapters (n=50,368; 0.06% of all reads). We next mapped all reads against our preliminary mitochondrial sequence with minimap2 v2.17-r941 and found that 382 K reads mapped to this chromosome for a mean depth coverage of 33,662× (vs about 32× across the nuclear genome). We filtered out both the putative chimeric and mitochondrial long-reads for the reference PromethION FASTQ read-pool, to reduce the risk of downstream interference in the main genome assembly.

RNA processing and transcriptome assembly
Long-read cDNA processing. The barcoded ONT PromethION and minION cDNA reads were de-multiplexed with ONT guppy_barcoder v4.4.1 and then processed with ONT Pychopper v2.5.0 (https://github.com/nanoporetech/pychopper) to find, orient and trim full-length cDNA reads using the SSP/VNP primer sequences (cDNA_SSP_VNP.fas). The resulting 8.7 M reads were combined into a single archive and later used for gene annotation (see below). We also produced a non-redundant set of sequences to represent potential genes in the krill genome. Here we extracted long and high-quality cDNA reads (length > 500 bp and Q > 12) from both runs and clustered them into putative genes with isONclust287 using the “-x sahlin” mode. We then error-corrected the clustered cDNA reads with isONcorrect v0.0.688. Lastly, for each error-corrected cluster with four or more cDNA reads we generated a consensus sequence with vsearch v2.14.189:
vsearch --cluster_fast <clusters.fa> --threads 40 --id 0.95 --clusterout_id --clusterout_sort --sizeorder --sizeout --consout <clusters.fa.cons>
This produced 24,632 consensus sequences with N50 = 1,810 bp. The non-redundant transcript sequences were later used for RNA-based scaffolding (see below).

Short-read RNA-seq processing and transcriptome assembly. The Illumina RNA-seq data from the reference specimen and two published RNA-seq archives, TI787-MN-M (adult male; SRR3657321) and TI787-MN-F (adult female; SRR3657320), were trimmed for adapters and base qualities using Trim Galore! v0.6.1 (https://github.com/FelixKrueger/TrimGalore/) and Cutadapt v2.390:
trim_galore -j 8 --gzip --length 50 --trim-n --paired
All read pairs and 89–99% of bases were retained. The trimmed reads were later used in genome assembly scaffolding and annotation (see below). We used Trinity v2.5.180 to assemble a transcriptome from all 459 M trimmed read-pairs from the reference specimen, producing 573,869 with N50 = 900 bp. Using BUSCO v3.0.2b91 with the Arthropoda odb9 lineage set, we detected almost all expected genes (C: 97.2% [S: 42.4%, D: 54.8%], F: 1.6%, M: 1.2%, n: 1066). A subset of 500 bp or longer transcripts were used in scaffolding (see below). We processed all transcripts with Trinotate v3.1.192 and identified 60,677 possibly coding transcripts with best hits against metazoan templates from parsing the “Trinotate.xls” output and used these transcripts for gene annotation (see below). A subset of 16,509 non-redundant transcripts with best hits against arthropods were used to measure assembly quality scores (see below).

Mitochondrial re-assembly
We re-assembled the mitochondrial chromosome for the K20 reference specimen similarly as described above. PromethION reads identified during the screening process and that were also between 1–16 kb long were used to assemble the chromosome in Canu v2.0-development81 with “genomeSize=17k” that produced a 26.1 kb contig that the assembler suggested was circular. We next mapped the 369,946 mitochondrial long-reads to the sequence with minimap2 v2.17-r94182 and polished it with Racon using the “variants” and “vcf2fasta” subcommands. We circularized the sequence manually, producing a mitochondrial sequence that was 17,840 bp. We then mapped a subset of the 10x Chromium data (656,771 read pairs) to the chromosome with BWA and performed one round of polishing with Pilon v1.23 with the “--fix all” setting, followed by a round of polishing with K20 RNA-seq data. The final mitochondrial sequence was 17,944 bp with %GC = 28. Automated annotation was again carried out with the online MITOS2 platform as before and features were visualized with DNA Plotter93 (Supplementary Fig. 3).

Genome assembly
The hybrid genome assembly underwent multiple steps of polishing, scaffolding and refinements, in which we prioritized achieving high gene completeness and gene contiguity. We mainly used genetic information from the reference specimen, but also included data from other specimens for specific tasks (such as identifying redundant haplotigs). Nanopore long-reads were used for assembly, polishing (reducing base-level errors) and scaffolding, while 10x barcoded linked-reads and RNA data were used for polishing (in highly expressed regions) and scaffolding.

Genome assembly. Scaffolding was performed with top priority on recovering gene bodies using RNA from the reference specimen and additional RNA-libraries94, followed by genome-wide scaffolding using Nanopore synthetic mate-pairs with FAST-SG95 and 10x Chromium barcodes with Scaff10x. Some inadvertently removed contigs with genes were manually re-inserted into the genome. Screens against the SILVA ribosomal database96 were used to isolate contigs from bacterial contaminants. The assembly was carried out in 13 steps and we used BUSCO97 and genecovr (https://github.com/NBISweden/genecovr) throughout the pipeline to assess gene detection, completeness and sequence error.

Assessment of genome assembly quality and completeness. We used BUSCO v3.0.2b with the Arthropoda odb9 lineage set or BUSCO v5.0.097 with the Arthropoda odb10 lineage set, respectively, to track improvements in gene detection rates between genome assembly versions, and to assess the overall completeness and levels of duplication in the final genome. In addition, we assessed overall gene completeness by mapping a set of 16,509 trinity transcripts to the assemblies using gmap version 2020-06-3098. We parsed the psl output and summarized various statistics related to mapping quality, including number of insertions, number of mismatches, number of distinct contig hits, alignment length and alignment depth, to obtain “gene coverage” statistics for each transcript. Compilation of the statistical summaries and plots were implemented in a custom R package genecovr:
• https://github.com/NBISweden/genecovr.

Assembly v1: Production of a preliminary assembly with the wtdbg2 assembler. To assemble the krill genome, we used the wtdbg2 v2.5 assembler99 and 6–300 kbp long reads with Q-scores ≥ 10 selected with NanoFilt v2.6.0100 and that had been pre-screened for internal adapters or mitochondrial sequence, corresponding to ~16× depth of coverage across the genome. We used a short kmer setting (“-p 17”; as recommended for noisy data by the developers) and sparse kmer subsampling (“-S 2”) for sensitive read analysis and binning. Assembling the genome using these settings (“-A -p 17 -S 2 -s 0.05 -L 6000 -g 18g -t 168”) took 16 days on a 168 core machine with 4.7 TB RAM. Consensus contigs were exported with wtdbg-cns, resulting in genome assembly v1, which spanned 21.8 Gbp across 797,490 contigs with N50 = 41 kbp.

Assembly v2: Long-read polishing with Racon. The PromethION long-reads were mapped in parallel chunks to the v1 genome assembly using minimap2 v2.17-r941 set to use a split prefix and without secondary alignments (“-L -ax map-ont -y -t 16 --secondary=no --split-prefix”). Data was first saved as gzipped SAM files and then sorted according to coordinates into BAM files with samtools v1.10101, followed by merging into a single BAM file with samtools. The genome assembly was split into 248 100 Mbp FASTA chunks with a custom Perl script and the mapped reads were partitioned into corresponding SAM and FASTQ files using samtools. We used one round of Racon v1.4.11 with the average base quality threshold set to 8 (“-q 8”) to polish all 248 chunks in parallel. The polished contigs were concatenated into genome assembly v2.

Assembly v3: Long-read polishing with Medaka. We followed the guidelines from the medaka source documentation (https://github.com/nanoporetech/medaka#improving-parallelism) that suggests running the medaka components independently and in parallel for large genomes. Briefly, the PromethION long-reads were mapped in 250 parallel chunks to the v2 genome assembly using minimap2 v2.17-r941 with parameters “-x map-ont --MD --secondary=no -t 1”. This was followed by the consensus algorithm medaka consensus with parameters “--model r941_flip235 --batch_size 100 --threads 1” to generate models over each assembly chunk. Finally, consensus sequences were generated with medaka stitch and merged to obtain assembly v3.

Assembly v4: Short-read polishing with Pilon. The 10x Chromium paired-end linked-reads were contained in 16 paired FASTQ libraries that were mapped to the v3 genome assembly with BWA v0.7.17102,103. First, the genome reference sequence was indexed (“bwa index -a bwtsw”) and the data was then mapped with the “bwa mem” algorithm, piping the output to samtools which sorted the data on the fly before saving a BAM to disk. The 16 BAMs were merged with samtools. As in the Racon step above, both the reference genome and the BAM were split according to chunks of contigs spanning 100 Mbp. Polishing was performed in parallel with Pilon v1.23 set to polish putative base-level errors and specifying a diploid dataset (“--fix snps,indels --diploid”). The short-read polished contigs were concatenated into genome assembly v4.

Assembly v5: Purging haplotigs with Purge_haplotigs. Since our genome assembly was about 3 Gbp larger than the expected genome size (21.8 Gbp vs ~19 Gbp), we next sought to identify and remove potential haplotigs from the assembly using Purge_haplotigs v1.1.0104. The tool uses sequence similarity among contigs as well as a bimodal signature of read mapping depths to isolate and remove alternative alleles of the same locus from highly heterozygous diploid genome assemblies. It also identifies “junk” contigs with extreme mapping depths (low or high) that may be assembly artifacts or otherwise intractable in downstream analyses.

We first mapped the ~32× PromethION long-reads using Minimap2 to the v4 assembly, as in previous steps, and then used the “purge_haplotigs hist” command to generate a depth-of-coverage histogram across the genome. However, this relatively low-coverage dataset did not result in a clearly bimodal mapping profile with 1x and 0.5x mapping peaks on its own. We therefore also mapped the 10x Chromium linked-reads and a preliminary set of ~140× population resequencing data to the genome using BWA, in order to use all available data to specify depth thresholds to flag offending contigs. The short-read BAMs were marked for duplicates with Picard MarkDuplicates. We then exported per-base depth-of-coverage data running “samtools depth” across all positions in the genome assembly (“-a”) and generated a depth histogram with a discernable “1x” peak at ~150×, presumably corresponding to properly haplotype-fused contigs. The expected position of the 0.5x peak at half of this depth was less clear and possibly intermixed with genomic regions of low mappability. Given the depth profile, we then configured Purge_haplotigs to use 20× and 240× as lower and upper cut-offs, respectively, and 125× as the mid-point. We specified that contigs with 100% or less of their depths at diploid level of coverage to be flagged as suspected haplotigs (i.e. to comprehensively evaluate all contigs) and that those that had depths outside of the upper/lower cut-offs across >50% of their lengths to be flagged as junk:
purge_haplotigs cov -l 20 -m 125 -h 240 -s 100 -j 50 -i <coverage_file> -l 20 -m 125 -h 240 -o <coverage_stats_output>
Purge_haplotigs infers contig similarities from the proportion of sequence that can be aligned between two contigs using minimap2, and can be set to disregard alignments across repeat sequence. We therefore produced a preliminary catalog of repeat sequence intervals in BED format using Red105, masking 11.7 Gb (54%) of the genome. We used the default alignment score cut-off (≥ 70%) to mark contigs for reassignment as haplotigs:
purge_haplotigs purge -align_cov 70 -repeats <repeats_file>
The purge subcommand executes the script “purge.pl”. We modified this script to increase performance in the handling of our large, repeated and relatively fragmented genome assembly by: implementing faster parsing of tabular data; a string-based instead of array based data structure for hit summaries; and faster internal repeat-handling instead of executing the external tool bedtools from within the script. This code is available at: https://github.com/andreaswallberg/Ecological-Genomics-Northern-Krill.

Using these settings, we purged 168,316 contigs (N50 = 13,074 bp) spanning 1.95 Gbp as potential haplotigs and 49,996 contigs (N50 = 14,236 bp) spanning 565 Mbp as junk/artifact contigs, respectively, removing 12% of the original genome assembly. The remaining main assembly had 586,516 contigs (N50 = 46,208 bp) spanning 19.18 Gbp and was taken as genome assembly v5.

Assembly v6: Mixed-read polishing with HyPo. We aimed to polish the genome further after purging haplotigs and mapped both the PromethION long-reads and 10x Chromium linked-reads from the reference specimen back to the genome as above. We next split the genome assembly, BAMs and FASTQ reads (using samtools) into ten chunks of contigs spanning 2 Gbp of the genome each and ran the polisher HyPo v1.0.3106 separately on each chunk, which uses both short and long reads to polish contigs in a single run. The polished contigs were taken as genome assembly v6.

Assembly v7: Scaffolding the contigs with Trinity transcripts using L_RNA_Scaffolder. At this stage, quality assessments indicated that we had many broken gene models with exons distributed across more than one contig. Preliminary trials suggested that scaffolding contigs together based on information in RNA-mappings before long-read or linked-read DNA scaffolding recovered more complete gene models for this genome than vice-versa. We started by scaffolding contigs with non-redundant full-length cDNA transcripts (n=24,632; see section above) and Trinity-assembled transcripts, i.e. contiguous RNA evidence that was assembled independently from the genome assembly itself, using L_RNA_Scaffolder107. We first filtered the Trinity transcripts to only include the longest isoform of each gene and kept only the longest isoforms that were >500 bp long (n=73,422; N50 = 1,482 bp). L_RNA_Scaffolder typically uses BLAT to infer query alignments, which does not support genomes larger than 4 Gbp. We therefore instead used NCBI BLAST v2.2.31+108 to build a database:
makeblastdb -dbtype nucl -in <reference>
We then queried both sets of transcripts using MEGABLAST:
blastn -query <transcripts.fa> -task megablast -db <reference> -outfmt 5 -perc_identity 80 -max_target_seqs 100 -num_threads 36
We then converted the XML output into the expected PSL format using UCSC blastXmlToPsl v412 (http://hgdownload.soe.ucsc.edu/admin/exe/). L_RNA_Scaffolder joined contigs with shared transcripts to generate 12,277 scaffolds. The scaffolds had N50 = 158,477 bp, compared to N50 = 44,758 bp among the unscaffolded contigs (n=553,073). Scaffolds and contigs were joined into genome assembly v7.

Assembly v8: Scaffolding the contigs with short-read RNA-seq data with a first pass of BESST_RNA. We next used BESST_RNA (https://github.com/ksahlin/BESST_RNA) to perform additional RNA-based scaffolding using paired-end RNA-seq data. We first mapped the 755 M processed Illumina short-read RNA-seq read-pairs derived from the reference specimen itself and two published specimens (one male with ID TI787-MN-M; one female with ID TI787-MN-M; Supplementary Table 1; Supplementary Fig. 2)94 against the assembly with HISAT2 v2.2.1109 (mapping rates 81–86%) and merge the data into a single BAM file using samtools. We then performed scaffolding with BESST_RNA requiring three or more witness links with a maximum of 200,000 bp distance (“-e 3 -T 200000 -k 500 -d 1 -z 10000 -g 1”), which resulted in a total of 15,519 RNA-based scaffolds (including those from the previous step) and reduced the number of sequences to 557,087 in genome assembly v8.

Assembly v9: Scaffolding the contigs with long-read DNA data using FAST-SG + ScaffMatch. We hypothesized that there could be residual unresolved long-range contact information in the PromethION DNA data and used the Fast-SG scaffolder95 with KMC v3.0.0110 for k-mer counting to accomplish further scaffolding. First, we extracted a subset of long-reads ≥ 5 kbp and with average Q-scores above 10 from the read pool. Next, we generated synthetic mate-pairs in Fast-SG with insert lengths ranging from 4,000 to 40,000 bp (n=10) using 32 bp kmer matches spaced along reads and the genome to enable downstream detection links between contigs. The mate-pair settings were specified in a single-line read configuration file:
longont <reads.fq> 4000,5000,6000,7000,8000,10000,15000,20000,30000,40000 1
Fast-SG was run as:
FAST-SG.pl -k 32 -l <read_configuration.txt> -r <genome_assembly.fa> -p results -c $KMC -t 40
This produced two SAM files per synthetic mate-pair group, containing perfect forward/reversed kmer matches along reads and contigs at the approximate insert lengths specified for each insert group. The empirical mean and standard deviation of the insert lengths was computed from the first 1 M observations in each group as per the FAST-SG manual (https://github.com/adigenova/fast-sg/wiki/Hybrid-scaffolding-of-NA12878).

The statistics and SAM mappings were next processed in ScaffMatch v0.9111 to generate and resolve a scaffolding graph in order to scaffold the genome assembly:
python 2.7 $PATH/scaffmatch.py -w ONT -K 32 -c <genome_assembly.fa> \
-s 1461,1734,1907,2137,2194,2646,3736,7711,10503,13425 #std-dev \
-i 3957,4964,5965,6958,7938,9921,14634,16513,26175,34772 #mean \
-p fr,fr,fr,fr,fr,fr,fr,fr,fr,fr \
-m 1 <ont.I4000.fwd.sam>,<ont.I5000.fwd.sam>,… \
-m 2 <ont.I4000.rev.sam>,<ont.I5000.rev.sam>,… \
2>&1 | tee scaffmatch.log
Fast-SG+ScaffMatch produced a scaffolded genome assembly spanning 247,751 sequences (N50 = 142,666 bp), out of which 131,888 were scaffolds and 115,863 were contigs, and was taken as genome assembly v9.

Assembly v10: Scaffolding with linked-read DNA data using Scaff10x. Having produced scaffolds with RNA data and long-reads and reduced the overall number of sequences in the genome assembly, we next explored the possibility to scaffold sequences based on signatures of shared 10x Chromium linked-read barcodes between the outer edge regions of scaffolds and contigs. Scaffolding was performed using Scaff10X v4.2 (https://github.com/wtsi-hpag/Scaff10X). First, we mapped the linked-reads to the genome assembly using BWA, as in previous steps. For performance and memory reasons, we next pre-filtered the BAM file to only include reads mapping to the outer 20 kbp regions of scaffolds or contigs using samtools and a custom Perl script. To reduce the influence of ambiguous mappings to highly repetitive sequences in the outer regions, we also filtered the reads to only include read-pair mappings with MAPQ ≥ 20 in the same script. We then ran one iteration of Scaff10x, and configured the program to require at least six reads per barcode (“-reads 6”), eight shared barcodes as evidence of links between sequences (“-link 8”), edge length to match our estimated linked-read lengths (“-edge 20000”), with others settings at their default values:
scaff10x -nodes 20 -gap 100 -longread 1 -reads 6 -link 8 -edge 20000 -plot <mappings.bam>.png -bam <mappings.bam> <genome_assembly.fa>
This resulted in scaffolded genome assembly v10, spanning 219,207 sequences with N50 = 214,496 bp, including 108,310 scaffolds.

Assembly v11: Scaffolding with short-short read RNA-seq data with a second pass of BESST_RNA. After scaffolding the genome assembly with long-read and linked-read DNA data, we applied a second round of RNA scaffolding using BESST_RNA and RNA-seq data, as above, to find gene-based links that were not resolved in previous steps. BESST_RNA reported that 2,294 new scaffolds had been formed. This round of scaffolding further reduced the number of sequences from 219,207 in v10 to 216,722 (N50 = 220,983 bp) in genome assembly v11, out of which 106,444 sequences were scaffolds.

Assembly v12: Short-read polishing with Pilon using high-coverage RNA-seq data. The Illumina RNA-seq data from the reference specimen was re-mapped to the genome assembly (88.3% mapping rate) and one final round of base-level polishing Pilon (“--fix snps,indels”) was applied across transcribed regions with at least 100× depth of coverage (“--mindepth 100”), resulting in sequences with minor adjustments that were taken as genome assembly v12.

Assembly v13: Finishing the assembly by removing contaminants and re-inserting contigs. Some qualitative adjustments were carried out to reassign sequences between the main genome assembly and other classes of sequence in order to maximize gene completeness and remove putative non-krill contaminations and mitochondrial artifacts. We screened the main assembly (v12), as well as the “artifact” (~565 Mb) and “haploid” (~1.95 Gb) sequences classified by purge haplotigs (produced in genome assembly v5) to identify candidate sequence for removal or inclusion.

Reintroduction of krill genes into the main assembly. We identified missing genes in the main assembly using both comparative and experimental evidence. By running BUSCO on the artifacts and haplotigs we identified two complete single-copy BUSCO genes that had no corresponding BUSCO tblastn hit in the main assembly. In addition, we used GMAP to identify ten contigs from the haplotig class that had unique Trinity transcripts mapped with high quality (>90% identity across >90% of the transcript) but with no hits in the main assembly. We considered these 12 contigs to have been inadvertently misclassified and transferred them back into the main assembly.

Ribosomal gene sequences occur in high copy numbers in many metazoans and may be difficult to target accurately with short-read data. We therefore used the nucleotide-nucleotide BLAST (blastn) and the extensive SILVA LSU and SSU ribosomal RNA databases (release 132)96,112 to identify artifact sequences with ribosomal loci from the krill. In the majority of contigs where hits were detected, the top scoring krill alignment was against the expected M. norvegica templates (LSU accession AY744900.1.2423 [2,423 bp]; SSU accessions AY781434.1.1853 [1,853 bp] or DQ900731.1.1884 [1,884 bp], respectively), although sometimes templates for other krill had better matches (typically from the Antarctic krill Euphausia superba), which occurred in particular when matches were short and fragmented. For example, without considering the quality and lengths of the reported hits against the main assembly and the artifacts, the M. norvegica LSU template had the best krill match in n = 70/96 [73%] of cases, while the SSU templates were the best krill matches in 15/27 [56%] of cases. Requiring hits to be at least 1,500 bp with 98% identity against the template raised the on-target M. norvegica best-hits to 91% (n = 10/11) and 100% (n = 8/8) for LSU and SSU. We found that most sequences with long motifs matching M. norvegica ribosomal templates (>1,500 bp; >90% identity) had been placed into the artifact class (LSU: n = 14 vs n = 5; SSU: n = 10 vs n = 5), likely due to aberrant mapping depths across these multi-copy loci. We transferred these artifact contigs (n = 18) with long matches back into the main assembly.

Detection of contaminants. We searched the SILVA LSU and SSU ribosomal RNA databases to identify possible eukaryotic or bacterial contaminants using blastn and considering the top-scoring alignment for each assembly sequence. Among contigs flagged as artifacts, we detected long and high-identity ribosomal matches (>1,500 bp and >90% identity) against Enterobacter ludwigii/cloacae, Lactobacillus reuteri and Delftia tsuruhatensis. Querying the Nanopore sequence data itself against the SILVA database (see below), we re-classified the Delftia matches as Delftia acidovorans. We also identified several sequences in the main assembly with partial SSU low-identity hits against Aphelidium desmodesmi (n = 204; average identity = 81%), a eukaryotic endoparasite of freshwater green algae. However, this SSU accession (KY249641.1.3521 [3,521 bp]) also contains a partial LSU sequence with relatively high identity towards euphausiids (~85% identity across 550–650 bp). A combined search using both SSU and LSU simultaneously completely masked the hits against this accession in favor of M. norvegica LSU hits.

We downloaded the corresponding bacterial genome sequences from NCBI Genome Database (CP027618.1 for Enterobacter cloacae [5.0 Mbp genome]; CP015408.2 Lactobacillus reuteri [2.0 Mbp] and CP019171.1 for Delftia acidovorans [6.6 Mbp], respectively). We then used blastn to screen the krill genome assembly for contigs and scaffolds matching these sequences. Among artifact contigs, we recovered multiple high-scoring alignments between assembly sequences and the Enterobacter (n = 68; average alignment length = 48 kb; average identity = 92.5%) and Delftia (n = 75; average alignment length = 44 kb; average identity = 92.4%) genomes, respectively. Among the sequences in the main assembly, we detected less Enterobacter (n = 7; average alignment length = 24 kb; average identity = 99.2%; genome coverage = 3%) and Delftia (n = 12; average alignment length = 88 bp; average identity 96.2%; genome coverage = 4%). Matches against the Lactobacillus tended to score lower, covering only 0.6% of the bacterial genome among artifact sequences (n = 7; average alignment length = 1.9 kbp; average identity 86%) and 2% among main assembly sequences (n = 17; average alignment length = 2.1 kbp; average identity = 97%).

We then queried up to 1 M Nanopore reads produced by each PromethION flow cell and sequencing-run (original or flushed) against each bacterial genome (in total, 14.8 M reads out of 84.4 M reads) using blastn and estimated that we had sequenced about 73 k reads from Enterobacter for a total of about 60 Mb or 12× average coverage across its genome, 66 k reads of Lactobacillus (~82 Mb; 40× coverage) and 14 k reads of Delftia (~180 Mbp; 27× coverage). Taken together (assuming no overlap between the datasets) we estimate that about 0.02% of reads and 0.05% of data was derived from these bacteria rather than the krill itself. Because of the low frequencies and sporadic appearance of reads among flow cells, we hypothesize that the bacteria were technical contaminants rather than being naturally associated with the krill itself. Delftia acidovorans has previously been associated with biofilm contamination of Blue Pippin size selection cassettes113. We used a basic blastn criterion to move main and artefact assembly sequences with at least one top-scoring alignment with >90% identity towards any of the bacteria across 500 bp into a “contamination” class (Enterobacter n = 51; Lactobacillus n = 9; Delftia n = 54), spanning about 98–99% of aligned sequence.

Flagging residual mitochondrial sequence. We had already screened the Nanopore dataset for mitochondrial reads against a draft assembly of the mitochondrial chromosome. Nevertheless, we probed the genome assembly for sequences matching the M. norvegica mitochondrion with >90% identity across at least 2 kb (longer than any individual mitochondrial gene) and detected 15 sequences harboring either mitochondrial genes or fragments of the AT-rich control region. Three of these were scaffolds, each of which with a single contig containing the long match but with short matches on one or more other contigs. We split these scaffolds and extracted only the offending contig while transferring the remaining contigs back into the assembly. We can not exclude that some of these sequences are nuclear mitochondrial DNA (NUMTs) rather than misassembled sequences and therefore put these sequences into their own mitochondrial artifact class.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "• What are the noteworthy results?\nThis is a notable addition to our understanding of the biology of Northern krill, providing a number\nof insights into the evolution and adaptive potential of this key species in the pelagic community of\nthe North Atlantic, Mediterranean and Arctic seas. The paper examined genome scale variation in\nspecimens taken from most of the extremities of the distributional range of this species. Although\nthere was homogenising gene flow, there was also evidence of adaptive divergence across a\nnumber of genes, principally involved in light perception, circadian regulation and moulting\namongst others. A lipid transporter gene was found to exhibit a high degree of variation indicating\nlocal adaptation – the authors speculated this may allow the species to respond to the early spring\nconditions in the Mediterranean. There was also a number of interesting insights into lineage and\nprevious population size although some results were equivocal.\n• Will the work be of significance to the field and related fields? How does it compare to the\nestablished literature? If the work is not original, please provide relevant references.\nI am not a molecular biology specialist although have a strong background in the biology and\necology of Northern krill. I found the paper to be adequately written in terms of drawing out the\nsignificance of the highly technical results on various genomic analyses. The fact that variation in\nthe molecular results related to features known to be important in the biology of euphausiids, such\nas photoreception, moulting and thermal tolerance, added a high level of significance to the\nfindings and contributes to a new level of understanding on the major drivers and limitations to\neuphausiid behaviour and population ecology. The work is highly original and makes a major step\nforward in our understanding of the molecular biology of this species.\n• Does the work support the conclusions and claims, or is additional evidence needed?\nAlthough my understanding of the molecular biology aspects of this work has limitations, I found\nthe work to be thoroughly implemented and rigorous. The conclusions are consistent with the\nresults presented throughout the work, with some caveats (see below).\n• Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit\npublication or require revision?\nI do not disagree with the conclusions but will point out to the authors that maybe further\nexplanation is required as to why effective population size (NE) appears to be such an\nunderestimate of the likely true population size (line 196) – although they provide some\nexplanatory text (lines 197-202), they do not fully explore to what degree this error calls\nassociated analyses into question.\nFig 2 cii – this is hard to interpret without a colour code – and why are the blocks for Spain without\nvariation (please include an explanation in the legend)\nLine 331 – from a biological perspective, it is not clear why moulting in euphausiids is so frequent\ncompared to other crustaceans. One potential driver is to reduce external parasitism causing drag,\nwhich is worth mentioning here (Tarling, Geraint A., and J. Cuzin-Roudy. \"External parasite\ninfestation depends on moult-frequency and age in Antarctic krill (Euphausia superba).\" Polar\nBiology 31 (2008): 121-13)\nLine 368 – Velsch et al 1994 found that Northern krill had a 23 h circadian rhythm in constantly\ndark conditions (Velsch, J. P., & Champalbert, G. (1994). Rytmes d'activité natatoire chez\nMeganyctiphanes norvegica (Crustacea, euphaausiacea). Comptes rendus de l'Académie des\nsciences. Série 3, Sciences de la vie, 317(9), 857-862.)\nLine 383 – it is not made clear why Scandinavian stocks in particular serve as a source of genetic\ndiversity – this argument needs to be elaborated on much more clearly within the main body of\nthe manuscript to be a concluding remark here\nLine 395 – it would be helpful for the reader if the authors could specify here which candidate\ngenes to focus on for such monitoring\n• Is the methodology sound? Does the work meet the expected standards in your field?\nI am not qualified to comment on all aspect of the molecular biological methods. I found the\nannotation of genes and their relationship to biological and ecology function to be sound\n• Is there enough detail provided in the methods for the work to be reproduced?\nThe manuscript is thoroughly written and supported by extensive Supplementary material", "reviewer_2": "The new availability of a reference genome for the euphausiid, Meganyctiphanes norvegica, is an\nextremely useful resource for the research community, including oceanographers and molecular\ngeneticists. The krill species has been a frequent focus of genetic, genomic, and transcriptomic\nstudies for exactly the reasons the authors describe: the species’ successful occupation and\nadaptation to an exceptional range of environmental / oceanographic conditions.\nThis study is a significant contribution and major step forward using long-read sequencing to\nproduce a new reference genome, which is both exceptionally large and uniquely structured, and\nhas drawn the attention of the research community for many years.\nNew evidence is reported regarding several persistent questions, including the role of present-day\nmigration (i.e., dispersal) and natural selection versus demographic history in determining\npatterns of population genetic/genomic diversity and structure of this unique species.\nThe detailed analysis and whole-genome population resequencing, transcriptomic and epigenetic\nanalyses provide solid evidence for novel conclusions about likely patterns of adaptation,\nevolution, and future responses to climate change. Of particular interest is identification and\nanalysis of “gene families” that may play significant roles in the successful adaptation of the\nspecies to diverse environmental conditions across the geographic range.\nThe methodology is sound and state-of-the-art, meeting and exceeding standards for comparative\ngenomic/transcriptomic analysis of non-model species. The methods are described in sufficient\ndetail with necessary terminology (albeit with inevitable use of jargon and acronyms) to allow\nresearchers in the field to repeat the analyses as desired, while still providing understandable\nsummaries to non-specialists.\nIn summary, this is an outstanding and invaluable contribution to a fascinating question in the\nfield of marine molecular evolution, focused on a zooplankton species that has captured the\nattention of researchers for many years."}, "author_response": "Reviewer #1 (Remarks to the Author):\n• What are the noteworthy results?\nThis is a notable addition to our understanding of the biology of Northern krill, providing a number\nof insights into the evolution and adaptive potential of this key species in the pelagic community of\nthe North Atlantic, Mediterranean and Arctic seas. The paper examined genome scale variation in\nspecimens taken from most of the extremities of the distributional range of this species. Although\nthere was homogenising gene flow, there was also evidence of adaptive divergence across a number\nof genes, principally involved in light perception, circadian regulation and moulting amongst others.\nA lipid transporter gene was found to exhibit a high degree of variation indicating local adaptation –\nthe authors speculated this may allow the species to respond to the early spring conditions in the\nMediterranean. There was also a number of interesting insights into lineage and previous population\nsize although some results were equivocal.\nThis is a good summary of the work and we hope to clear up the equivocal results below.\n• Will the work be of significance to the field and related fields? How does it compare to the\nestablished literature? If the work is not original, please provide relevant references.\nI am not a molecular biology specialist although have a strong background in the biology and ecology\nof Northern krill. I found the paper to be adequately written in terms of drawing out the significance\nof the highly technical results on various genomic analyses. The fact that variation in the molecular\nresults related to features known to be important in the biology of euphausiids, such as\nphotoreception, moulting and thermal tolerance, added a high level of significance to the findings\nand contributes to a new level of understanding on the major drivers and limitations to euphausiid\nbehaviour and population ecology. The work is highly original and makes a major step forward in our\nunderstanding of the molecular biology of this species.\nThank you! We appreciate this analysis and view.\n• Does the work support the conclusions and claims, or is additional evidence needed?\nAlthough my understanding of the molecular biology aspects of this work has limitations, I found the\nwork to be thoroughly implemented and rigorous. The conclusions are consistent with the results\npresented throughout the work, with some caveats (see below).\nAgain, we appreciate this opinion and hope to resolve the caveats below.\n• Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit\npublication or require revision?\nexplanation is required as to why effective population size (NE) appears to be such an underestimate\nof the likely true population size (line 196) – although they provide some explanatory text (lines 197-\n202), they do not fully explore to what degree this error calls associated analyses into question.\nThis is a fair point and we agree it deserves some more attention. While we would indeed expect\nlarge census population sizes (N) to be associated with large effective population sizes (N ), it is well\nc e\nrecognized that N as inferred from genetic data is often much smaller than N. In the neutralist\ne c\npopulation genetic framework, the effective population size is the size of an ideal population (that\nmeets all Hardy-Weinberg assumptions) of reproductive individuals that over time would lose\ngenetic variation through genetic drift at the same rate as the observed population. In this sense, N\ne\nis typically a long-term measure and does not reflect the current-day number of reproductives. For\nexample,N is often said to be ~10,000 in humans, many orders of magnitude lower than our current\ne\nN. A common explanation for low N /N is that populations have expanded from severe bottlenecks\nc e c\ncharacterized by low N and loss of genetic variation due to increased genetic drift. Under cycling\ne\nconditions between large and small populations, long-term N equals the harmonic mean of N\ne e\nacross generations, making it strongly influenced by the smallest population sizes over the time\nperiod. However, many other factors could contribute (see below). A low N /N ratio is common in\ne c\nmarine species, see for example Hedgecock & Pudovkin 2011. Even the superabundant marine\nphytoplankton species Emiliania huxleyi appears to have a modest N , estimated to only N 2.7\ne e\nmillion (Krasovec et al. 2020).\n∼\nWe provided two inferences of N (Fig 2b):\ne\n1. a long-term inference under mutation-drift equilibrium for an idealized population, based on\nthe mutation rate and observed levels of genetic variation in the population, which produces\na single estimate of about 1.53 million under these analytical conditions\n2. an inference using coalescent modeling with PSMC/MSMC, showing signatures of\ndemographic fluctuations around the long-term estimate that is consistent with a recent\nexpansion from a historical bottleneck, and which seem to increase beyond 2 million closer\nto the present.\nFor our analyses, we assume a mutation rate from snapping shrimp and use our empirical SNP\ndataset to estimate levels of genetic variation and infer N . We note that our estimates of N are\ne e\ncomparable to those made from WGS data from the Antarctic krill by Shao et al. (2023) (~2 million)\nand those made by Choquet et al. (2023) for multiple species using more limited transcriptome data,\nincluding for M. norvegica (~1.8 million). Those works used a similar analytical framework as here.\nShao et al. suggested that the ratio of N /N in the Antarctic krill could be as low as 5×10−9\ne c\n(2×106/4×1014). Assuming N of trillions also in the Northern krill, it could be very low in this species\nc\ntoo. In theory, there is a risk that our analyses underestimate N , for example if the mutation rate\ne\nused is too high and/or the observed levels of variation are too low. Moreover, simulations suggest\nthat estimation of N in marine species with very large populations is difficult from limited point\ne\nsamples and likely to provide artificially low N /N ratios (Waples 2016). Lastly, because mutations\ne c\naccumulate slowly, neither the equilibrium approach or coalescent approach model effective\npopulation size close to the present, but provide long-term/historical estimates.\nShao et al. did) that the low N /N ratio could also reflect biological processes. Charlesworth &\ne c\nJensen (2022) list a number of potential factors that may shape genetic variation and influence its\nrelationship with N: neutral processes such as mutational biases and a negative association\nc\nbetween mutation rate and population size (the drift-barrier hypothesis). Additionally, the effects of\ndemography could manifest from skewed distributions of offspring numbers (i.e. “sweepstakes\nreproductive success”), extinction and recolonization of local populations, and fluctuations in\npopulation size (e.g. bottlenecks). Moreover, background selection and recurrent selective sweeps\ncould also drive genetic changes at linked sites within populations.\nWe observe clearly reduced levels of genetic variation up to about 50-100 kbp away from genes on\naverage, which could be an indication of the effects of linked selection. Assuming ~30 kbp protein\ncoding genes and an effect extending on average 75 kbp upstream and downstream of each gene\nlocus, the effect of reduced genetic variation through linked selection could span almost 5 Gbp\nacross the 19 Gbp genome. As stated above, we also detect signatures consistent with a population\nbottleneck.\nThe simple equation N = θ /4µ could be unrealistic for this kind of species. It “assumes mutation-\nE W\ndrift equilibrium, which is only reached after a number of generations of the order of N . As N\ne e\nincreases, the assumption that the considered population has been devoid of bottlenecks and\nsweeps during the last N generations becomes less and less plausible” (Galtier & Rousselle 2020).\ne\nWith a 1-year or so generation time, it is possible that genetic variation in species like krill is\ninfluenced by many of the factors above, and that populations never reach the equilibrium levels of\ndiversity before evolutionary and demographic processes perturb the variation. The challenges of\naccurately estimating effective population size from genetic variation and uncovering the processes\nand life history traits that influence it, and reconciling observed levels of genetic variation with\ncensus population sizes, are active research fields with many perspectives and methodologies (see\nfor example Buffalo 2021 and Waples 2023).\nTo accommodate the reviewer’s concern we have considered the following:\n1. Patterns of linkage disequilibrium (LD) and chromosomal segments of identity by descent\n(IBD) can also offer insights into population history. Long blocks of IBD suggest recent\ninbreeding or low N , and short blocks indicate large contemporary N . Likewise, the\ne e\nspectrum of the rates of decay of LD between SNPs at different genetic distances is\ninfluenced by how N has fluctuated in the recent past. While mutations accumulate slowly,\ne\nrecombination events are frequent and can be used to estimate recent N (e.g. Santiago et\ne\nal. 2020).\nWe have now performed an inference of near-contemporary N using LD patterns with the\ne\nprogram GONE, looking at three regional gene pools. This tool suggests a recent N of 4-5\ne\nmillion, consistent with our inferences of population expansion from a bottleneck (but still\nfar below the census population size of trillions). We have written a short addition to the\n3\nmethods section (lines 1829-1838; adding one reference), put the results into\nSupplementary Figure 19 and briefly expanded the results section:\nLines 191-192: “The rate of decay of linkage disequilibrium between SNPs in turn suggested\nN may recently have reached 4–5 million (Supplementary Fig. 19B)”\nE\n2. Population bottlenecks followed by an expansion may generate an excess of rare alleles,\nsimilar to selective sweeps. If populations have very different demographic histories, this\ncould complicate interpretation of evolutionary processes. Our inferences of population\nmutation rate (θ ) and nucleotide diversity (π) in each population (Supplementary Figure\nW\n19A) are very similar, suggesting they may have a shared recent evolutionary history.\nWe have now added per-population Tajima’s D to these analyses, which indicates a similar\nlevels of excess of low-frequency variants compared to neutral expectation among all\npopulations, for example due to similar demographic histories. Tajima’s D is slightly less\nnegative in the Mediterranean sample, but the difference is small and we argue more\nsamples are needed to distinguish between Atlantic and Mediterranean patterns. The results\nof these analyses have been added to Supplementary Figure 19A.\n3. Although we estimate N in this work, not so many downstream analyses actually directly\ne\ndepend on the numerical estimate itself. The exception is the estimation of the\nrecombination rate (using iSMC) followed by the estimation of the age distribution of minor\nalleles associated with divergent loci (using GEVA), which explicitly use estimates of N to\ne\nderive a recombination frequency or sets of allele ages, respectively.\nIn the case of recombination rate, we estimate it to 0.32 cM/Mb using the re-arranged\nequation: r = ρ / ( 4 * N ), in which we inferred ρ to be 0.013 in iSMC and long-term N to be\ne e\n1.02 million from the heterozygous genotypes in\nthe deeply sequenced reference individual alone.\nSince this recombination rate was inferred from a\nsingle individual, we re-estimated θ and N for\nW e\nthis analysis from that individual alone, in order to\navoid using estimates not well-represented by that\nsubset of data. In the case of GEVA, we used the\nglobal long-term estimate of N (1.53M).\ne\nIn principle, a different N could be used, such as 4\ne\nmillion, which makes the recombination rate per\nbase lower, i.e. r = 0.085 cM/Mb (an extremely\nlow recombination rate). To test the effect that a\nlarger N could have, we re-estimated the age\ne\ndistribution of minor alleles at divergent loci\nspecifying both a lower recombination rate (r = 0.085 cM/Mb) and a larger effective\npopulation size (N =4 million) under the joint method in GEVA, combining estimates based\ne\non both mutation rate and recombination rate. Instead of the age of minor alleles being\n4\ncentered around 1-2 million years (Fig 3e), they are inferred to be older, being centered\naround 4 million years (attached fig). However, since this very large N could be recent and\ne\ndoes not appear to be supported by the data as representative of historical N , we argue\ne\nthat the long-term estimate as originally used is better suited for this analysis, and we\nchoose to keep it as is. Moreover, our statement that most variants in the putatively\nadaptive loci appear to predate multiple glacial cycles appears to be true under both\nscenaria.\nWe acknowledge that estimating N is difficult in species like krill. We have applied several methods\ne\naimed to estimate N over different time periods. We do not think our estimates of much lower N\ne e\nthan N are erroneous per-se. They could reflect evolutionary and demographic history and possibly\nc\nalso specific aspects of krill biology.\nIn addition to the edits above, we have added a section in the discussion bringing up why we may\nsee low N /N ratio in krill, citing previously used krill references as well as the Charlesworth &\ne c\nJensen 2022 review paper. We speculate that the mode of reproduction in krill could limit N due to\ne\nskewed reproduction, but acknowledge that more research is needed to understand krill population\ndynamics. This text is written across lines 299-308.\nWe hope our response here and revised text provides enough context and explanation to be\nacceptable.\n Hedgecock & Pudovkin 2011:\nhttps://www.ingentaconnect.com/content/umrsmas/bullmar/2011/00000087/00000004/ar\nt00013\n Krasovec et al 2020: https://academic.oup.com/gbe/article/12/7/1051/5869440\n Waples 2016: https://onlinelibrary.wiley.com/doi/10.1111/jfb.13143\n Charlesworth & Jensen 2022: https://academic.oup.com/gbe/article/14/7/evac096/6615374\n Galtier & Rousselle 2020: https://academic.oup.com/genetics/article/216/2/559/6066183\n Buffalo 2021: https://elifesciences.org/articles/67509\n Waples 2023: https://onlinelibrary.wiley.com/doi/full/10.1111/1755-0998.13879\n Santiago et al 2020:\nhttps://academic.oup.com/mbe/article/37/12/3642/5869049?login=true\nvariation (please include an explanation in the legend)\nIn Fig 2 cii, we show an admixture plot. In this analysis, a model-based framework is used to assign\nthe observed genetic variation to one or more hypothetical ancestral pools or populations (given by\nthe K number). This method simultaneously estimates and optimizes allele frequencies for ancestral\npopulations and ancestry proportions for individuals from genotype data. It is rerun with different\nvalues of K, which is used to reveal genetic structure at different levels of resolution. Each bar\nreflects the overall proportions to which different hypothetical ancestors have contributed to the\ngenetic fingerprint of an individual.\nTo perform this analysis, we ran the tool Admixture 50 times to analyze a subsample of LD-pruned\nSNPs, essentially producing 50 pseudoreplicate runs with different starting conditions for some\nparameters but using the same data. We then summarized the results using CLUMPAK, which\nidentifies the major clustering modes across the 50 replicate runs (which individually may have had\nslightly more ambiguity in the ancestry assignment than the summary output).\nWe used the summary output for the figure and have avoided re-using the same color scheme for\nthis plot that we used for other figures, such as Fig 2 ci & Fig 2 ciii, as these hypothetical ancestors\ndo not necessarily reflect geographic structure. Samples from Spain stand out for having genetic\nvariation that is comparably easily attributed to one of the ancestral populations. Their distinct\ngenetic variation is also reflected in figures 2 ci and 2 ciii. We see no correspondingly strong genetic\nstructuring between SW and NE Atlantic samples (consistent with the PCA overlap in i and low F\nST\nseparating them in ii).\nWe have updated the legend to better explain the meaning of the admixture plot in Fig 2 cii and in\nthe corresponding supporting figure (Supplementary Figure 20).\nLine 331 – from a biological perspective, it is not clear why moulting in euphausiids is so frequent\ncompared to other crustaceans. One potential driver is to reduce external parasitism causing drag,\nwhich is worth mentioning here (Tarling, Geraint A., and J. Cuzin-Roudy. \"External parasite\ninfestation depends on molt-frequency and age in Antarctic krill (Euphausia superba).\" Polar\nBiology 31 (2008): 121-13)\nWe agree this deserves an explanation. We now state that the frequent molting in euphausiids is\nalso a means to reduce load and drag from parasites and epibionts in the discussion. We quote both\nTarling et al. as suggested and the chapter on parasites and disease by Gómez-Gutiérrez & Morales-\nÁvila in Biology and Ecology of Antarctic Krill, in order to link molting physiology to drag.\nIn the discussion we thus write:\n“Molting is a crucial process in krill, being interlinked with growth and reproduction and a means to\nreduce load and drag from parasites and epibionts57,58” (lines 291-292)\ndark conditions (Velsch, J. P., & Champalbert, G. (1994). Rytmes d'activité natatoire chez\nMeganyctiphanes norvegica (Crustacea, euphaausiacea). Comptes rendus de l'Académie des\nsciences. Série 3, Sciences de la vie, 317(9), 857-862.)\nThank you for this remark. We have included this reference. Together with observations in the\nspecies we already brought up, this observation suggests that the short circadian rhythm is an\nancestral feature in krill. We have edited this paragraph to include the reference and reflect that this\nis a shared trait among krill.\n“M. norvegica, along with other krill species such as E. superba and Thysanoessa inermis, exhibit\nendogenous circadian rhythms (ECR) that cycle at rates faster than 24 hours in the absence of\nlight70,71, or respond to minute irradiance.” (lines 334-336)\nIn this revision, we are replacing the original Biscontin et al. 2017 reference with Biscontin et al.\n2019, which addressed circadian rhythm in E. superba more comprehensively, at the transcriptome\nlevel.\nLine 383 – it is not made clear why Scandinavian stocks in particular serve as a source of genetic\ndiversity – this argument needs to be elaborated on much more clearly within the main body of\nthe manuscript to be a concluding remark here\nWe appreciate the issue and have clarified these observations. When examining divergent loci, we\nfound that Mediterranean haplotypes were relatively common in the Norwegian sample (mean\nfrequency of 13%), while South-Western haplotypes were common in the Swedish sample (mean\nfrequency of 17%). In contrast, the other Northern stocks do not harbor haplotypes from Southern\nregions to the same extent. These results are shown in Figs. 3diii and Figs. 3dvi, respectively. We\nhave now described these results more clearly in the text:\n“We found that Mediterranean haplotypes were comparably common in the Norwegian sample\n(mean frequency of 13%, Fig. 3diii), while South-Western haplotypes were frequent in the Swedish\nsample (mean frequency of 17%, Fig. 3dvi), suggesting these Scandinavian stocks in particular could\ncontain genetic material that is otherwise rare at high latitudes.” (lines 229-232)\nLine 395 – it would be helpful for the reader if the authors could specify here which candidate\ngenes to focus on for such monitoring\nWe agree. Recommended markers for monitoring would be those that we observe to be\ngeographically divergent and ecophysiologically relevant. As an example, we have now highlighted\nsix genes linked to thermal stress or reproduction towards the end of the discussion.\nWe now write:\n“The many candidate genes reported here can be used as biomarkers to diagnose and monitor\nchange of adaptive variation in response to changing conditions, including CCT8, HIP-R and Hsp83\n7\nassociated with thermal stress response and nrf-6, Kr-h1 and Hr3 involved in reproductive\nphysiology, as well as those potentially coordinating photoperiodic regulation. These markers could\nalso be used to monitor other species, in order to better forecast the future distributions of krill and\nestimate risks of the great many species that depend on them.” (lines 362-367)\n• Is the methodology sound? Does the work meet the expected standards in your field?\nI am not qualified to comment on all aspect of the molecular biological methods. I found the\nannotation of genes and their relationship to biological and ecology function to be sound\n• Is there enough detail provided in the methods for the work to be reproduced?\nThe manuscript is thoroughly written and supported by extensive Supplementary material\nThank you for the careful and insightful review of this manuscript!\nWe thank the reviewer for the thorough and constructive review. We believe that the issues that\nthey underscored have been directly relevant to central topics in our manuscript, and that\naddressing them have helped improve its quality. We hope that our replies and revisions have\nsatisfyingly resolved each of the highlighted issues.\nThe new availability of a reference genome for the euphausiid, Meganyctiphanes norvegica, is an\nextremely useful resource for the research community, including oceanographers and molecular\ngeneticists. The krill species has been a frequent focus of genetic, genomic, and transcriptomic\nstudies for exactly the reasons the authors describe: the species’ successful occupation and\nadaptation to an exceptional range of environmental / oceanographic conditions.\nThis study is a significant contribution and major step forward using long-read sequencing to\nproduce a new reference genome, which is both exceptionally large and uniquely structured, and\nhas drawn the attention of the research community for many years.\nNew evidence is reported regarding several persistent questions, including the role of present-day\nmigration (i.e., dispersal) and natural selection versus demographic history in determining\npatterns of population genetic/genomic diversity and structure of this unique species.\nThe detailed analysis and whole-genome population resequencing, transcriptomic and epigenetic\nanalyses provide solid evidence for novel conclusions about likely patterns of adaptation,\nevolution, and future responses to climate change. Of particular interest is identification and\nanalysis of “gene families” that may play significant roles in the successful adaptation of the\nspecies to diverse environmental conditions across the geographic range.\nThe methodology is sound and state-of-the-art, meeting and exceeding standards for comparative\ngenomic/transcriptomic analysis of non-model species. The methods are described in sufficient\ndetail with necessary terminology (albeit with inevitable use of jargon and acronyms) to allow\nresearchers in the field to repeat the analyses as desired, while still providing understandable\nsummaries to non-specialists.\nIn summary, this is an outstanding and invaluable contribution to a fascinating question in the\nfield of marine molecular evolution, focused on a zooplankton species that has captured the\nattention of researchers for many years.\nThank you for the thoughtful and positive review of our manuscript!\nWe thank the reviewer for this feedback and hope they find our response and edits to have\nimproved it further.\nBelow we list (A) updates, (B) corrections of technical issues and (C) editorial changes.\nA. Updated analyses or reporting.\n1. Results, Main figure 1f (update) in the context of DNA\nmethylation. We have enlarged Fig 1f and put the mean DNA\nmethylation frequencies into the figure. We have also put\nthese values into the legend for Supplementary Figure 17A,\nwhich has the same mean methylation frequencies presented\nas a bar plot. In Supp Fig 17A one bar (“non-repeats”) was\nmisdrawn and has been corrected.\n2. Results, Supplementary Fig 19D (update) in the context of divergence.\nWe originally performed an genetic isolation-by-distance test using the\nd estimator against physical distance. We have realized such a test\nXY\ntypically uses the F estimator and specifically linearized F (see\nST ST\nRousset 1997, 4,000+ citations). To improve our test, we have updated\nit to instead use this standard estimator for this context.\nAs in the original analysis, genetic distance only increases marginally with physical distance.\nThe association between genetic and physical distance is weak when comparing North\nAtlantic Ocean vs. Mediterranean Sea populations (R2=0.47; P=0.09) but strong when\ncomparing North Atlantic Ocean populations (R2=0.69; P=3.346e-06), i.e. within the same\nocean basin. These results match the F -tree in Fig 2ciii well. The results have been\nST\nreworded (line 197-199), a methods section has been written (lines 1900-1911) (which was\noriginally missing) and values are reported in Supplementary Data 8.\n Rousett: https://pubmed.ncbi.nlm.nih.gov/9093870/\n3. Results, Supplementary Data 8 (reporting) in the context of divergent gene haplotypes. In\nthe results, we address the haplotype distribution of 660 and 34 genes segregating among\npopulations in two major contrasts (Fig. 3d-e), respectively. Information about the genes\nwere missing. We have now added two gene lists to Supplementary Data 8 with per-gene\ninformation about the distribution of haplotypes in each population.\nWe have identified a few technical errors when preparing the revisions. We believe they do not\nmeaningfully impact results or conclusions, but propose to correct them. We apologize for making\nthese mistakes and hope the corrections are acceptable.\n1. Results, in Supplementary Fig 17C (formerly D) on DNA methylation. Bootstrap confidence\nintervals (CIs) for the P90 statistic (frequency of CpGs with >90% of reads methylated) were\nincorrectly estimated from mean methylation frequencies. We now use the correct P90 95%\nCIs, and longer whiskers to improve visibility. Methylation is structured as before (more\nmethylation in genes and repeats), but has a little more variability within regions with this\nstatistic.\nBefore: After:\n2. Results, in Fig 3bi and Supplementary Fig 19C in the context of divergence. These two\nsubpanels were based on incomplete SNP data. Compared to the correct SNP dataset (~760\nmillion SNPs) underlying all other analyses, simulations and results, this data spanned about\nbased on the complete dataset, as in the rest of the manuscript.\nIn Fig 3bi, the corrected plot is somewhat right-shifted: 1.2% instead of 0.9% of the genome\nhas F >0.4 in the Atlantic vs Mediterranean contrast (at/me), and 0.08% vs. 0.05% in the\nST\nAtlantic contrast (ea/we).\nBefore: After:\nWe have updated the main text in one place, from:\n“Divergent regions (F >0.4) spanned <1% of the genome and were about 2× enriched for\nST\ngene sequences and extended haplotypes”; to:\n“Divergent regions (F >0.4) spanned ~1% or less of the genome and were about 2× enriched\nST\nfor gene sequences and extended haplotypes” (lines 215-216)\ndistances. There are no meaningful changes to interrelationships among individuals. As\nbefore, the tree lacks structure and is star-like overall, with the Mediterranean samples\nforming the only recognizable group on a short separate branch.\nBefore: After:\n3. Results, in Figure 4 in the context of showing example candidate genes for adaptation. We\nhave updated Figure 4bi. The subpanel highlights divergence across an rdgB homolog, a\ngene involved in retinal photopigment metabolism.rdgB genes exist in several copies in the\nNorthern krill, and more than one is divergent among populations. We accidentally\nhighlighted the second most divergent gene (“REF_STRG_1_54916_XLOC_034835”) instead\nof the intended and most divergent one (“COM_SPAL_1_013733_XLOC_048553”). We have\nnow replaced the subpanel with the intended gene.\nBefore: After:\nIn 4biii, the x-axis scales were different for the upper F plot and the lower XP-nSL plot for\nST\nthe Acer gene. The two plots have now been corrected to have the same range.\nBefore: After:\nWe have also made editorial changes such as correcting references or adding explanations to\nlegends, or sample sizes, following instructions from the editors.\n1. In Results, references. Fig 1a, we use a map with Sea Surface Temperatures from Climate\nReanalyzer. As per their instructions:\nTo reference Climate Reanalyzer in a publication, please use the website citation below in\naddition to source citations for any datasets used:\nClimate Reanalyzer (n.d.). [Title of specific page]. Climate Change Institute, University of\nMaine. Retrieved [Month Day, Year], from https://climatereanalyzer.org\nWe have updated the reference as per these instructions. We have also followed the\nguidelines here: https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres.html to add\nHuang et al. (2021) to the reference list for the most recent OISST v2 citation, (adding it to\nthe top of the methods section, which is slightly revised to accommodate this change) (lines\n358-360). We have also clarified in the legend that circles indicate sample locations.\nWe have also replaced the original photo of a Northern krill (that was used with permission)\nwith one recently taken by us.\n2. In Results, Fig 3. We added the values of the respective F percentiles indicated in C i and ii\nST\nto the legend, as these were neither given in the text or the figure.\n3. In Results, Fig 4. We added the values of the respective F percentiles indicated in A, B and\nST\nC to the legend, as these were neither given in the text or the figure.\n4. In Discussion, reference. We have revised the former “Russel 1969” reference. This\nreference is for a book with the two main authors being Mauchline and Fisher (Russel is the\neditor). We have updated this reference to name the authors rather than the editor.\n5. In Results, Supplementary Fig 1. Corrected a misalignment glitch in 1C about read lengths\nwith not all N50 dots showing.\n6. In Results, Supplementary Fig 5. A minor correction: “Assembly completeness is assessed by\nlooking at the number of gene bodies covered to a certain extent at 75% identity (E).” It\noriginally said 90% identity.\n7. In Results, Supplementary Fig 17. A legend has been added to indicate the meaning of mean\n(black) and median (blue) circles in A (methylation rates in genomic regions). Sample size\nand explanation has been added to D (splice isoforms vs. DNA methylation):\n“Only gene models supported by RNA-seq data and with at least one CpG site with\n13\nmethylation data were included in the analysis (n=19,777 genes)”.\nThe number of CpG sites were added to E. The section on sample size now says:\n“Number of repeats per identity interval (left to right): n=7; n=91; n=1397; n=211. Number\nof CpG sites per identity interval (left to right): n=549; n=6,008; n=110,085; n=14,021.”\n8. In Results, Supplementary Fig 23. Renumbered subpanels (from a, b ... to i, ii …).\n9. In Results, Supplementary Fig 25. Added numbering to subpanels (i, ii …).\n10. In Results, Supplementary Data 4. Minor changes to table titles/headers in file\nSupplementary Data 4 for clarification.\n11. In Results, Supplementary Data 7. This table originally contained accessions used to find\nand characterize DNA methylation genes through comparative analyses. It has now been\nexpanded to also contain the source for the homeodomain sequences used to characterize\nHox genes. The title and headers of the file have been updated accordingly.\n12. In Methods. Updates to database version specifications. Change of “NCBI Genbank” to\n“NCBI Genome Database” for the source of the bacterial genomes used to screen the\nassembly. Specified RefSeq release version (RefSeq release 204). KrillDB specified as\nKrillDB2. Dfam specified as Dfam_3.5. HOMEODB specified as HomeoDB2. Public databases\nused are now listed under ‘Data Availability’.\n13. Supplementary Notes. The supplementary text sections about whole genome duplication\nand geographical regions of population samples are now written as Supplementary Note 1\nand Supplementary Note 2.\n14"}]

Decision letter