Extensive secondary contact among three glacial lineages of Arctic Char (Salvelinus alpinus) in Labrador and Newfoundland

Abstract Aim The Pleistocene glaciation event prompted the allopatric divergence of multiple glacial lineages of Arctic char (Salvelinus alpinus), some of which have come into secondary contact upon their recolonization of the Holarctic. While three glacial lineages (Arctic, Atlantic, and Acadian) are known to have recolonized the western Atlantic, the degree of overlap of these three lineages is largely unknown. We sought to determine the distribution of these three glacial lineages in Labrador and Newfoundland at a fine spatial scale to assess their potential for introgression and their relative contribution to local fisheries. Location Labrador and Newfoundland, Canada. Methods We sequenced a portion of the D‐loop region in over 1,000 Arctic char (S. alpinus) samples from 67 locations across Labrador and Newfoundland. Results Within Labrador, the Arctic and Atlantic lineages were widespread. Two locations (one landlocked and one with access to the sea) also contained individuals of the Acadian lineage, constituting the first record of this lineage in Labrador. Atlantic and Acadian lineage individuals were found in both eastern and western Newfoundland. Multiple sampling locations in Labrador and Newfoundland contained fish of two or more different glacial lineages, implying their introgression. Glacial lineage did not appear to dictate contemporary genetic divergence between the pale and dark morph of char present in Gander Lake, Newfoundland. Both were predominately of the Atlantic lineage, suggesting the potential for their divergence in sympatry. Main conclusions Our study reveals Labrador and Newfoundland to be a unique junction of three glacial lineages which have likely hybridized extensively in this region.


| INTRODUC TI ON
Glaciation events are a significant driver of evolution, physically isolating species into separate glacial refugia which can undergo allopatric divergence for thousands of years (Fraser, Nikula, Ruzzante, & Waters, 2012;Hewitt, 2000Hewitt, , 2004. During allopatry, populations may experience differential selection and drift resulting in the formation of genetically distinct glacial lineages (Hewitt, 2003;Moore, Bajno, Reist, & Taylor, 2015;Ruzzante et al., 2008). Retreating glaciers allowed access to new environments, sometimes facilitating secondary contact (Hewitt, 2000;Soltis, Morris, McLachlan, Manos, & Soltis, 2006;Swenson & Howard, 2005). Upon secondary contact, glacial lineages may demonstrate: (a) extensive gene flow, leading to complete genomic introgression; (b) complete reproductive isolation and an absence of gene flow; and (c) some intermediate level of gene flow (Hewitt, 1988;Noor, 1999;Schluter, 2001). The degree of hybridization is likely to depend on the accumulated genetic divergence among lineages and on the adaptive quality of each gene (Hewitt, 1988). The amount of genetic divergence accumulated between glacial lineages and the degree of erosion of this divergence in secondary contact zones can significantly influence the contemporary genetic structure of a species Hewitt, 2000Hewitt, , 2004. These areas of secondary contact and hybridization therefore not only inform conservation management but also offer natural experiments for the study of the factors driving speciation (Hewitt, 1988).
Some of the lakes colonized by anadromous char in Labrador subsequently lost their access to the sea resulting in contemporarily landlocked char populations (Scott & Crossman, 1998). Other anadromous char populations (particularly in Newfoundland) lost their anadromous lifestyle and remain lacustrine residents year-round (Scott & Crossman, 1998).
The glacial lineages present in anadromous versus landlocked populations remain largely unknown in this region. Landlocked and anadromous populations might have been founded by different lineages, for example, if one was better adapted to a particular environment or life history. Alternatively, landlocked populations could have been founded only by lineages that were present before access to these lakes was lost. Investigation of which lineages are present in these two types of populations may therefore give an indication of the timing of clonization by different glacial lineages (Moore et al., 2015).
Within-lake genetic structure, previously found in Labrador and Newfoundland char populations, may also be influenced by glacial lineage. Glacial lineage has been suggested as the origin of the substantial genetic divergence observed between a pale and a dark morph documented for Gander Lake in Newfoundland (Gomez-Uchida, Dunphy, O'Connell, & Ruzzante, 2008). Salisbury et al. (2018) alternatively found that the genetic structure in two landlocked lakes and one sea-accessible lake in Labrador were unrelated to glacial lineage.
Anadromous Arctic char populations are economically significant and form the basis of a commercial, recreational, and subsistence fishery in Labrador (DFO, 2001;Dempson, Shears, Furey, & Bloom, 2008). Historic allopatry may be an important underlying influence on genetic structure if char of different glacial lineages contribute to the fishery but remain reproductively isolated.
While it is currently unknown which glacial lineages contribute to the Labrador fishery, this knowledge is potentially critical for its management.
Here, we investigated the consequences of secondary contact of the Arctic, Atlantic, and Acadian glacial lineages across Labrador and Newfoundland at a fine spatial scale. We predicted that the Arctic lineage would be more prevalent in northern populations than the Atlantic lineage based on the hypothesis that Labrador was colonized from the north by the Arctic lineage and from the south by the Atlantic lineage. Hybridization among lineages was expected to be prevalent and evidenced by multiple lineages co-existing in single populations. In Labrador, we anticipated that lineages which colonized more recently would be present only in the sea-accessible but not the landlocked populations. Finally, we hypothesized that the large divergence between the pale and dark morphs within Gander Lake was due to their founding by different glacial lineages. To test these hypotheses, we employed mtDNA to identify the glacial lineage of hundreds of fish across Labrador and Newfoundland.
Collections from eastern Newfoundland originate from two locations containing only freshwater residents: Gander Lake (including samples of the pale and dark morphs described for this lake (Gomez-Uchida et al., 2008)) and Wing Pond.
Samples from Labrador were collected using electrofishing in the rivers (sea-accessible sites) and variably sized standardized nylon monofilament gillnets (1.27-8.89 cm diagonal) at the landlocked and sea-accessible lake sites. Samples were collected from anadromous char populations in the Okak and Voisey regions as well as from the Fraser River, Anaktalik River, and Tikkoatokak River. These populations contribute to the three stock complexes of the commercial Labrador char fishery (Okak, Voisey, Nain) (DFO, 2001). Gander Lake was sampled using Lundgren multimesh gillnets (Hammar & Filipsson, 1985) (bar length from 0.625 cm to 7.5 cm) (see Gomez-Uchida et al., 2008, for more details). The Upper Humber River was sampled using fyke nets and electrofishing (see Gomez-Uchida et al., 2009, for more details). Wing Pond was sampled with gillnets. Fish were weighed, measured for fork length (FL) in mm, and assessed for sex and maturity. Tissue samples (fin or gill) were obtained and immediately stored in 95% ethanol; alternatively, some fin clip samples were stored dry. All samples were collected in collaboration with the Department of Environment and Conservation for Labrador and Newfoundland and/or Parks Canada and in accordance with Dalhousie University's Animal Ethics Guidelines.

| DNA extraction, amplification, and genotyping
Tissue samples were digested at 55°C for approximately eight hours using Proteinase K (Bio Basic Inc., Markham, ON, Canada). DNA was then extracted using a Multiprobe II plus liquid handling system (Perkin Elmer, Waltham, MA, USA) using a glassmilk protocol modified from Elphinstone, Hinten, Anderson and Nock (2003).
The left domain region of the mitochondrial control region was amplified and sequenced following Moore et al. (2015). In brief, the primers Tpro2 (Brunner et al., 2001) and SalpcrR (Power, Power, Reist, & Bajno, 2009) were used to amplify the entire control region using the thermocycler program and PCR outlined in Brown Gladden, Postma Maiers, Carmichael, and Reist (1995). A shorter fragment was amplified using Char3 instead of Tpro2 for a minority of samples which had poor quality as determined from visual inspection of a 1% agarose gel. For all samples, a total of ~500 bp of the left domain was sequenced using Char3 (Power et al., 2009) at MacrogenUSA (Rockville, MD). Each unique haplotype detected was validated by resequencing a representative sample for each haplotype using Tpro2.

| Analyses
Our sequences were trimmed, validated, and aligned using GENEIOUS (10.0.9, Auckland, NZ, www.geneious.com). Using default alignment parameters, our sequences were aligned to a reference haplotype set (control region haplotypes verified by Moore et al. (2015) and Salisbury et al. (2018)), and control region sequences for three other salmonid species present in the region (brook trout (Salvelinus fontinalis), lake trout (Salvelinus namaycush), and Atlantic salmon (Salmo salar); for accession numbers see Supporting Information Table S1). Sequences verified as Arctic char were ascribed to the reference haplotype(s) for which they had 0 basepair differences. Non-char sequences were ascribed to the brook trout, lake trout, or Atlantic salmon haplotype to which they had the minimum number of basepair differences.
A representative forward sequence for each unique haplotype (i.e., those sequences which contained one or more basepair differences from those haplotypes verified by Moore et al. (2015)) was aligned with its reverse complement (sequenced with Tpro2) using a pairwise GENEIOUS alignment and default parameters to create a consensus sequence. The consensus sequences for these unique haplotypes were then aligned with the reference haplotype set using a GENIOUS alignment. A gap penalty of 7 was used, and all other parameters were kept at default values. A maximum-likelihood tree was constructed based on this alignment using the PhyML (Guindon & Gascuel, 2003) plugin in GENEIOUS to compare the phylogenetic relationships among these unique consensus sequences with those haplotypes verified by Moore et al. (2015) and those of an outgroup species (brook trout). The Nearest Neighbour Interchange topology search algorithm and the HKY85 + I + G model were used to calculate 1,000 bootstraps for each node following Moore et al. (2015).
A haplotype map based on all unique haplotypes found in this study along with all haplotypes verified by Moore et al. (2015) was created using PopArt version 1.7 (Leigh & Bryant, 2015). Haplotypes were trimmed to 501 bp, the length for which all haplotypes had no missing base pairs, since PopArt masks missing base pairs. This meant that haplotype ATL04 and a unique haplotype ATL31 were indistinguishable in this analysis since the SNP differentiating these haplotypes lies outside of this 501 bp region. The haplotype map was created using a Median-Joining network (Bandelt, Forster, & Röhl, 1999) with an Epsilon value of 0.
A spatial analysis of molecular variance (SAMOVA 2.0) (Dupanloup, Schneider, & Excoffier, 2002) was employed to detect groups of sampling locations whose F CT were maximally differentiated based on mtDNA sequences. All sequences were aligned using GENEIOUS alignment and default parameters and trimmed to 482 bp to include all relevant SNPs differentiating haplotypes.
Sampling locations with fewer than 10 sequences were excluded from the analysis to minimize the probability of biased groupings due to small sampling size. F CT values were estimated using a simulated annealing optimization process for K = 2-10 groups for all sampling locations and for only the Labrador sampling locations, and for K = 2-4 groups for only the Newfoundland sampling locations. For each K-value, molecular distance was calculated using Tamura and Nei distance between all sampling locations and between only those sampling locations connected using a Delaunay network (Delaunay, 1934) based on the latitude and longitude of each sampling location.
The use of a Delaunay network limits groupings to geographically proximate sampling locations. Simulations were run for 10,000 steps from 100 initial configurations using a missing data value of 1 (such that the entire 482 bp was included in the analysis).
Linear regressions between latitude and the number of lineages present in each location as well as binomial logistic regressions of latitude on the presence or absence (coded as 1 and 0, respectively) of each of the relevant glacial lineages were conducted using R (R Core Team, 2013) for all locations, only Labrador locations, and only Newfoundland locations.

| Species distribution
Samples from 59 locations in Labrador (of which 43 are sea-accessible and 16 are landlocked (Anderson, 1985)) and eight locations in Newfoundland (all containing lacustrine residents) for a total of 67 locations overall were successfully sequenced (Table 1). Five locations in Labrador were excluded from analyses due to poor sequence quality (N = 109 individuals). A further 20 individuals from across the remaining 67 locations were excluded from analyses due to poor sequence quality. A total of N = 1,296 individuals were successfully sequenced across all locations in Labrador and Newfoundland (Table 1). Of these, 1,133 had haplotypes consistent with the Arctic char species. The remaining 163 individuals were identified as brook trout, lake trout, or Atlantic salmon. All locations contained at least one Arctic char haplotype except G06 which only contained Atlantic salmon. In the remaining locations, we sampled between 1 and 48 Arctic char (median = 18.5), which made up between 2% and 100% of the haplotypes in each sample.

| Glacial lineage distribution
The Arctic and Atlantic glacial lineages were ubiquitous across Labrador, and both lineages were present in all 10 drainages Across all sampling locations, the presence of the Atlantic lineage was unrelated with latitude (p ≥ 0.435). However, across Labrador and Newfoundland, the probability of Arctic lineage presence increased with latitude (p ≤ 6.46 × 10 −3 ). Similarly, the presence of the Acadian lineage was inversely related to latitude (p ≤ 1.08 × 10 −4 ).

| Haplotype distribution
The most common haplotype within a lineage coincided with the most common haplotypes reported in Moore et al. (2015) (haplotypes at each location: see Supporting Information Table S2). A total of 86% of Atlantic lineage individuals exhibited haplotype ATL01, while 78% of Acadian lineage individuals exhibited haplotype ACD9.
All Labrador samples of the Acadian lineage had this haplotype.
Lastly, over 99% of Arctic lineage individuals exhibited haplotypes ARC19 or ARC24. These two haplotypes were distinguished by a single SNP outside of the region sequenced using the Char3 primer.
However, the reverse complement of 29 samples from 28 locations and nine drainages with either the ARC19 or ARC24 haplotype was sequenced using Tpro2 and all were found to have the haplotype ARC19. Therefore, unambiguously ARC19 sequences were grouped with sequences that could be either ARC19 or ARC24 when counting the number of haplotypes present in a given site. The ARC19, ATL01, and ACD9 haplotypes were found across the modern distributions of the Arctic, Atlantic, and Acadian lineages, respectively, by Moore et al. (2015).
Other detected haplotypes which had been previously described by Moore et al. (2015) and Salisbury et al. (2018) include ACD11, ARC20, ARC22, ATL19, ATL23, ATL24, and ATL25. A single sample had the ATL19 haplotype, a dark morph char from Gander. This was also the only Atlantic haplotype other than ATL01 detected in Newfoundland. The ATL23, ATL24, and ATL25 haplotypes were only observed in S03 and S04 as described in Salisbury et al. (2018) except for one individual with ATL23 found in S02. This individual may have been washed downstream from the immediately upstream landlocked S03 and S04. This is supported by its identification as a putative migrant from S03 based on GENECLASS2 (Piry et al., 2004) results as reported in Salisbury et al. (2018).
Five samples from three landlocked sites in the Voisey drainage had shortened sequences that prevented their differentiation TA B L E 1 Number of Arctic char (Salvelinus alpinus) samples, glacial lineages, and haplotypes as well as number of brook trout (Salvelinus fontinalis), lake trout (Salvelinus namaycush), and Atlantic salmon (Salmo salar) samples verified by mtDNA sequencing at sampling locations across Labrador and Newfoundland. Accessibility of locations (A for sea-accessible, L for landlocked).

Drainage Watershed
Latitude, longitude

S. alpinus
between ATL01 and ATL04. Since these sites also contained individuals unambiguously identified as ATL01, these shortened sequences were considered to be ATL01 when counting the number of haplotypes present in these lakes. The ATL04 haplotype was also found in 12 sea-accessible sampling locations across seven drainages in Labrador. The reverse complement of nine of these samples from six drainages was sequenced using Tpro2 and all were found to contain a consistent SNP in the consensus sequence, differentiating this haplotype from ATL04. One of these samples was from R01, previously mistakenly identified as ATL04 in Salisbury et al. (2018). Given the consistency of this SNP across samples from multiple drainages, sampling locations, and studies, we denoted this as a new haplotype ATL31. We considered all ATL04 haplotypes and verified ATL31

| Landlocked versus sea-accessible sampling locations
The average number of lineages observed in anadromous sites was higher than in landlocked sites ( Gander Lake and Wing Pond maintain sea access but contain lacustrine residents, and these lakes were therefore categorized as "landlocked" for analyses.

| SAMOVA
SAMOVA results were similar across all samples and when considering Labrador and Newfoundland sampling locations separately.
Results were also similar with and without the use of a Delaunay network to take into account geographic proximity of locations. When considering all sampling locations, F CT was maximized for K = 6 ( Figure 4). However, the difference in F CT between K = 6 and K = 5 was small (i.e., 0.07) and a plot of F CT versus K revealed that F CT leveled off at K = 5 (Supporting Information Figure S1a). Given
Unlike the Atlantic lineage, the Arctic lineage does not appear to have invaded Newfoundland. The absence of the Arctic lineage from Newfoundland may be due to our low sample sizes F I G U R E 2 Maximum-likelihood phylogenetic tree of Arctic char (Salvelinus alpinus) haplotypes of the mtDNA control region. Tree was generated using PhyML (Guindon & Gascuel, 2003) with 1,000 bootstrap replicates. Those bootstrap values greater than 50% are shown on the tree. Haplotypes are color-coordinated by lineage as designated in Moore et al. (2015): blue-Arctic, red-Bering, orange-Siberia, purple-Atlantic, green,-Acadian. New haplotypes identified in this study and Salisbury et al., 2018 are starred and the fewer number of locations sampled. However, our study confirmed the previous observations of the Atlantic and Acadian lineage in eastern Newfoundland (Brunner et al., 2001;Moore et al., 2015). We also demonstrated that the contemporary range of both the Atlantic and Acadian lineages extends to western Newfoundland.
Our extension of the Acadian lineage's contemporary presence into Labrador counters previous suggestions of its relatively conserved range from its putative refugium near the northeastern United States (e.g., Brunner et al., 2001;Esin & Markevich, 2018).
This brings the scale of the contemporary range of the Acadian lineage in line with those observed in other char lineages.

| Evidence for introgression
Our results here suggest extensive secondary contact but also the introgression of the Arctic, Atlantic, and Acadian lineages in Labrador and Newfoundland. Many sampled locations contained at least two glacial lineages suggesting the potential for hybridization among lineages. Furthermore, we found six sea-accessible locations in Nachvak and Saglek fjords (N01-N04, S01-S02) contained both Arctic and Atlantic lineage individuals based on mtDNA, yet no genetic structuring was found within each of these same locations based on 11 microsatellite markers in Salisbury et al. (2018). This suggests that these lineages have fully introgressed.
Hybridization between these lineages may seem surprising given that Arctic lineage is thought to have split off from all other lineages between 716,000 and 1,432,000 years BP based on mtDNA (Moore et al., 2015). Alternatively, Esin and Markevich (2018) estimate the divergence of the Arctic lineage at 400,000-700,000 years BP during the Nebraskan-Kansan cooling. During this time, the Canadian Arctic archipelago (the putative refugium for the Arctic lineage) was separated from a refugium in the Bering Sea (Esin & Markevich, 2018). Many species have demonstrated reproductive isolation between different glacial lineages upon secondary contact within such a time scale Hewitt, 2003). However, our results support previous research suggesting hybridization among Arctic char glacial lineages. Atlantic lineage nuclear DNA has been found in Nunavut populations of Arctic lineage individuals (Moore et al., 2015). A similar lack of a relationship between mtDNA and nuclear DNA has also been observed in three-spine stickleback (Lescak et al., 2017). Many Salvelinus species are known to readily hybridize (Taylor, 2004), and there is evidence for Arctic char having hybridized with brook trout in Quebec (Bernatchez, Glémet, Wilson, & Danzmann, 1995;Glémet, Blier, & Bernatchez, 1998) and Labrador (Fraser River) (Hammar, Dempson, & Verspoor, 1991), and with lake trout in Nunavut (Wilson & Hebert, 1993) and Quebec . These hybridizations overcome a much older allopatric divergence than that among Arctic char glacial lineages. Hybridization among species does F I G U R E 3 Haplotype map of Arctic char (Salvelinus alpinus) haplotypes created with PopArt version 1.7 (Leigh & Bryant, 2015) using a Median-Joining network (Bandelt et al., 1999) and an Epsilon value of 0. New haplotypes identified in this study and Salisbury et al., 2018 are starred not necessitate the capacity for hybridization among intraspecific glacial lineages. However, given the relatively short duration of allopatric divergence, the lack of reproductive isolation among glacial lineages is unsurprising.
Some of the brook trout and lake trout mtDNA haplotypes detected in our samples may therefore reflect hybridization or backcrosses between these species and Arctic char. This would require further validation using nuclear markers but was beyond the scope of this study. An open area for future investigation is the degree to which genes from lake trout and brook trout have introgressed into Arctic char genomes within this region.

| Colonization history
We detected several rare haplotypes that were previously found in other populations within each lineage's respective range allowing for insight into the origins of the three glacial lineages in this region.
The ARC20 and ARC22 haplotypes we detected in Labrador were previously observed in geographically distinct locations across the high Canadian Arctic (Moore et al., 2015). The Arctic lineage may have therefore colonized Labrador multiple times from geographically distant populations. The ATL19 haplotype we observed in a single dark char morph in Gander Lake was previously observed in an unspecified morph in this lake as well as in a resident lacustrine population from Scotland (Moore et al., 2015). Lastly, the ATL31 haplotype we found in multiple anadromous populations was also found in a landlocked, Swedish population (Oleinik et al., 2017).
The appearance of these Atlantic haplotypes on opposite sides of the Atlantic Ocean suggests extensive colonization throughout the Atlantic from the Atlantic refugium. While our study area demonstrates a high diversity of Atlantic lineage haplotypes, this diversity is no doubt due to our intensive sampling. Whether the Atlantic refugium was located on the western or eastern side of the Atlantic therefore requires further investigation. Though they did not share a common lineage, most landlocked populations contained a single lineage and low haplotypic diversity.
This could be due to a founder-take-all scenario, where the lineage that first colonized a lake rapidly expanded to fill available habitat, preventing subsequent incursions from other lineages (Orsini et al., 2013;Waters, 2011;Waters et al., 2013). Also, landlocked populations are more isolated and tend to exhibit smaller effective sizes (Salisbury et al., 2018) and thus experience more drift than anadromous populations potentially leading to a greater loss of mtDNA haplotypes.
Several landlocked lakes countered this trend of reduced diversity. Landlocked lakes within the Kogaluk River system (i.e., V10, V11, V15, V16) had Arctic and Atlantic lineage char co-occurring.
Access to this watershed may have been enhanced by significant runoff from the paleolake Naskaupi, which drained through the Kogaluk between 7,500 and 6,000 years BP (Barnett & Peterson, 1964;Jansson & Kleman, 2004

| Glacial lineage and contemporary morph divergence in Gander Lake
Previous work has suggested that the high degree of neutral genetic differences observed between pale and dark morph char could be ascribed to differential glacial origins (Gomez-Uchida et al., 2008).
Our results, indicating that most char in Gander Lake were of the  Dempson, 2002;Power, O'Connell, & Dempson, 2005) may have arisen in sympatry in Gander Lake within the last ~10,000 years since its deglaciation (Bryson et al., 1969;Dyke, 2004;Shaw et al., 2006). This is consistent with the presumed sympatric divergence of other lacustrine Arctic char morphs (Gíslason, Ferguson, Skúlason, & Snorrason, 1999;Magnusson & Ferguson, 1987;Volpe & Ferguson, 1996). The large genetic divergence among pale and dark morph char in Gander suggests substantial genetic differences can accumulate between morphs within a short period of time, potentially fueled by divergent selection (Taylor, 2004)

| Management implications and the utility of intensive mtDNA sampling
The likely introgression among glacial lineages in Labrador has important implications for the char fishery in Labrador. There was evidence of Arctic, Atlantic, and even Acadian lineage fish in sea-accessible locations in the Notakwonan, Voisey, Anaktalik, Nain, and Okak drainages. These populations probably contribute to the commercial fishery stock complexes (DFO, 2001;Dempson et al., 2008).
The expected introgression between lineages suggests that there is likely no need to manage them separately; however, this should be further validated by investigating the relative lineage makeup of commercially caught char.
Our results verify the utility of intensive mtDNA sampling across many populations, particularly within a secondary contact zone. This approach facilitated the detection of a number of new haplotypes for the Arctic, Atlantic, and Acadian lineages (Figures 2 and 3)

DATA ACCE SS I B I LIT Y
Newly identified Arctic char mtDNA D-loop haplotypes were submitted for archival with GenBank (Accession Numbers: MK208868-