Introduction
Hybridization is a common occurrence among animal species (Mallet 2005). Weak isolating mechanisms and characteristics of fish biology such as external fertilization (Campton 1987; Leary et al. 1995), rarity of one or both parental species (Frisch and van Herweden 2006; Montanari et al. 2012), and habitat overlap of parental species (Marie et al. 2007; Yaakub et al. 2007) are likely causes of the high incidence in fishes. Hybridization is often viewed as having negative consequences on biodiversity (Edmands 2007), given its potential to break up coadapted gene complexes for local adaptations (Barton and Hewitt 1989); however, introgression can also provide benefits by creating genetic and phenotypic variation that may give individuals a fitness advantage (Pardo-Diaz et al. 2012; Dowling et al. 2016; Hanemaaijer et al. 2018). Such benefits become especially relevant in regions such as the Arctic Ocean, which are facing wide-scale habitat modifications (Cronin and Cronin 2015) that are expected to change reproductive interactions among species (Garroway et al. 2010; Chunco 2014) and ultimately impact community-wide dynamics (Farkas et al. 2015; Canestrelli et al. 2017). However, predicting the outcomes and attributing drivers of hybridization (e.g., ocean warming) is difficult (Potts et al. 2014), primarily due to the lack of baseline assessments. Therefore, evaluating levels of current genetic exchange is a necessary first step toward understanding potential future outcomes, not only in terms of genomic vulnerability for species but also its potential community and ecosystem-wide impacts.
Gadids (cod; Gadidae) are a diverse group of teleost fish that inhabit a variety of marine ecosystems. This group’s ability to successfully flourish in a broad range of habitats may be due to the acquisition of genetic material through adaptive introgression or hybrid speciation (Árnason and Halldórsdóttir 2019). However, apart from Atlantic cod (Gadus morhua, e.g., Bradbury et al. 2014; Weist et al. 2019), little attention has been focused on estimating the influence of contemporary hybridization on genetic diversity within gadids at the population level. In the Arctic Ocean, the Polar cod (Boreogadus saida Lepechin, 1774; hereafter referred as Boreogadus) and Arctic cod (Arctogadus glacialis Peters, 1872; hereafter referred as Arctogadus) play an essential role in the Arctic food web (Finley and Gibb 1982; Welch et al. 1992) as they provide an energy-rich food source upon which many upper trophic level predators rely (Harwood et al. 2015). Both species have a pan-Arctic distribution and are often associated with sea-ice as both are well-adapted to cold-water environments. Assuming both species are stenothermic (evident in Boreogadus; Laurel et al. 2016), warming temperatures will likely negatively impact these species, in particular at early life stages, and their ability to compete with invading temperate species such as the eurythermal walleye pollock (Gadus chalcogrammus Pallas, 1814) (Bouchard et al. 2017; Majewski et al. 2017; LeBlanc et al. 2020; Marsh and Mueter 2020). Although Boreogadus has been extensively studied, many information gaps remain unanswered (Mueter et al. 2016). Understanding the genetic and phenotypic composition of these Arctic dwelling species and levels of interspecies genetic exchange can provide valuable insight into factors, such as habitat use and other life-history characteristics, influencing population dynamics.
Local segregation, either through temporal or ecological separation (i.e., microhabitat isolation), is thought to limit hybridization opportunities in some species. Despite a general lack of data on both species about many life-history characteristics, there appears to be some spatial segregation on the summer feeding grounds at the adult life stage. For instance, Boreogadus is more associated with the pelagic and Arctogadus with the benthic food web (Christiansen et al. 2012). However, the trophic relationships (e.g., pelagic or benthic) of these species may change markedly across seasons with shifts in sea-ice concentrating common prey in certain areas (Christiansen et al. 2012), potentially increasing opportunities for species interactions during fall and winter, when spawning most likely occurs (Craig et al. 1982; Süfke et al. 1998). Despite spatial segregation of adults during summer, there is considerable overlap in ecological niche and local distribution during early life stages (Bouchard et al. 2016). Notwithstanding, the potential sympatry on reproduction sites, there are no current reports of hybridization between these two species. However, spawning locality and behavior are not well known (Jordan et al. 2003). Regardless, Arctic marine environments are rapidly changing and may lead to changes in community-wide species interactions. Habitat disturbance is a main driving force for hybridization (Ma et al. 2019), but changes in local distributions for these two species may also result in the formation of new areas of secondary contact or the breakdown of local spatial structure, which may currently restrict hybridization events.
To assess changes in distribution and level of genetic exchange, accurate species determination is needed as this may contribute to the lack of hybridization reports. For instance, Arctogadus and Boreogadus are often misidentified, as larval and juvenile stages are almost phenotypically indistinguishable (Teletcha et al. 2006). Misidentification during field observations has likely led to the less abundant Arctogadus being mistaken for the numerically dominant Boreogadus. Thus, a suite of genetic markers that aids in the accurate determination of species is needed to monitor changes in species abundances, distributions and frequency of hybridization. Although mitochondrial DNA is sufficiently distinct, assessments of hybridization require markers located in the nuclear genome—relatively lacking for Arctogadus—since mitochondrial DNA is solely passed through the maternal lineage. To date, only a single genetic marker located in the nuclear genome (a microsatellite GMO8) has been proposed to distinguish these two species with high accuracy (Madsen et al. 2009). However, the allele size distribution of this microsatellite is not mutually exclusive as both species possess alleles in the lower size range, limiting its use in the detection of hybridization and restricting species identification to only those individuals that possess alleles in the higher fragment size range. Using double digest restriction site-associated (ddRAD) sequences from Boreogadus and Arctogadus sampled from the Alaska and western Canada Beaufort Sea and adjacent seas, we aimed to discover a set of species-specific loci to investigate the degree of genetic divergence and the current level of potential genetic exchange across species in an area of sympatry. Further for a single hybrid individual, we used whole-genome sequences to explore its genomic composition using previously published genomic data from other gadid species (Atlantic cod; Pacific cod, Gadus macrocephalus; and walleye pollock; Árnason and Halldórsdóttir 2019). As Boreogadus and Arctogadus are often found within the remote waters of the Arctic region, genetic information can provide valuable insight into reproductive segregation that may not be easily obtained through direct observations.
Methods
Sampling and DNA extraction
Genomic DNA was extracted using a DNeasy Blood and Tissue kit and following the manufacturer’s protocols (Qiagen, Valencia, CA, USA) from fin clips or tissue of 127 fish morphologically identified as Boreogadus (n = 115) and Arctogadus (n = 12; Fig. 1). Boreogadus samples collected from the Beaufort and Chukchi seas were provided by researchers conducting other research projects and samples from Bering Sea were provided by University of Washington Burke Museum Fish Collection (Wilson et al. 2018). Arctogadus samples from western Canadian Beaufort Sea were obtained from Fisheries and Oceans Canada Specimen Collection. As species determination based on morphology can be difficult, initial species identity based on mouth and general body morphology in the field was verified using mtDNA (Wilson et al. 2017). Extractions were quantified using a Modulus Microplate (Turner BioSystems, Inc.) using a Broad Range Quant-iT dsDNA Assay Kit (Thermo Fisher Scientific, Inc.) to ensure a minimum concentration of 4 ng/µL. Detailed sample information including GPS coordinates is in Wilson et al. (2018).
[Image omitted: See PDF]
ddRAD library preparation and bioinformatics
Sample preparation for ddRAD sequencing followed the double-digest protocol outlined in DaCosta and Sorenson (2014). Genomic DNA (up to 1 µg) was digested with high fidelity versions of SbfI and MspI restriction enzymes (New England Biolabs, Ipswich, MA, USA). Amplification and sequencing adapters containing unique barcode or index sequences were ligated to the sticky ends generated by the restriction enzymes. The products were then electrophoresed on 2% low-melt agarose gel and DNA fragments of size between 300 and 450 base pair (178–328 base pair (bp), excluding adapters) were selected. DNA was extracted from the gel using a MinElute Gel Extraction Kit (Qiagen) following the manufacturer’s protocol. Size-selected fragments were then PCR-amplified with Phusion high-fidelity DNA polymerase (Thermo Scientific, Pittsburgh, PA, USA) for 24 cycles, and amplified products were cleaned using magnetic AMPure XP beads (Beckman Coulter, Inc., Indianapolis, IN, USA). Quantitative PCR using an Illumina library quantification kit (KAPA Biosystems, Wilmington, MA, USA) was used to measure the concentration of purified PCR products, and samples were pooled in equimolar concentrations. A multiplexed library was sequenced as a single-end 150-bp run on an Illumina HiSeq 2500 at the Tufts University Core Genomics Facility. Raw Illumina reads have been accessioned on NCBI’s Sequence Read Archive (SRA; Bioproject: PRJNA436585; Biosample: SAMN18039800-SAMN18039926).
Raw Illumina reads were processed using a computational pipeline described by DaCosta and Sorenson (2014; http://github.com/BU-RAD-seq/ddRAD-seq-Pipeline). First, reads were assigned to individual samples based on barcode/index sequences using bcl2fastq-1.8.4 software (Illumina Inc.) by Tufts University Core Genomics facility. Preprocessing of reads was done using a custom python script (ddRAD_fastq_qc.py, Jeff DaCosta, Boston College) that removed chimera sequences (SbfI-SbfI or MspI-MspI) or reads containing two or more mismatches in SbfI recognition site and trimmed reads to either MspI concatemer or P2 adapter. Reads per sample then were collapsed into identical clusters using the CondenseSequences.py script with low-quality reads (i.e., sequences that failed to cluster with any other reads (−id setting of 0.90) and an average per-base Phred score (<20) filtered out with the FilterSequences.py script and the uclust function in usearch v.5 (Edgar 2010). Condensed and filtered reads from all samples were then concatenated and clustered with an −id setting of 0.85 in uclust as well as clusters with identical or nearly identical BLAST hit to reference genome of G. morhua (Tørresen et al. 2017) were combined. muscle v.3 (Edgar 2004) was used to align reads within each cluster, and samples within each aligned cluster were genotyped using the RADGenotypes.py script. Homozygotes and heterozygotes were identified based on thresholds outlined in DaCosta and Sorenson (2014), with individual genotypes falling into three categories: “missing” (no data), “good” (unambiguously genotyped), and “flagged” (recovered heterozygous genotype, but with haplotype counts outside of acceptable thresholds or with >2 alleles detected). Clusters were also “flagged” if the number of single-nucleotide polymorphisms (SNPs) were >8 and if >3 SNPs showed strong linkage from output produced in the genotyping step (“clustersummary.out” file). We used Geneious v10 (Biomaters Inc., San Francisco, CA, USA) to manually check and edit loci. To limit any biases due to sequencing error and (or) allelic dropout, alleles with less than five times coverage were scored as missing, such that a minimum of ten reads was required to score a locus as heterozygous. Final output files were generated with custom python scripts (available at https://github.com/BU-RAD-seq/Out-Conversions and see Lavretksky et al. 2016).
Whole-genome sequencing and bioinformatics
Whole-genome sequencing, conducted independently of the ddRAD-seq analyses, followed Árnason and Halldórsdóttir (2019). Briefly, genomic DNA was isolated from the tissue samples from the Árnason and Halldórsdóttir DNAfish collection using the E-Z 96 Tissue DNA Kit (Omega469Biotek), according to the manufacturer’s recommendation. The DNA was normalized with elution buffer to 10 ng/µL. The normalized DNA was analyzed at the Bauer Core of Harvard University. The Bauer Core used the Kapa HyperPrep Plus kit (Kapa Biosystems) for enzymatic DNA fragmentation and adapter ligation according to the manufacturer’s recommendation except that the reaction volume was one-fourth of the recommended volume. The target insert size was 350 bp with a resulting average of 487 bp. The libraries were indexed using IDT (Integrated DNA Technologies) unique dual 8 bp indexes for Illumina. The Core uses Tapestation (Agilent Technologies) and Picogreen qPCR for sizing and quality checks. Multiplexed libraries were sequenced on NovaSeq (Illumina) S4 lanes at the Broad Institute with a 2 × 150 bp read format, and demultiplexed fastq files of each individual were returned. Data processing followed the Genome Analysis Toolkit (GATK) best practices () as implemented in the fastq_to_vcf pipeline of Alison Shultz (github.com/ajshultz/comp-pop-gen). Using the pipeline, the raw reads were adapter trimmed using Ngmerge (Gaspar 2018), the trimmed fastq files aligned to the gadMor3.0 chromosome-level reference genome assembly (NCBI accession ID: GCF_902 167 405.1) using bwa mem (Li and Durbin 2009), and the resulting bam files deduplicated, sorted, and indexed with gatk (). The deduplicated bam files were used for population genetic analysis with ANGSD (Korneliussen et al. 2014) using filtering flags [-uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 -trim 0 -C 50 -baq 1 -minMapQ 20 -minQ 20 -skipTriallelic 1 -SNP_pval 1e-3].
Genotype calling and covariance matrix was done with Angsd (Korneliussen et al. 2014) and its relative ngsTools (Fumagalli et al. 2014), using two to three individuals of the species Atlantic cod, Pacific cod, walleye pollock, Polar cod, Arctic cod, and our hybrid. Positions where all individuals had at least one in coverage depth were included and cut-off where coverage depth of the combined individuals was more than 500 (e.g., if all individuals had ca. 30 in coverage depth).
mtDNA sequencing
We leveraged previously published cytochrome b and mitogenome sequences for Boreogadus (Wilson et al. 2017, 2019, 2020). For the 12 Artogadus samples, we amplified either an 818 or 1,266 bp fragment of cytochrome b using primers and protocols listed in Wilson et al. (2017 and 2019). All new (GenBank accession numbers: OK143300-OK143312) and publicly available sequences (see Wilson et al. 2018 for GenBank accession numbers) were trimmed to a common length of 707 bp.
Levels of divergence
Pair-wise levels of divergence (ϕST) and nucleotide diversity (π) for each ddRAD-seq loci and overall were calculated using the custom script out2phistA.py (available at https://github.com/BU-RAD-seq/Out-Conversions). For mtDNA cytb, we calculated the number of fixed differences and genetic distance (Dxy) between species and nucleotide diversity within species in DnaSP v5.1 (Librado and Rozas 2009).
Genetic admixture
We used three approaches to assess the level of genetic admixture and infer the presence of gene flow between species. First, we performed a Principal Component Analysis (PCA) using allelic states extracted with a custom Python script (out2structureA.py; available at https://github.com/BU-RAD-seq/Out-Conversions) implemented in the adegenet R package (i.e., “dudi.pca”; Dray and Dufour 2007; Jombart 2008). For PCAs, we plotted the first two principal components. Second, we obtained maximum likelihood estimates of population assignments across individuals using ADMIXTURE v.1.3 (Alexander et al. 2009; Alexander and Lange 2011). We used biallelic SNPs, with singletons (i.e., rare SNPs observed in only one individual) excluded and without a priori assignment of individuals to populations or species. SNPs were first formatted for analyses using plink (Purcell et al. 2007), following steps outlined in Alexander et al. (2012). ADMIXTURE analysis was run with a ten-fold cross-validation (CV), and a quasi-Newton algorithm employed to accelerate convergence (Zhou et al. 2011). To limit any possible stochastic effects from single analyses, we ran 100 iterations for each value of K (number of groups; from K of 1–10). Each analysis used a block relaxation algorithm for point estimation and terminated once the change (i.e., delta) in the log-likelihood of the point estimations increased by <0.0001. The optimum K was based on the average of CV errors across the 100 analyses per K; however, additional Ks were analyzed for further population structure resolution. We then used the program CLUMPP v.1.1 (Jakobsson and Rosenberg 2007) to determine the robustness of the assignments of individuals to populations at each K. The R program PopHelper (Francis 2017) was first used to efficiently convert ADMIXTURE outputs into CLUMPP input files at each K. In CLUMPP, we employed the Large Greedy algorithm and 1000 random permutations. Final admixture proportions for each K and per sample assignment probabilities (Q estimates; the log likelihood of group assignment) were based on CLUMPP analyses of all 100 replicates per K.
Third, we used the program fineRADstructure (Malinsky et al. 2018) to infer population structure via shared ancestry. FineRADstructure differs from ADMIXTURE in that it focuses on the most recent coalescent events to provide information on recent relatedness, which is informative in situations of contemporary gene flow between species. Due to the vast differences in sample sizes between species, we ran fineRADstructure using the entire data set as well as randomly selected 15 Boreogadus, including sample ID 20 015 632, which was designated as a potential hybrid (see results below). This was repeated with three sets of 15 randomly chosen Boreogadus samples. Samples were assigned to populations using 1000 000 iterations sampled every 1000 steps with a burn-in of 100 000. The tree was built using 1000 000 iterations and output visualized using the R scripts fineradstructureplot.r and finestructurelibrary.r (available at https://github.com/millanek/fineRADstructure; Malinsky et al. 2018).
Testing for generational hybrids
To better assess whether small interspecific admixture assignments reflect hybridization or retained ancestral variation, we employed methods outlined in Lavretsky et al. (2016) to simulate expected assignment probabilities for first-generation hybrids (F1) and nine generations of backcrosses (F2–F10) into either Arctogadus or Boreogadus for all ddRAD-seq markers. In short, a total of ten F1 hybrids were first generated by randomly sampling an allele from Arctogadus and Boreogadus gene pool across biallelic SNP positions—each position was randomly sampled based on a probability proportional to the allelic frequency in each respective gene pool. Five hybrids were then backcrossed to either the Arctogadus or Boreagadus for nine generations. To limit potential biases in simulations, hybrid indices were reconstructed using only individuals with ADMIXTURE-based probabilities of ≥95% assignment to either Arctogadus or Boreogadus. We ran a total of ten independent simulations, with data subsequently inputted into ADMIXTURE to estimate assignment probabilities for a K of 2. At each K, 25 iterations were run per simulation for a total of 250 ADMIXTURE outputs generated per K, which were then combined and converted in PopHelper (Francis 2017) into CLUMPP input files. We employed the Large Greedy algorithm and 1000 random permutations with final admixture proportions for each K and per sample assignment probabilities based on CLUMPP analyses of all 250 replicates per K. Per generation expected assignment probabilities were based on the average of either all ten (F1) or each of the five (F2–F10) backcrosses, along with each lower and upper limit.
Using the whole-genome data set, a PCA was done as described with ngsTools (Fumagalli et al. 2013) and using PCAngsd (Meisner and Albrechtsen 2018) to illustrate genomic affinity to other parental species. In addition, genetic distance was calculated between the hybrid individual and potential parental species using ngsDist (Vieira et al. 2015).
Results
We obtained 138 598 600 raw sequencing reads with a maximum 150 bp length using single-end sequencing on an Illumina HiSeq 2500. The number of reads per individual ranged from 484 801 to 3 325 291 with an average of 989 990. To select loci for analysis, we first retained loci with at least median depth of 10 within only the Arctogadus data set which resulted in 3332 loci. From that set of 3332 loci, we further filtered the data set by selecting loci with a depth of at least 10 and <15% missing genotypes per locus within the Boreogadus sample set. This resulted in a final data set of 219 monomorphic and 2792 polymorphic loci (3011 total) used for downstream analyses, of which each Arctogadus sample had at least 90% good (i.e., met all filtering thresholds) genotypes (range 90%–94%: average depth 58) and Boreogadus had on average 92% portion of good genotypes (range 75%–95%; average depth 56).
Species divergence—mtDNA
We found 46 fixed bp differences (6.5% of total 707) in the mtDNA cytochrome b gene between Boreogadus and Arctogadus, with an average of 0.074 nucleotide substitutions per site per population (Dxy; see Fig. 2). Of those, 25 polymorphic sites within Boreogadus and three within Arctogadus were monomorphic in the other species. Nucleotide diversity was higher in Boreogadus (0.0036) than in Arctogadus (0.0012).
[Image omitted: See PDF]
Species divergence—ddRAD loci
Overall level of divergence was high (overall ϕST = 0.902; range 0.0–1.0, Table S1) across all linkage groups (LG01–LG23; Fig. 3). Of the 3011 loci analyzed, at least 2219 had a ϕST > 0.1, 1829 loci had a level of divergence nearing fixation (>0.75; Fig. S1), and 401 loci essentially at fixation (between 0.99 and 1.0). This overall pattern of divergence was also reflected when loci were separated by linkage groups based on Atlantic cod genome (LG01–23; Fig. 3). Unlike mtDNA, overall nucleotide diversity within each species was similar: 0.00 316 (range: 0.0–0.1592) for Boreogadus and 0.00 292 (range: 0.0–0.2168; Fig. S1 and Table S1) for Arctogadus. However, 1721 loci were monomorphic within Arctogadus samples, compared to 422 within Boreogadus; 211 loci were monomorphic in both.
[Image omitted: See PDF]
Genetic admixture
An interspecific mtDNA haplotype found in one individual (sample 20015 632) was recovered for one Boreogadus sampled from the Canadian Beaufort Sea (Arctogadus mtDNA haplotype and Boreogadus ddRAD genotype; Fig. 2). Based on PCA and ADMIXTURE for ddRAD loci, all individuals were assigned to their respective species with >>99.9% posterior probabilities (Fig. 2). However, the posterior probability for sample 20 015 632 was slightly lower at 99.5%. The cluster coancestry matrix and cladogram from the fineRADstructure analysis confirmed the ADMIXTURE results; the two species formed distinct clusters. Additionally, we observed more shared coancestry from Arctogadus into Boreogadus than vice versa. Interestingly, sample 20 015632 as well as sample 117 268 in one subset from the Bering Sea had slightly higher levels of shared coancestry in both directions, as donor and recipient, than other samples, potentially indicating late-generation hybridization event. All three subsampled data sets of randomly selected individuals showed similar patterns.
Hybrid identification
Simulations identify an exponential decay in the number of generational backcrosses, and specifically, we can distinguish between F1 through F6 generations, with ≥F6 being genetically indistinguishable from the parental species they backcross into (Fig. 2). In fact, hybrid indices follow expected proportions for generational backcrosses (Table 1). Given simulated results, we are unable to positively detect hybrid or generational backcross individuals among Arctogadus samples (but see above), and only a single Boreogadus (20 015 632), which we categorize as a late-generational backcross into Boreogadus (F5 or F6 generation; Figs. 2 and 4).
[Image omitted: See PDF][Image omitted: See PDF]
Genetic distance based on whole-genome sequencing between hybrid sample 20015632 and potential parental species ranged from 0.137 (Arctogadus) to 0.250 (Greenland Boreogadus; Fig. S2). Despite the closer genetic distance to Arctogadus, a PCA showed sample 20015632 does not cluster with either Boreogadus or Arctogadus (Fig. 4), further indicating that the sample is likely not a “pure” individual of either species or any other gadid. However, Boreogadus could not be confirmed as a major parental species with whole-genome sequencing as it was with the ddRAD sequences.
Discussion
Although they are nearly morphologically identical during certain life stages and partially sympatric, Beaufort Sea Boreogadus and Arctogadus possess distinct nuclear genotypes. These two species are thought to have diverged between 5 and 10 million years ago (Malmstrøm et al. 2016; Baalsrud et al. 2017), which is reflected by the substantial number of loci recovered (>50%) nearing fixation (ϕST > 0.750) and skewed ϕST distribution with a more pronounced tail of higher ϕST values as well as loci with no differentiation (see Fig. 3). Both of these characteristics are indicative of species that have likely undergone allopatric speciation following established reproductive isolation with little or no ongoing geneflow (Via 2001; Savolainen et al. 2006). Therefore, and although our sample size for Arctogadus is low, the absence of early generational hybrids (F1 or F2) is not surprising. The presence of the single-potential late-generational (F5 or F6) hybrid suggests some F1 hybrids could survive to maturity and successfully reproduce. However, it is possible that the genomic signature of mixed ancestry in this one sample could be due to a past introgression event as historical hybridization has been shown to be important in the evolution of gadids (Árnason and Halldórsdóttir 2019). We also showed that introgression or gene flow may be asymmetric from Arctogadus into Boreogadus, perhaps resulting from numerous factors such as differences in mating behavior and timing (Konkle and Philipp 1992), environmental variables influencing fitness (Ostberg et al. 2004; Sotola et al. 2019), or relative abundance of parental species (Lepais et al. 2009).
The single hybrid observed was sampled near the Mackenzie River Delta, second largest Arctic delta, which discharges freshwater that drives productivity on the Beaufort Shelf (Emmerton et al. 2008). Our prior population genetics analyses suggest that this area is an intermixing zone between Boreogadus populations (Arctic and Pacific lineages), and at the regional level, intermixing appears restricted to this area (Wilson et al. 2019). The dynamic and heterogeneous environment in this area, such as the large influx of freshwater discharge, and therefore warmer temperatures in this region, may increase the potential for hybridization to occur. Bouchard and Fortier (2011) have shown that the hatching dates of Boreogadus are correlated with freshwater input with earlier hatching and prolonged hatching season associated with large river discharge, such as the Mackenzie River, relative to other areas, such as the Canadian Archipelago. Although the spawning season for Arctogadus remains uncertain throughout much of the species’ range (Jordan et al. 2003), there is significant overlap in hatching date between the species within the Mackenzie River region, with Boreogadus showing a more flexible and extended hatching season than Arctogadus (Bouchard et al. 2016). If this hatching date overlap equates to overlap in spawning season, this could provide an avenue for greater frequency of hybridization in regions with more environmental heterogeneity that allow for a more flexible spawning/hatching season. Alternatively, although environmental heterogeneity may allow for a greater degree of cross-species interaction, environmentally dependent selection against hybrids in other regions may constrain gene flow as observed between Arctic char (Salvelinus alpinus) and Dolly varden (S. malma; May-McNally et al. 2015). Although larval niches are similar, parental ecological niches of Boreogadus and Arctogadus are quite different (Christiansen et al. 2012); this suggests the presumed low density of hybrids observed could be due to that intermediate morphology performs poorly in parental niches (e.g., Hatfield and Schulter 1999). This may explain the discrepancy between the whole-genome sequencing and ddRAD results on which species are the major parental species. Ecological selection on hybrids is predicted to result in a pattern of biased ancestry, in particularly around functionally important genomic regions (Moran et al. 2021). The use of methylation-sensitive enzymes (MspI in this study) has been shown in plants to recover a high proportion of SNPs located in or close to coding regions of genes (Pootakham et al. 2016). Therefore, our ddRAD results could be exhibiting a signature of the removal of the minor parental species (Arctogadus) and retention of the major parental species (Boreogadus) genomic legacy at functionally important genes as we uncovered a high proportion of loci near fixation. Further, the rate of removal of foreign ancestry is not consistent across the neutral and coding regions of the genome (Moran et al. 2021); therefore, this could also explain why the whole-genome data still retain a higher proportion of the genomic legacy of Arctogadus illustrated by the smaller genetic distance with the hybrid individual. If the whole-genome data contain a higher proportion of noncoding regions, those regions would be purged at a slower rate. A survey of hybridization frequency in other areas of sympatry with no or little river discharge would be needed to test for ecological selection against hybrids and more samples with mixed ancestry would be needed to determine variability in the rate of purging one parental legacy across the genome.
Skewed proportions of parental species could influence the evolutionary outcome of hybridization events (i.e., “Hubb’s effect”; Hubbs 1955). Boreogadus is one of the most abundant fishes in the Arctic with population size estimates ranging in the millions (Rand and Logerwell 2011; De Robertis et al. 2017), while Arctogadus is thought to be in much lower abundance where species co-occur (Cobb et al. 2008, ∼8%–9% of catch Bouchard et al. 2016). Furthermore, there are few reports of Arctogadus west of the Mackenzie River Delta into the Alaska Beaufort Sea (Thorsteinson and Love 2016; Wilson et al. 2019), suggesting this area represents the western edge the species’ range within the Beaufort Sea. If the ratio of population size is skewed as suggested, then hybrids would be expected to mate more frequently with the locally abundant species (Boreogadus in this case), resulting in asymmetric introgression as observed in other fishes (e.g., pufferfish Takifugu; Takahashi et al. 2017). Indeed, our fineRADstructure analysis, which shows a more contemporary genetic structure signature, uncovered a signature of asymmetric geneflow from Arctogadus into the more common Boreogadus. A signature of asymmetric introgression can also explain differences in within-species genetic diversity. Although overall genetic diversity is similar between species, approximately 57% of loci were not variable within the less common Arctogadus. This difference could be due to the low sample of Arctogadus; however, the lower level has been noted in other studies using mitochondrial DNA (Pálsson et al. 2009; Wilson et al. 2020).
Conclusions
Seascape changes, whether naturally or anthropogenically induced, can influence species demography and dispersal patterns, which can also change the rates of gene flow across species. As water temperatures increase in the Arctic, fish communities are already being altered (borealization; Mueter and Litzow 2008; Frainer et al. 2017), which will significantly alter community dynamics. Hybridization in fish has been found to be associated with a decrease in habitat availability, predicted to occur for both Boreogadus and Arctogadus along with other population-limiting environmental factors (e.g., temperature; Marie et al. 2012). Thus, to better understand the mechanisms at both genomic and ecological scales and to inform conservation decision-making to preserve biodiversity, it is increasingly important to identify the degree to which introgressive hybridization and reproductive isolation are occurring. This study is the first step to assessing hybridization between two sympatric cold-water specialists. We have identified 401 loci essentially at fixation (Table S1), and similarly fixed loci relative to other cod species (e.g., walleye pollock), which can be used to monitor the density of hybrids. In addition, these loci can be used to assess spatial and temporal distribution changes at different life stages that may cause an increase in species interactions as shifts have already been reported for Boreogadus (Huserbråten et al. 2019) as well as provide more information on the lesser-known species Arctogadus. A more broad-scale genetic analysis would be required in the Canadian Arctic, where Arctogadus is more abundant and opportunities for hybridization may greater, to determine the extent of potential hybridization and introgression and its potential role in species diversification and future adaptation to changing oceanic conditions in these Arctic specialists.
Acknowledgements
The authors thank B. Norcross, L. Edenfield (University of Alaska Fairbanks), and University of Washington Burke Museum Fish Collection for providing samples for this study, and G. DeGange (USGS) for laboratory support. Canadian samples were collected by Fisheries and Oceans Canada as part of the Beaufort Regional Environmental Assessment–Marine Fishes Project. Study collaboration and funding were provided by the Bureau of Ocean Energy Management (BOEM), Environmental Studies Program, Anchorage, AK under Agreement Number M14PG0008. This research used resources provided by the Core Science Analytics, Synthesis, and Libraries (CSASL) Advanced Research Computing (ARC) group at the U.S. Geological Survey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Alexander D.H., Lange K. 2011. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinf. 12: 246.
Alexander D.H., Novembre J., Lange K. 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19: 1655–1664.
Alexander D.H., Novembre J., Lange K. 2012. Admixture 1.22 Software Manual.
Árnason E., Halldórsdóttir K. 2019. Codweb: whole genome sequencing uncovers extensive reticulations fueling adaptation among Atlantic, Arctic, and Pacific gadids. Sci. Adv. 5: eaat8788.
Baalsrud H.T., Voje K.L., Tørresen O.K., Solbakken M.O., Matschiner M., Malmstrøm M., 2017. Evolution of hemoglobin genes in codfishes influenced by ocean depth. Sci. Rep. 7: 7956.
Barton N.H., Hewitt G.M. 1989. Adaptation, speciation and hybrid zones. Nature, 341: 497–502.
Bouchard C., Fourtier L. 2011. Circum-Arctic comparison of the hatching season of polar cod Boreogadus saida: a test of the freshwater winter refuge hypothesis. Prog. Oceanogr. 90: 105–116.
Bouchard C., Geoffroy M., LeBlanc M., Majewski A., Gauthier S., Walkusz W., 2017. Climate warming enhances polar cod recruitment, at least transiently. Prog. Oceanogr. 156: 121–129.
Bouchard C., Mollard S., Suzuki K., Robert D., Fortier L. 2016. Contrasting the early life histories of sympatric arctic gadids Boreogadus saida and Arctogadus glacialis in the Canadian Beaufort Sea. Polar Biol. 39: 1005–1022.
Bradbury I.R., Bowman S., Borza T., Snelgrove P.V.R., Hutchings J.A., Berg P.R., 2014. Long distance linkage disequilibrium and limited hybridization suggest cryptic speciation in Atlantic cod. PLoS One, 9:e106380.
Campton D.E. 1987. Natural hybridization and introgression in fishes. In: Population genetics and fishery management, N. Ryman, F. Utter(Eds.). University of Washington: Seattle, WA.
Canestrelli D., Bisconti R., Chiocchio A., Maiorano L., Zampiglia M., Nascetti G. 2017. Climate change promotes hybridization between deeply divergent species. PeerJ. 5 : e3072.
Christiansen J.S., Hop H., Nilssen E.M., Joensen J. 2012. Trophic ecology of sympatric arctic gadoids, Arctogadus glacialis (Peters, 1872) and Boreogadus saida (Lepechin, 1774), in NE greenland. Polar Biol. 35: 1247–1257.
Chunco A.J. 2014. Hybridization in a warmer world. Ecol. Evol. 4: 2019–2031.
Cobb D., Fast H., Papst M.H., Rosenberg D., Rutherford R., Sareault J.E. 2008. Beaufort Sea large ocean management area: ecosystem overview and assessment report. Can. Tech. Rep. Fish. Aquat. Sci. 2780: ii–ix + 188 p.
Craig P.C., Griffiths W.B., Haldorson L., McElderry H. 1982. Ecological studies of arctic cod (Boreogadus saida) in Beaufort Sea coastal waters, Alaska, Can. J. Fish. Aquat. Sci. 39: 395–406.
Cronin T.M., Cronin M.A. 2015. Biological response to climate change in the arctic ocean: the view from the past. Arktos, 1: 4.
DaCosta J.M., Sorenson M.D. 2014. Amplification biases and consistent recovery of loci in a double-digest RAD-seq protocol. PLoS One, 9: e106713.
De Robertis A., Taylor K., Wilson C.D., Farley E.V. 2017. Abundance and distribution of arctic cod (Boreogadus saida) and other pelagic fishes over the U.S. continental shelf of the northern Bering and Chukchi seas. Deep-Sea Res. II 135: 51–65.
Dowling T.E., Markle D.F., Tranah G.J., Carson E.W., Wagman D.W., May B.P. 2016. Introgressive hybridization and the evolution of lake-adapted catostomid fishes. PLoS One 11: e0149884.
Dray S., Dufour A-B 2007. The ade4 package: implementing the duality diagram for ecologists. J. Stat. Softw. 22: 1–20.
Edgar R.C. 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32: 1792–1797.
Edgar R.C. 2010. Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 26: 2460–2461.
Edmands S. 2007. Between a rock and a hard place: evaluating the relative risks of inbreeding and outbreeding for conservation and management. Mol. Ecol. 16: 463–475.
Emmerton C.A., Lesack L.F.W., Vincent W.F. 2008. Mackenzie river nutrient delivery to the arctic ocean and effects of the mackenzie delta during open water conditions. Global Biogeochem. Cycles, 22: GB1024.
Farkas T.E., Hendry A.P., Nosil P., Beckerman A.P. 2015. How maladaptation can structure biodiversity: eco-evolutionary island biogeography. Trends Ecol. Evol. 30: 154–160.
Finley K., Gibb E. 1982. Summer diet of the narwhal (Monodon monoceros) in pond inlet, northern Baffin Island. Can. J. Zool. 60: 3353–3363.
Frainer A., Primicerio R., Kortsch S., Aune M., Dolgov A.V., Fossheim M., Aschan M.M. 2017. Climate-driven changes in functional biogeography of arctic marine fish communities. Proc. Natl. Acad. Sci. U.S.A. 114: 12202–12207.
Francis R.M. 2017. Pophelper: an r package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 17: 27–32.
Frisch A.J., van Herwerden L. 2006. Field and experimental studies of hybridization between coral trouts, Plectropomus leopardus and Plectropomus maculatus (Serranidae), on the Great Barrier Reef, Australia. J. Fish Biol. 68L: 1013–1025.
Fumagalli M., Vieira F.G., Korneliussen T.S., Linderoth T., Huerta-Sánchez E., Albrechtsen A., Nielsen R. 2013.Quantifying population genetic differentiation from next-generation sequencing data. Genetics, 195: 979–992.
Fumagalli M., Vieira F.G., Linderoth T., Nielsen R. 2014. ngsTools: methods for population genetics analyses from next-generation sequencing data. Bioinformatics, 30: 1486–1487.
Garroway C.L., Bowman J., Cascaden T.J., Holloway G.L., Mahan C.G., Malcolm J.R., 2010. Climate change induced hybridization in flying squirrels. Global Change Biol. 16: 113–121.
Gaspar J.M. 2018. NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors. BMC Bioinf. 19: 536.
Hanemaaijer M.J., Collier T.C., Chang A., Shott C.C., Houston P.D., Schmidt H., 2018. The fate of genes that cross species boundaries after a major hybridization event in a natural mosquito population. Mol. Ecol. 27: 4978–4990.
Harwood L.A., Smith T.G., George J.C., Sandstrom S.J., Walkusz W., Divoky G.J. 2015. Change in the Beaufort Sea ecosystem: diverging trends in body condition and/or production in five marine vertebrate species. Prog. Oceanogr. 136: 263–273.
Hatfield T., Schluter D. 1999. Ecological speciation in sticklebacks: environment-dependent hybrid fitness. Evolution, 53: 866–873.
Hubbs C. L. 1955. Hybridization between fish species in nature. Syst. Zool. 4: 1–20.
Huberbråten M.B.O., Eriksen E., Gjæsæter, Vikebø F. 2019. Polar cod in jeopardy under the retreating arctic sea ice. Commun. Biol. 2: 407.
Jakobsson M., Rosenberg N.A. 2007. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23: 1801–1806.
Jombart T. 2008. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24: 1403–1405.
Jordan A.D., Møller P.R., Nielsen J.G. 2003. Revision of the Arctic cod genus Arctogadus. J. Fish Biol. 62: 1339–1352.
Konkle B. R, Philipp D. P. 1992. Asymmetric hybridization between two species of sunfishes (Lepomis:Centrarchidae). Mol. Ecol. 1: 215–222.
Korneliussen T. S., Albrechtsen A., Nielsen R. 2014. ANGSD: analysis of next generation sequencing data. BMC Bioinf. 15: 356.
Laurel B.J., Spencer M., Iseri P., Copeman L.A. 2016. Temperature-dependent growth and behavior of juvenile arctic cod (Boreogadus saida) and co-occurring North Pacific gadids. Polar Biol. 39: 1127–1135.
Lavretsky P., Peters J.L., Winker K., Bahn V., Kulikova I., Zhuravlev Y.N., 2016. Becoming pure: identifying generational classes of admixed individuals within lesser and greater scaup populations. Mol. Ecol. 25: 661–674.
Leary R.F., Allendorf F.W., Sage G.K. 1995. Hybridization and introgression between introduced and native fish. Am. Fish. Soc. Symp. 15: 91–101.
LeBlanc M., Geoffroy M., Bouchard C., Gauthier S., Majewski A., Reist J.D., Fortier L. 2020. Pelagic production and the recruitment of juvenile polar cod Boreogadus saida in Canadian arctic seas. Polar Biol. 43: 1043–1054.
Lepais O, Petit R. J, Guichoux E, Lavabre J. E, Alberto F, Kremer A, Gerber S. 2009. Species relative abundance and direction of introgression in oaks. Mol. Ecol. 18: 2228–2242.
Li H., Durbin R. 2009. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 25: 1754–1760.
Librado P., Rozas J. 2009. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452.
Ma Y., Marczewski T., Xue D., Wu Z., Laio R., Sun W., Marczewski J. 2019. Conservation implications of asymmetric introgression and reproductive barriers in a rare primrose species. BMC Plant Biol. 19: 286.
Madsen M. L., Fevolden S.-E., Christiansen J. S. 2009. A single molecular approach to distinguish between two arctic gadoid fishes Arctogadus glacialis (Peters, 1874) and Boreogadus saida (Lepechin, 1774). Polar Biol. 32: 937–939.
Majewski A., Atchinson S., MacPhee S., Eert J., Niemi A., Michel C., Reist J.D. 2017. Marine fish community structure and habitat associations on the Canadian Beaufort shelf and slope. Deep Sea Res. Part I 121: 169–182.
Malinksy M., Trucchi E., Lawson D.J., Falush D. 2018. RADpainter and fineRADstructure: population inference from RADseq data. Mol. Biol. Evol. 35: 1284–1290.
Mallet J. 2005. Hybridization as an invasion of the genome. Trends Ecol. Evol. 20: 229–237.
Malmstrøm M., Matschiner M., Tørresen O., Star B., Snipen L.G., Hansen T.F., 2016. Evolution of the immune system influences speciation rates in teleost fishes. Nat. Genet. 48: 1204–1210.
Marie A., van Herwerden L., Choat J., Hobbs J.-P. A. 2007. Hybridization of reef fishes at the indo-pacific biogeographic barrier: a case study. Coral Reefs, 26: 841–850.
Marie A.D., Bernatchez L., Garant D. 2012. Environmental factors correlate with hybridization in stocked brook charr (Salvelinus fontinalis). Can. J. Fish. Aquat.Sci. 69: 884–893.
Marsh J.M., Mueter F.J. 2020. Influence of temperature, predators, and competitors of polar cod (Boreogadus saida) at the southern margin of their distribution. Polar Biol. 43: 995–1014.
May-McNally S., Quinn T.P., Taylor E.B. 2015. Low levels of hybridization between sympatric arctic char (Salvelinus alpinus) and dolly varden char (Salvelinus malma) highlights their genetic distinctiveness and ecological separation. Ecol. Evol. 5: 3031–3045.
Meisner J., Albrechtsen A. 2018. Inferring population structure and admixture proportions in low depth NGS data. Genetics, 210: 719–731.
Montanari S.R., van Herwerden L., Pratchett M.S., Hobbs J-P. A., Fugedi A. 2012. Reef fish hybridization: lessons learnt from butterflyfishes (genus Chaetodon). Ecol. Evol. 12: 310–328.
Moran B.M., Payne C., Langdon Q., Powell D.L., Brandvain Y., Schumer M. 2021. The genomic consequences of hybridization. eLife, 10: e69016.
Mueter F.J., Litzow M.A. 2008. Sea ice retreat alters the biogeography of the Bering Sea continental shelf. Ecol. Appl. 18: 309–320.
Mueter F.J., Nahrgang J., John Nelson R., Berge J. 2016. The ecology of gadid fishes in the circumpolar arctic with a special emphasis on the polar cod (Boreogadus saida). Polar Biol. 39: 961–967.
Ostberg C. O., Slatton S.L., Rodriguez R.J. 2004. Spatial partitioning and asymmetric hybridization among sympatric coastal steelhead trout (Oncorhynchus mykiss irideus), coastal cutthroat trout (O. clarki clarki) and interspecific hybrids. Mol. Ecol. 13: 2773–2788.
Pálsson S., Källman T., Paulsen J., Árnason E. 2009. An assessment of mtDNA variation in arctic gadoids. Polar Biol. 32: 471–479.
Pardo-Diaz C., Salazar C., Baxter S.W., Merot C., Figueiredo-Ready W., Joron M., 2012. Adaptive introgression across species boundaries in Heliconius butterflies. PLos Genet. 8: e1002752.
Pootakham W., Sonthirod C., Naktang C., Jomchai N., Sangsraku D., Tangphatsornruang S. 2016. Effects of methylation-sensitive enzymes on the enrichment of genic SNPs and the degree of genome complexity reduction in a two-enzyme genotyping-by-sequencing (GBS) approach: a case study in oil palm (Elaeis guineensis). Mol. Breed. 36: 154.
Potts W.M., Henqiques R., Santos C.V., Munnik K., Ansorge I., Dufois F., 2014. Ocean warming, a rapid distributional shift, and the hybridization of a coastal fish species. Global Change Biol. 20: 2765–2777.
Purcell S., Neale B., Todd-Brown K., Thomas L., Ferreira M.A.R., Bender D., 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81: 559–575.
Rand K.M., Logerwell E.A. 2011. The first demersal trawl survey of benthic fish and invertebrates in the Beaufort Sea since the late 1970s. Polar Biol. 34: 475–488.
Savolainen V., Anstett M.C., Lexer C., Hutton I., Clarkson J.J., Norup M.V., 2006. Sympatric speciation in palms on an oceanic island. Nature, 441: 210–213.
Sotola V.A., Ruppel D.S., Bonner T.H., Nice C.C., Martin N.H. 2019. Asymmetric introgression between fishes in the Red River basin of Texas is associated with variation in water quality. Ecol. Evol. 9: 2083–2095.
Süfke L., Piepenburg D., von Dorrien C.F. 1998. Body size, sex ratio and diet composition of Arctogadus glacialis (Peters, 1874) (Pisces:Gadidae) in the northeast water polynya (Greenland). Polar Biol. 20: 357–363.
Takahashi H., Toyoda A., Yamazaki T., Narita S., Mashiko T., Yamazaki Y. 2017. Asymmetric hybridization and introgression between sibling species of the pufferfish Takifugu that have undergone explosive speciation. Mar. Biol. 164: 90.
Teletchea F., Laudet V., Hänni C. 2006. Phylogeny of the Gadidae (sensu Svetovidov, 1948) based on their morphology and two mitochondrial genes. Mol. Phylogenet. Evol. 38: 189–199.
Thorsteinson L.K., Love M.S., eds, 2016. Alaska Arctic marine fish ecology catalog: U.S. Geological Survey Scientific Investigations Report 2016–5038 (OCS Study, BOEM 2016-048), 768p.
Tørresen O.K., Star B., Jentoft S., Reinar W.B., Grove H., Miller J.R., 2017. An improved genome assembly uncovers prolific tandem repeats in atlantic cod. BMC Genomics, 18: 95.
Via S. 2001. Sympatric speciation in animals: the ugly duckling grows up. Trends Ecol. Evol. 16: 381–390.
Vieira F.G., Lassalle F., Korneliussen T.S., Fumagalli M. 2015. Improving the estimation of genetic distances from Next-generation sequence data. Biol. J. Linn. Soc. 117: 139–149.
Weist P., Schade F.M., Damerau M., Barth J.M.I., Dierking J., Andre C., 2019. Assessing SNP-markers to study population mixing and ecological adaptation in baltic cod. PLoS One, 14: e0218127.
Welch H.E., Bergmann M.A., Siferd T.D., Martin K.A., Curtis M.F., Crawford R.E., 1992. Energy flow through the marine ecosystem of the Lancaster Sound region, Arctic Canada. Arctic, 45: 343–357.
Wilson R.E, Sage G.K., Sonsthagen S.A., Gravley M.C., Menning D.M., Talbot S.L. 2017. Genomics of Arctic Cod. Anchorage, AK: US Dept. of the Interior, Bureau of Ocean Energy Management, Alaska OCS Region. OCS Study BOEM 2017-066. 92pp.
Wilson R.E., Pierson B.J., Sage G.K., Sonsthagen S.A., Gravley M.C., Menning D.M., Talbot S.L. 2018. Genetic data from arctic, polar, and saffron cod and walleye pollock, Alaska and Canada, 2011–2017: U.S. Geological Survey data release. Available from https://doi.org/10.5066/F7CF9P23.
Wilson R.E., Sage G.K., Wedemeyer K., Sonsthagen S.A., Menning D.M., Gravley M.C., 2019. Micro-geographic population genetic structure within arctic cod (Boreogadus saida) in Beaufort Sea of Alaska. ICES J. Mar. Sci. 76: 1713–1721.
Wilson R.E., Sonsthagen S.A., Sme N., Gharrett A.J., Majewski A.R., Wedemeyer K., 2020. Mitochondrial genome diversity and population mitogenomics of polar cod (Boreogadus saida) and arctic dwelling gadoids. Polar Biol. 43: 979–994.
Yaakub S., Bellwood D., van Herwerden L. 2007. A rare hybridization event in two common caribbean wrasses (genus Halichoeres; family Labridae). Coral Reefs, 26: 597–602.
Zhou H., Alexander D., Lange K. 2011. A quasi-Newton acceleration for high-dimensional optimization algorithms. Stat. Comput. 21: 261–273.
Robert E. Wilson https://orcid.org/0000-0003-1800-0183 [email protected]
School of Natural Resources, University of Nebraska-Lincoln, Lincoln, NE 68583, USA
University of Nebraska State Museum, University of Nebraska-Lincoln, Lincoln, NE 68588, USA
Author Contributions: Conceptualization, Data curation, Methodology, Formal analysis, Investigation, Writing – original draft, and Writing – review & editing.
Sarah A. Sonsthagen
U. S. Geological Survey–Nebraska Cooperative Fish and Wildlife Research Unit, School of Natural Resources, University of Nebraska-Lincoln, Lincoln, NE 68583, USA
Author Contributions: Conceptualization, Investigation, Formal analysis, and Writing – review & editing.
Philip Lavretsky
Department of Biological Sciences, University of Texas at El Paso, El Paso, TX 79668, USA
Author Contributions: Formal analysis and Writing – review & editing.
Andrew Majewski
Freshwater Institute, Fisheries and Oceans Canada, Winnipeg, MB R3T 2N6, Canada
Author Contributions: Resources and Writing – review & editing.
Einar Árnason
Institute of Life and Environmental Sciences, University of Iceland, Reykjavik, IS102, Iceland
Author Contributions: Investigation, Formal analysis, and Writing – review & editing.
Katrín Halldórsdóttir
Institute of Life and Environmental Sciences, University of Iceland, Reykjavik, IS102, Iceland
Author Contributions: Investigation, Formal analysis, and Writing – review & editing.
Axel W. Einarsson
Institute of Life and Environmental Sciences, University of Iceland, Reykjavik, IS102, Iceland
Author Contributions: Investigation, Formal analysis, and Writing – review & editing.
Kate Wedemeyer
Bureau of Ocean Energy Management, Anchorage, AK 99503, USA
Author Contributions: Funding acquisition, Project administration, and Writing – review & editing.
Kate Wedermeyer is retired from the Bureau of Ocean Energy Management, Anchorage, Alaska, USA
Sandra L. Talbot
Far Northwestern Institute of Art and Science, 427 D Street, Anchorage, AK 99501, USA
Author Contributions: Conceptualization, Funding acquisition, and Writing – review & editing.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2022. This work is published under https://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
As marine ecosystems respond to climate change and other stressors, it is necessary to evaluate current and past hybridization events to gain insight on the outcomes and drivers of such events. Ancestral introgression within the gadids has been suggested to allow cod to inhabit a variety of habitats. Little attention has been given to contemporary hybridization, especially within cold-water-adapted cod (Boreogadus saida Lepechin, 1774 and Arctogadus glacialis Peters, 1872). We used whole-genome, restriction-site associated, and mitochondrial sequence data to explore the degree and direction of hybridization between these species where previous hybridization had not been reported. Although nearly identical morphologically at certain life stages, we detected very distinct nuclear and mitochondrial lineages. We detected one potential hybrid with a Arctogadus mitochondrial haplotype and Boreogadus nuclear genotype, but no early generational hybrids. The presence of a late generation hybrid suggests that at least some hybrids survive to maturity and reproduce. However, a historical introgression event could not be excluded. Contemporary gene flow appears asymmetrical from Arctogadus into Boreogadus, which may be due to overlap in timing of spawning, environmental heterogeneity, or differences in population size. This study provides important baseline information for the degree of potential hybridization between these species within Alaska marine environments.
