Reduced-representation sequencing detects trans-Arctic connectivity and local adaptation in polar cod (Boreogadus saida)

Information on connectivity and genetic structure of marine organisms remains sparse in frontier ecosystems such as the Arctic Ocean. Filling these knowledge gaps becomes increasingly urgent, as the Arctic is undergoing rapid physical, ecological, and socio-economic changes. The abundant and widely distributed polar cod (Boreogadus saida) is highly adapted to Arctic waters, and its larvae and juveniles live in close association with sea ice. Through a reduced-representation sequencing approach, this study explored the spatial genetic structure of polar cod at a circum-Arctic scale. Genomic variation was partitioned into neutral and adaptive components to respectively investigate genetic connectivity and local adaptation. Based on 922 high-quality single nucleotide polymorphism (SNP) markers genotyped in 611 polar cod, broad-scale differentiation was detected among three groups: (i) Beaufort-Chukchi seas, (ii) all regions connected by the Transpolar Drift, ranging from the Laptev Sea to Iceland, including the European Arctic, and (iii) West Greenland. Patterns of neutral genetic structure suggested broadscale oceanographic and sea ice drift features (i.e. Beaufort Gyre and Transpolar Drift) as important drivers of connectivity. Genomic variation at 35 outlier loci indicated adaptive divergence of the West Greenland, and the Beaufort-Chukchi Seas populations, possibly driven by environmental conditions. Sea ice decline and changing ocean currents can alter or disrupt connectivity between polar cod from the three genetic groups, potentially undermining their resilience to climate change, even in putative refugia, such as the Central Arctic Ocean and the Arctic Archipelago.

Methods
A total of 652 polar cod samples were collected in the Central Arctic Ocean, Fram Strait and the Arctic Ocean shelves bordering Alaska, Canada, Russia, Greenland and Northern Europe, during several expeditions between 2003 and 2021. Fish were collected using bottom trawl, bongo net, Young fish trawl, Surface and Under-Ice Trawl, ROVnet, Multpelt 832 pelagic trawl, or zooplankton net. We collected fin clips of all fish sampled and stored them in 96 % ethanol. Genomic DNA was extracted from fin clips using the NucleoSpin® Tissue kit (Macherey-Nagel), following manufacturer’s instructions. A modified version of the Elshire et al. (2011) genotyping-by-sequencing (GBS) method was used, with a single restriction enzyme (PstI) and size selection (320-720 bp).Three GBS libraries were paired-end sequenced on an Illumina Novaseq platform 6000 (PE100) at the Genomics Core Leuven (www.genomicscore.be). Sequences were quality checked using FastQC v0.11.8 and Stacks v2.5 was used to process the GBS data. Reads were aligned to a draft Boreogadus saida reference genome (GenBank accession number GCA_900302515.1) using Bowtie2 v2.3.4.3. SNPs were called using default parameters in Stacks, which employs a Bayesian genotype caller, and filtered using VCFtools and various R packages (SNPRelate v1.28.0, poppr v.2.9.3, hierfstat v.0.5-11, pegas v1.1), resulting in SNP datasets for 611 individuals. For the Redundancy analysis, we employed a set of 16 climate and biogeochemical variables, obtained from a georeferenced dataset describing seawater conditions from the Copernicus Marine Environment Monitoring Service.