Demography and evolutionary history of grey wolf populations around the Bering Strait.

Glacial and interglacial periods throughout the Pleistocene have been substantial drivers of change in species distributions. Earlier analyses suggested that modern grey wolves (Canis lupus) trace their origin to a single Late Pleistocene Beringian population that expanded east and westwards, starting c. 25,000 years ago (ya). Here, we examined the demographic and phylogeographic histories of extant populations around the Bering Strait with wolves from two inland regions of the Russian Far East (RFE) and one coastal and two inland regions of North-western North America (NNA), genotyped for 91,327 single nucleotide polymorphisms. Our results indicated that RFE and NNA wolves had a common ancestry until c. 34,400 ya, suggesting that these populations started to diverge before the previously proposed expansion out of Beringia. Coastal and inland NNA populations diverged c. 16,000 ya, concordant with the minimum proposed date for the ecological viability of the migration route along the Pacific Northwest coast. Demographic reconstructions for inland RFE and NNA populations reveal spatial and temporal synchrony, with large historical effective population sizes that declined throughout the Pleistocene, possibly reflecting the influence of broadscale climatic changes across continents. In contrast, coastal NNA wolves displayed a consistently lower effective population size than the inland populations. Differences between the demographic history of inland and coastal wolves may have been driven by multiple ecological factors, including historical gene flow patterns, natural landscape fragmentation, and more recent anthropogenic disturbance.

INTRODUCTION

Climatic‐driven environmental changes and geological events profoundly influence species ranges and population demographic processes (Hewitt, 2000). During Pleistocene glacial periods, the expansion of ice sheets, coupled with a global drop in sea levels, led to a substantial increase in the exposure of the continental shelf, connecting previously isolated landmasses. One emblematic example with major consequences for Holarctic species was the exposure of the Bering Land Bridge (BLB), expanding the Beringia region from Eastern Siberia (Russia) to Western Yukon (Canada; Brigham‐Grette, 2001; Hopkins et al., 1982). The BLB remained predominantly ice‐free during glacial advances across its Eurasian and North American extensions, acting for many taxa as a glacial refugium and as an important source and route of post‐glacial colonization. Examples can be found among plants (Brubaker et al., 2005), insects (Elias & Crocker, 2008), fish (Campbell et al., 2015), and mammals (Anijalg et al., 2018; Heintzman et al., 2016), including humans (Watson, 2017). Yet, habitat and climate heterogeneity across the BLB shaped species ranges and promoted population structure (Dale Guthrie, 2001; Elias & Crocker, 2008; Kuzmina et al., 2011; Vershinina et al., 2021). The decline in sea level that exposed the BLB also reshaped the coastline of Pacific Northwest North America, connecting several islands to the mainland (Baichtal, 2010; Shugar et al., 2014). Studies combining phylogeography and reconstruction of Pleistocene habitats have suggested that South‐eastern Alaska (SEAk) may have supported a Pleistocene refugium, as synthesized in the Coastal Refugium Hypothesis (Carrara et al., 2007). This hypothesis highlights the diversity, divergence, and high levels of endemism of SEAk taxa, observed, for instance, in several species of mammals [e.g., Martes americana (Small et al., 2003), Mustela erminea (Colella et al., 2018) or Ursus americanus (Puckett et al., 2015)]. Following the Last Glacial Maximum (LGM), the Pacific Northwest coast may have provided an unglaciated and viable route connecting Beringia to central North America (Lesnek et al., 2018).

Whereas historical climate fluctuations have played a predominant role in shaping current phylogeographic patterns, the intensification of human‐related activities over the last centuries has increasingly influenced species distributions, abundance, and ecological interactions (Scheffers et al., 2016; Yackulic et al., 2011). Land use activities such as habitat conversion and harvesting are among the highest direct pressures on species (Newbold et al., 2015). The SEAk region is an example of how historical climatic fluctuations have produced major phylogenetic patterns (Sawyer et al., 2019; Stone et al., 2002), but recent anthropogenic changes to the landscape are promoting fine‐scale spatial shifts in species patterns. An example is the reported association between timber harvesting practices and local changes in the density of the Sitka black‐tailed deer (Odocoileus hemionus sitkensis) population (Albert & Schoen, 2013; Brinkman et al., 2011).

Grey wolves (Canis lupus) are among the few top predators that survived the Pleistocene's rapid environmental changes and the last centuries' intense anthropogenic pressures (Barnosky et al., 2004; Dufresnes et al., 2018). The species has a Holarctic distribution covering a wide range of environments, from tundra to deserts. Based on mitogenomes, Loog et al. (2020) suggested that extant wolf populations trace their origin to an expansion of the Beringian population starting c. 25 thousand years ago (ka) towards Eurasia and later, c. 15 ka, towards North America. This expansion led to the replacement of indigenous Pleistocene wolf populations (Loog et al., 2020). An alternative scenario for North America is that extant wolves on the continent descend from a pre‐LGM expansion out of Beringia followed by a post‐LGM northward expansion from a glacial refugium south of the Laurentide/Cordilleran ice sheet complex (Koblmüller et al., 2016; Meachen et al., 2020). Overall, mitochondrial signatures illustrate a complex phylogeographic history for extant North American and Eurasian populations and high genetic similarity between the wolves of Alaska and Eastern Russia (Ersmark et al., 2016; Koblmüller et al., 2016; Loog et al., 2020). Nuclear data, however, show two reciprocally monophyletic wolf clades (Fan et al., 2016; Pilot et al., 2019; Sinding et al., 2018), with signs of intercontinental admixture between East Siberian and Alaskan populations (Sinding et al., 2018), consistent with their higher geographic proximity to the area of common origin. Whereas nuclear data have been extensively explored for North American populations, few genetic studies have focused on wolves from the Russian Far East (hereafter RFE; but see Talala et al., 2020 and Vorobyevskaya & Baldina, 2011), limiting our understanding of demographic and phylogeographic patterns in this region. Similarly, inferences based on mitochondrial data by Loog et al. (2020) did not include modern RFE samples.

In contrast to the mitonuclear dichotomy between continents, a coherent pattern has emerged within North America. Both nuclear and mitochondrial data have indicated that wolves from the Pacific Northwest coast, including SEAk, represent a phylogeographic clade distinct from inland populations (Carmichael et al., 2001; Cronin et al., 2015; Schweizer et al., 2016; Weckworth et al., 2005, 2010, 2011). This differentiation was first hypothesized to represent the remnants of pre‐LGM colonization of North America that expanded into SEAk by the early Holocene, whereas wolves expanding from the Beringian refugium after the LGM would have contributed to the extant North‐western North American (hereafter NNA) populations (Weckworth et al., 2010). However, the recent study by Loog et al. (2020) comparing ancient and contemporary mitogenomes found that extant SEAk wolves carry the same mtDNA lineage as the wolves that expanded out of Beringia after the LGM. Additionally SEAk wolves have very low levels of genetic diversity compared with their inland counterparts, suggesting either reduced gene flow with the nearest inland populations or extensive genetic drift due to the small founding population size (Weckworth et al., 2011). However, previous demographic inferences based on microsatellite data showed no signs of recent or historic bottlenecks in the SEAk population (Weckworth et al., 2005).

In this study, we present a demographic and population genetic analysis of extant wolf populations from both sides of the Bering Strait. Our goal was to assess genetic diversity and structure patterns for these populations and understand the contribution of historical and contemporary processes shaping their demographic history. We hypothesized that (1) RFE and NNA populations diverged before the expansion out of Beringia (Loog et al., 2020), in accordance with the extensive climatic and geomorphological heterogeneity across Beringia during this period. Furthermore, we hypothesized that (2) the known differentiation between inland and coastal NNA wolves is more ancient than would be anticipated by their shared mitochondrial lineage (Loog et al., 2020), given the high molecular and ecological support for divergence (e.g., Weckworth et al., 2011). Finally, we hypothesized that (3) wolves from RFE, and inland NNA have higher genetic diversity than those from SEAk, which are geographically limited by and have probably suffered negative impacts from recent anthropogenic changes to landscapes. To test our hypotheses, we used a comprehensive sample of wolves from populations around the Bering Strait, and genome‐wide single nucleotide polymorphism (SNP) profiles, to (a) assess the structure and genetic diversity of these populations, (b) estimate divergence times among populations, (c) infer past and recent fluctuations in effective population size (N e), and (d) identify genetic signatures of recent demographic changes.

MATERIALS AND METHODS

Sampling and genotyping

We generated original data for 65 grey wolves from two regions of the RFE, 42 individuals from Yakutia and 23 from Chukotka (Figure 1; Table S1). Muscle tissue samples were collected from wolves found dead or legally harvested for purposes other than research in the study area, where the species is classified as a game animal. CIBIO/InBIO is a CITES certified institution (PT008). DNA extracts from these samples were obtained using the DNeasy Blood and Tissue Kit (Qiagen). Subsequently, samples were genotyped with the CanineHD BeadChip microarray (Illumina, Inc.) comprising more than 170,000 genome‐wide SNPs. Genotype calling was performed using GenomeStudio software (Illumina) following Illumina's guidelines (GenomeStudio, Genotyping Module version 2.0 Software Guide, 2016), and exported to plink format (Purcell et al., 2007) with a final total genotyping rate of 0.961. Additionally, genome‐wide profiles from 105 NNA wolves genotyped with the same SNP panel by Cronin et al. (2015) and Medrano et al. (2014) were included in our data set. Samples from NNA comprised two inland regions, Alaska (n = 37) and British Columbia (BC, n = 31), and the coastal region of SEAk (n = 37, coastal region plus Kopreanof Island, Figure 1; Table S1). The data sets were homogenized based on the “TOP/BOT” SNP calling method (Illumina, Inc., 2006), and merged using plink 1.9 (Purcell et al., 2007). Monomorphic loci and loci with greater than 5% missing data were excluded from the final merged data set. Relatedness among individuals within each of the five sampling regions was assessed using identity‐by‐descent distances calculated in plink 1.9. Individuals with pairwise relatedness higher than r = 0.30 and more than 10% missing genotypes were excluded. To ensure that wolves included in this study were genetically distinct from dogs, we also performed a preliminary principal component analyses (PCA) using the adegenet package (Jombart, 2008) in R 4.0.3 (R Core Team, 2020), including all wolves used in this study plus 36 mixed breed dog genotypes for the same SNP panel from Cronin et al. (2015) and Medrano et al. (2014) (Figure S1). The final data set included 147 wolves genotyped for 91,327 autosomal SNPs (henceforth 91K).

FIGURE 1

Map of the Bering Strait and contiguous areas in Eurasia and North America showing the approximate geographical locations of the samples and sample sizes of grey wolf populations used in this study: Yakutia (n = 42) and Chukotka (n = 23; Russian Far East; RFE); Alaska (n = 37) and British Columbia (n = 31; North‐western North America; NNA). The inset represents an enlargement of the study area in South‐eastern Alaska (n = 37; NNA), where the two stars represent the coastal region and Kopreanof Island. The pale blue indicates regions with sea levels shallow enough to have been exposed land during the Last Glacial Maximum (sea depth <120 m).

Population structure and differentiation

To determine the number and composition of the discrete populations present across our study area, we used a combination of likelihood clustering analyses and multivariate methods. For both approaches, we pruned the 91K data set for loci in high linkage disequilibrium (LD, r 2 ≥ .5) using plink with a 50 SNP sliding window and a step size of 10 SNPs, resulting in a final data set of 55,994 SNPs (henceforth 56K). We used the maximum likelihood approach implemented in admixture (Alexander et al., 2009) to evaluate distinct models for up to 10 populations (K) with a bootstrap of 5000 and 10 independent replicates per K‐value. To determine the best value of K, we performed 10 cross‐validations and selected the K‐value with the lowest validation error. The fit of the inferred model was assessed based on the estimated correlation of the residual difference between individuals using the software evaladmix (Garcia‐Erill & Albrechtsen, 2020). Additionally, we performed a PCA using the adegenet package (Jombart, 2008) in R 4.0.3 (R Core Team, 2020). Finally, the eigensoft software (Patterson et al., 2006) was used to estimate F ST between pairs of populations, and the average divergence between and within populations, estimated based on the normalized distance between each pair of individual eigenvectors.

Population genetic diversity

We used the 91K data set to calculate the observed (H o) and expected heterozygosity (H e) per population in plink and performed Wilcoxon–Mann–Whitney (WMW) tests comparing population heterozygosity between all population pairs, with statistical significance set at p < .05. The same software was used to estimate the LD decay for each population by calculating the genome‐wide pairwise genotype association coefficient (r 2) between all SNP pairs. Since this statistic is affected by the sample size, we randomly selected 10 individuals from each population. Additionally, we characterized the individual and population autozygosity (F HBD) using the R package RZooRoH (Bertrand et al., 2019). This package implements a model‐based approach to identify homozygosity‐by‐descent (HBD, obtained from runs of homozygosity; RoH) and non‐HBD segments across the genome. RZooRoH categorizes the length of HBD segments into generation classes, providing information on the timing of past inbreeding events. Longer segments correspond to recent inbreeding events (few generations ago) whereas shorter segments correspond to ancestral inbreeding events (many generations ago). Each class is associated with a rate (R k) equal to the size of the inbreeding loop in generations (Druet & Gautier, 2017). The rate of the class is approximately equal to twice the number of generations to the ancestor or inbreeding event. Based on the density of our SNP panel, we applied a MixKR model with 11 HBD classes (with R k equal to 2, 4, 8, 16, 32, 64, 128, 256, 512 and 1024) and 1 non‐HBD class. To convert from recombination distance to physical distance, we used the average recombination rate of 0.97 cM/Mb for the dog genome (Campbell et al., 2016). We performed WMW tests comparing population inbreeding coefficients between all population pairs, with statistical significance set at p < .05.

Phylogenetic relationships and divergence dating

We inferred the phylogenetic relationship among populations and associated divergence times in a coalescent framework using the snapp 1.3 package (Bryant et al., 2012) implemented in beast2 2.4 (Bryant et al., 2012), following Stange et al. (2018). To ensure independence among loci, only loci distanced by 100 kb (establish based on the LD decay pattern) were used, resulting in a data set of 15,340 SNPs. Due to computational constraints, each population was represented in this analysis by two individuals randomly selected among those satisfying the following three criteria: (1) highest assignment to the respective populations in the admixture analysis, to avoid violation in the species tree assumptions, namely the effect of gene flow; (2) lowest pairwise‐relatedness; and (3) lowest inbreeding coefficient, using the population average as a reference value. As a calibration point, we used the divergence between wolves from Europe and Asia previously estimated to c. 29 ka (Fan et al., 2016). To this end, we added two genotypes of Croatian wolves genotyped by Stronen et al. (2013) for the same SNP panel. We placed an age constraint on the node that separates the RFE and Croatian wolves using a lognormal distribution with an offset of 0, a mean of 29 ka, and a standard deviation (in real space) of 0.04 (Fan et al., 2016). We performed three independent runs in snapp each with a length of 10 million MCMC generations. Runs were assessed using tracer 1.6 (Rambaut et al., 2018) to examine convergence, and the tree topology and node heights were visualized using densitree (Bouckaert, 2010). We examined the effect of sample selection on the estimated phylogeny by repeating the analysis with a distinct set of individuals selected based on the same criteria.

Changes in effective population size over time

To infer the demographic history of the populations around the Bering strait, we used the model‐free approach implemented in stairway plot 2.1 (Liu & Fu, 2015; Liu & Fu, 2020) to reconstruct historical changes in the N e of each population. According to the software recommendations, potential confounding effects from linkage disequilibrium were minimized by performing an LD‐based variant pruning of the 91K data set using plink with 50 SNP sliding windows, a step size of 10 SNPs, and a threshold of r 2 = .8. For each population, the input folded site frequency spectrum (SFS) was estimated using the script easySFS.py (https://github.com/isaacovercast/easySFS; Figure S2). The output was scaled to time and population sizes assuming a mean generation time of 4.5 years (Mech et al., 2016) and a mutation rate of 4 × 10−9 per site per generation (Skoglund et al., 2015). Demographic changes were estimated using the median of simulations based on 200 bootstrap replicate analyses, and the precision of estimates was evaluated using 95% confidence intervals (CI).

To reconstruct the more recent demographic history (<1 ka) of these wolf populations, we used the software gone (Santiago et al., 2020). gone implements a genetic algorithm (Mitchell, 1998) to optimize the inference of the historical effective population size series that best fits the observed LD between pairs of SNPs with different recombination rates. Forty independent replicates of the gone analysis were run using standard input parameters, with the average recombination rate set to 0.90 cM/Mb (Campbell et al., 2016). N e estimates were obtained for each population individually using a data set with no missing data and no minor allele pruning. Notwithstanding, we note that the ascertainment scheme of the SNP chip used may affect estimates of the absolute N e, but not its trend (Lachance & Tishkoff, 2013). As the method does not generate parametric CIs of the estimated N e, and our sample size per population was not large enough to subsample individuals, we implemented a chromosome subsampling approach to calculate empirical CIs. For each individual, we took 50 random subsamples of half the number of chromosomes (only considering autosomes; n/2 = 19), to create genomic subsampled data sets of each population, containing the original number of individuals but half the number of chromosomes. All resulting subsampled data sets were analysed independently as described above.

RESULTS

Population inferences

The maximum‐likelihood clustering analysis based on 56K SNPs exhibited very similar low cross‐validation errors for the models with three, four and five populations. The model for K = 3 separated the RFE, inland NNA and SEAk populations. Models for K = 5 further subdivided the RFE and inland NNA into four distinct clusters, in accordance with the geographical distribution of the samples (Figures 1 and 2a). Results from evaladmix showed that the five‐population model had the lowest level of residuals (Figure S3). The PCA based on the same data set indicated a clear distinction between continents along the first axis while the second axis further differentiated populations within NNA (Figure 2b). Independent PCA for each continent corroborated the differentiation among the five populations, consistent with the admixture analysis (Figure 2c,d). Despite the partition of samples in different populations according to geographic origin, the results show a clear signature of gene flow between neighbouring populations within each continent, and residual signs of intercontinental shared ancestry. However, as downstream genetic inferences require the assumption of discrete populations, we hereafter considered the five most likely populations: Yakutia, Chukotka, Alaska, BC, and SEAk.

FIGURE 2

Analyses of population structure of grey wolves around the Bering Strait based on 56K SNPs. (a) admixture results for K = 5. (b) Principal component analysis (PCA) using all sampled individuals. (c) PCA using only North American individuals; (d) PCA using only Russian Far East individuals. SNP, single nucleotide polymorphism.

The highest pairwise F ST values were observed between RFE populations and SEAk, whereas comparisons within RFE and inland NNA populations showed the lowest values (Table 1). Pairwise estimates between RFE and inland NNA populations are very similar to the values observed between inland and coastal NNA populations. Average genetic divergence between populations showed the highest values between RFE and NNA populations and, unlike F ST estimates, similar pairwise values between the three NNA populations (Table 1). Average divergence statistics as estimated by eigensoft do not correct for intrapopulation divergence, which is probably the reason for the observed differences between these estimates and F ST. The highest average divergence within a population was observed in Yakutia, whereas the lowest was observed for the SEAk population (Table 1). The remaining populations showed intermediate values, closer to 1.

TABLE 1

Pairwise genetic differentiation (F ST) and divergence between wolf populations calculated in EIGENSOFT for the five sampling regions

	Russian Far East		North‐western North America
	Yakutia	Chukotka	Alaska	British Columbia	Southeast Alaska
Yakutia	1.28	0.04	0.15	0.14	0.24
Chukotka	1.27	1.09	0.17	0.16	0.26
Alaska	1.42	1.36	0.95	0.05	0.17
British Columbia	1.43	1.37	1.06	0.99	0.14
Southeast Alaska	1.46	1.40	1.11	1.06	0.69

Genetic diversity

Individual observed heterozygosity ranged from 0.12 to 0.31 [SD]; Figure 3a), with the Yakutia population showing the highest values and the SEAk population the lowest. Wilcoxon‐Mann‐Whitney tests between population pairs were significant (p < .05) for all except for Alaska and BC. The four inland populations (Yakutia, Chukotka, Alaska, and BC) showed LD decay curves reaching r 2 = .3 after 25 kb, whereas the LD decay of SEAk was much slower, only reaching values of r 2 = .3 after 200 kb (Figure 3b).

FIGURE 3

Genetic diversity patterns for the two Russian Far East (RFE): Yakutia (YAK) and Chukotka (CHU); and the three North‐western North American (NNA) wolf populations: Alaska (Ak), British Columbia (BC) and Southeast Alaska (SEAk). (a) Average and distribution of observed heterozygosity; (b) Extent of linkage disequilibrium, represented by the average genotypic association coefficient r 2 as a function of inter‐SNP distance; (c) Average and distribution of individual F HBD. Significance was calculated performing a Wilcoxon–Mann–Whitney test: ns, not significant, *significance at p < .05. HBD, homozygosity‐by‐descent; SNP, single nucleotide polymorphism.

SEAk wolves showed the highest average inbreeding coefficients (F HBD), significantly different from all other populations (Wilcoxon‐Mann‐Whitney test, p < .05), as well as the highest variance for this parameter, ranging from 0.15 to 0.57 (; Figure 3c). In contrast, the four inland populations showed very reduced levels of F HBD. To explore the influence of inbreeding (measured by F HBD values) on population differentiation in SEAk wolves, we compared the PCA coordinates of each individual with its respective F HBD value following Smeds et al. (2020). Some of the most divergent SEAk individuals exhibited high F HBD values, but there was no clear association between an individual's F HBD and its position on the PCA (Figure S4). This suggests that the observed population differentiation is not a direct result of genetic drift associated with high levels of inbreeding.

Phylogenetic relationships and divergence times

The phylogenetic relationship among populations inferred through Bayesian coalescent analysis showed a primary division between NNA and RFE populations (Figure 4; Figures S5 and S6). This divergence was dated to c. 34.4 ka (95% highest posterior densities [HPD] 31.4–37.7), around 10 ka prior to the LGM. The relationship among NNA populations showed a first divergence between coastal and inland populations, followed by the divergence between inland Alaska and BC. The split between SEAk and inland populations occurred around 16 ka, not long after the LGM (Clark et al., 2009). The time of divergence for mainland populations within NNA was very similar to that observed for populations in RFE, happening around 10 ka, and time estimates for the two events have overlapping CIs (Figure 4). For both NNA and RFE, this population divergence was also broadly consistent with the Pleistocene–Holocene transition (Figure 4).

FIGURE 4

Chronogram of wolf populations from Bayesian coalescent analysis of SNP data using snapp. Median ages in thousand years ago are provided above nodes, with the 95% highest posterior densities (HPD) represented by grey bars in the nodes. The x‐axis corresponds to time before present in thousand years. The vertical dashed lines indicate the intervals of the last glacial maximum (LGM) and the estimated time of flooding of the Bering Land Bridge (BLB). SNP, single nucleotide polymorphism.

Demographic history

For RFE and inland NNA populations, stairway plot reconstructed continuous processes of demographic decline until the mid‐Holocene, shifting to a stable demographic trend with constant N e between 5 and 4 ka (Figure 5a; Figure S7). During the last thousand years, the four inland populations showed a signature of low N e (Figure 5a; Figure S7). For the SEAk population, stairway plot inferred consistently smaller N e values compared with the other four inland populations, and a stable demographic trend with constant N e until c. 15 ka when the population started to decline (Figure 5a; Figure S7).

FIGURE 5

Reconstruction of demographic history for each wolf population with historic and recent estimates of effective population size (N e). (a) Historic estimates of temporal N e inferred using stairway plot and (b) recent (<1 ka) estimates of temporal N e inferred using gone. Lines indicate the median estimated N e, plotted from the past to the present. Full lines depict Russian Far East (RFE) populations and dashed lines North‐western North American (NNA) populations. The 95% confidence intervals associated with each temporal estimate are shown in Figures S7 and S8.

Results from the GONE analysis increased the resolution of the most recent (<1 ka) N e fluctuations, revealing an overall declining pattern across populations, particularly over the past 300 years (Figure 5b; Figure S8). Inferences for the two RFE populations showed that both suffered a sharp decline c. 150 years ago (ya). For Yakutia this decline was more pronounced, as this population had reached a considerably larger size than Chukotka in the preceding centuries. Among inland NNA populations, BC exhibited a more recent population decline, between 350 and 150 ya, than Alaska, which declined between 600 and 400 ya (Figure 5b; Figure S8). However, the latter exhibited a subsequent decline c. 100 ya. SEAk wolves showed a less sharp and more continuous population decline beginning c. 600 ya and lasting until the present (Figure 5b; Figure S8).

Characterization of individual levels of inbreeding

To further explore possible genetic consequences of the demographic history in these populations, we characterized the individual observed HBD segments based on segment lengths (Figure 6). This analysis revealed similar patterns between neighbouring populations (Yakutia/Chukotka and Alaska/BC), and a very different signature for SEAk. In RFE populations, average HBD segment length was 1.6 ± 3.3 Mb for Yakutia and 1.9 ± 3.1 Mb for Chukotka, with most fragments being smaller than 1 Mb (median of 7 and 9 kb respectively, Figure 6a). Accordingly, individual proportions of HBD segments were mostly associated with HBD classes with high rates, R k = 256 for Chukotka and R k = 512 for Yakutia (Figure 6b), meaning that most of the HBD segments observed in these populations trace to ancestors c. 128 and 256 generations ago for Chukotka and Yakutia, respectively.

FIGURE 6

Characterization of grey wolf individual and population inbreeding levels. (a) Distribution of the length in Mb of the homozygosity‐by‐descent (HBD) segments in each population. The inner part of each violin plot represents the mean and the standard deviation of the distribution. (b) Genomic proportion of each population associated with different HBD classes. Each bar represents the average individual genomic proportion per population associated with a specific HBD class. The rate (R k) of a HBD class is approximately equal to twice the number of generations since the common ancestor. For instance, a class with a R k equal to eight correspond to ancestors living approximately four generations ago. (c) Partitioning of individual genomes in different HBD classes. Each bar represents an individual, and the height of the different stacks represents the proportion of the genome associated with each HBD class. The colours of the bars represent the rates (R k) associated with each HBD class.

In NNA, both Alaska and BC exhibited an average HBD length of 2.8 ± 5.4 Mb, but as for RFE populations, most of the segments were smaller than 1 Mb (median = 1 Mb, Figure 6a). However, unlike the RFE populations, both Alaska and BC showed a wider range of ancestral inbreeding events (Figure 6b,c). Most individuals showed evidence of inbreeding events occurring 128 or more generations ago (R k ≥ 256, Figure 6b), but a small number of individuals presented longer HBD fragments tracing their origin to ancestors 16 or fewer generations ago (R k ≤ 32). SEAk wolves showed the highest genome‐wide proportion of HBD segments, which reflects an average segment length of 4.5 Mb ± 2.1 and a median of 2.1 Mb (Figure 6a). Unlike the other populations, there were few signs of inbreeding events occurring more than 64 generations ago (R k = 128). Instead, all SEAk individuals presented a high proportion of longer HBD segments throughout the genome, originating from ancestors 32 or fewer generations ago (R k ≤ 64, Figure 6b,c).

DISCUSSION

The main goal of this study was to provide a comprehensive overview of the diversity, differentiation, and demographic history of the extant wolf populations around the Bering Strait and expand our knowledge about the RFE populations. By describing the past fluctuations in N e for each population, we found a temporal synchronism in the demography of inland populations from both continents. In contrast, the coastal population of SEAk showed a distinct demographic pattern with a consistently lower population size compared to its continental counterparts. Divergence time estimates showed that NNA and RFE populations split prior to the LGM (c. 34 ka). The oldest divergence event within the continents corresponded to the split between SEAk and inland NNA populations dating to c. 16 ka. This estimate indicates that an ancient divergence associated with multiple expansions out of Beringia may have contributed to the deep differentiation between inland and coastal NNA populations. We also found a reduction in genetic diversity across a geographical gradient from RFE towards NNA, with SEAk wolves showing the lowest genetic diversity and highest genome‐wide homozygosity, associated with a recent decline in population size.

The history of population divergence and demographic trends

Our results show that RFE and NNA populations shared a common ancestry up to 34.4 (95% HPD 31.4–37.7) ka, and that their divergence followed a demographic decline that started c. 60 ka and lasted until the mid‐Holocene. This is concordant with our first hypothesis stating that divergence between wolf populations in Asia and North America predates the suggested 20 to 15 ka expansion out of Beringia (Loog et al., 2020). The difference in divergence time estimates may suggest the existence of genetic structure in the Beringian population prior to expansion, a pattern that has also been suggested for the Beringian moose population (Alces alces, Dussex et al., 2020). Such a pattern could be driven by habitat and climate heterogeneity between Western and Eastern Beringia during the late Pleistocene, which over time shaped species ranges and densities across Beringia (Dale Guthrie, 2001; Elias & Crocker, 2008; Kuzmina et al., 2011). This landscape heterogeneity has also been associated with a decline in the migration rate between Western and Eastern Beringian populations of caballine horses (Equus spp.) during the Middle and Late Pleistocene (Vershinina et al., 2021). Our results seem to agree, as RFE and inland NNA populations exhibit continuous trends of population decline in the period before and after the LGM. Although wolves are highly mobile, genetic differentiation of subpopulations living in close proximity can occur owing to habitat features, prey specialization or other ecological and geographic interactions (Schweizer et al., 2016; Silva et al., 2018; Stronen et al., 2014). Notably, inferences based on stable isotope analysis of ancient Beringian wolves have suggested some level of individual specialization associated with local prey availability (Fox‐Dobbs et al., 2008; Leonard et al., 2007; Meachen et al., 2020). Therefore, fluctuations in prey availability and habitat heterogeneity could have triggered the observed decrease in N e through local extinctions and subsequent divergence between ancestral populations of RFE and NNA wolves within Beringia. Nevertheless, it is also possible that the observed divergence between RFE and NNA populations resulted from a pre‐LGM expansion out of Beringia through the still open Yukon corridor towards central North America (Batchelor et al., 2019). This scenario is in line with Koblmüller et al. (2016) findings based on mitochondrial data, which suggested that extant North American wolves descend from wolves that persisted south of the ice sheets during the LGM. Additionally, fossil records from Wyoming (USA) suggest that Beringian wolves were present in this area prior to the LGM (Meachen et al., 2016). We note, however, that this alternative hypothesis is difficult to reconcile with recent findings of mitochondrial divergence observed for extant wolves (Loog et al., 2020).

Following the divergence between continents, NNA wolves differentiated further into distinct coastal and continental populations, with their common ancestor dating to c. 16 (95% HPD 14.4–18) ka. In contrast to our second hypothesis, this estimate places the arrival of modern wolves into the SEAk coastal region after the LGM. Assuming a within‐Beringia divergence scenario, the estimated arrival is in synchrony with the previously suggested eastward expansion out of Beringia estimated to c. 16 ka based on mitogenomes (Loog et al., 2020). After the LGM, the continental route between the Laurentide and Cordilleran ice sheets connecting Beringia to continental unglaciated North America became ecologically viable c. 13 ka (Heintzman et al., 2016; Pedersen et al., 2016). However, the retraction of the eastern front of the Cordilleran ice sheet was much faster, creating an ecologically viable coastal corridor along the Pacific Northwest coast as early as 17 ka (Darvill et al., 2018; Lesnek et al., 2018). Archaeological and genetic evidence support the hypothesis that this coastal route was used by the first humans migrating into the Americas around 17 ka (Becerra‐Valdivia & Higham, 2020; Moreno‐Mayar et al., 2018; Pedersen et al., 2016). The same pattern is also being suggested for the early lineages of domesticated dogs, which diverged from the Eastern Siberian lineage c. 16 ka (da Silva Coelho et al., 2021). Other large carnivores, such as brown bears (U. arctos), could also have used this coastal route to reach the SEAk region c. 17 ka (Salis et al., 2021). Based on the similarity between our estimated divergence time for inland and coastal NNA populations and those previously described for other species, and the concordance with geomorphological evidence of a viable coastal route, wolves may also have taken such a route when expanding out of Beringia towards continental North America. Notwithstanding, assuming a northward expansion from a southern refugium, the same unglaciated coastal corridor could have been used prior to the opening of the continental route between the Laurentide and Cordilleran ice sheets (Lesnek et al., 2018; Weckworth et al., 2010).

Demographic inferences over the last millennium revealed a common negative trend across populations. For both RFE populations, we observed a synchronized decline c. 150 ya, which coincides with an increase in historical hunting records and the implementation of a country‐wide bounty system in the beginning of the 19th century aiming for a considerable reduction in wolf numbers (Baskin, 2016; Graves, 2007). In contrast, the earlier decline observed for NNA populations could have been associated with the climatic oscillations during the Little Ice Age, an unusual cold period that lasted from the early 14th century through the mid‐19th century (Forbes et al., 2020). During this period, expansion of the ice sheets and glaciers along with a decrease in temperatures most probably influenced species ranges, movements, and resource availability (e.g., caribou; Gigleux et al., 2019). Additionally, anthropogenic land use, such as fire practices by Indigenous communities, particularly along the coast (Hoffman et al., 2017), and extensive logging, agriculture, and livestock grazing associated with the arrival of Europeans into the region towards the end of the 18th and 19th centuries (Turner et al., 2013), could have contributed directly or indirectly to the observed population declines. However, the extent to which environmental and/or anthropogenic effects may have influenced wolf numbers in NNA populations is unclear and would require further investigation.

Genetic diversity and the consequences of population decline

RFE and inland NNA wolves exhibited similar heterozygosity and inbreeding coefficients, comparable with previously described values for other extant populations in Eastern Europe and North America (Pilot et al., 2014; vonHoldt et al., 2011). In accordance with our third hypothesis, SEAk wolves showed the lowest levels of genetic diversity and higher inbreeding coefficients, revealing a genetic signature of persistent decline in effective population size, particularly over the past 64 generations (c. 280 years). The same dichotomy between inland and coastal populations was observed for the LD decay patterns. Inland populations showed a rapid decay with physical distance, suggestive of a historically large N e comparable with the patterns described for wolf populations from Eastern Eurasia and the north‐western United States (Pilot et al., 2014, 2019). In contrast, the coastal population showed a much stronger LD and considerably slower LD decay, suggestive of a small and inbred population (Hedrick et al., 2014).

Observed differences in diversity and LD for inland and coastal populations may, at least partially, result from the landscape. Inland populations are not geographically isolated, and gene flow between neighbouring populations within each continent is not limited, as suggested by the low pairwise F ST values. In these populations, the presence of gene flow, along with a reduction in persecution over the last century (Bragina et al., 2015; Hayes & Harestad, 2000; Mech, 1970), may have buffered the effects of a past demographic decline (Hayes & Harestad, 2000) and promoted stabilization at a lower population size. In contrast, coastal SEAk wolves inhabit a naturally fragmented landscape, bordered to the east by the Coast Mountains, which limits gene flow with neighbouring populations (Person et al., 1996; Schweizer et al., 2016). Over the last decades, this region has also experienced severe changes due to human activities such as logging and road construction, reducing habitat and prey availability and limiting connectivity in an already naturally fragmented habitat (Albert & Schoen, 2013; Roffler et al., 2018). Unlike other populations, SEAk revealed a continuous population decline, starting c. 600 ya and lasting until the present. Together with landscape and ecological factors, this decline could have triggered the observed increase in inbreeding events over the past 280 years, resulting in increased homozygosity. Such trends of genetic erosion and reduced gene flow are associated with an increased risk of inbreeding depression (Newman & Pilson, 1997; Roelke et al., 1993), as earlier described for Scandinavian and Isle Royale wolves, where inbreeding depression threatens population persistence (Kardos et al., 2018; Robinson et al., 2019).

Conclusions and future perspectives

Large predators play a key role in structuring terrestrial ecosystems through an interplay of top‐down and bottom‐up forces along trophic associations. Such a role has been described across modern and historical ecosystems, from the regulation of Pleistocene megaherbivore population sizes preventing habitat destruction associated with extensive grazing (e.g., Van Valkenburgh et al., 2016) to the trophic cascade effect of wolf reintroduction in Yellowstone National Park at the end of the 20th century after a 70‐year absence (Ripple & Beschta, 2012). Illuminating the demographic history of large carnivore populations at different temporal and spatial scales can help infer past evolutionary changes in species, communities, and ecosystems and inform conservation and management actions.

Evaluations of demographic fluctuations of grey wolf populations have mainly been focused on global or broad‐scale regions and historical periods, making regional evaluations and more recent periods difficult to extrapolate. Our results offer a population‐level view of the evolutionary history of the extant wolves living around the Bering Strait. We observed similar demographic trends among the populations of NNA and RFE wolves throughout the Pleistocene and Holocene, probably reflecting ancient large‐scale climatic fluctuations and more recent anthropogenic factors. Furthermore, we ascertained the differentiation of SEAk wolves, which is concordant with the colonization of this region through the Pacifc North‐western coastal route post‐LGM. Finally, we demonstrated that SEAk wolves present signs of genetic erosion associated with high inbreeding and a continuous population decline in recent times.

Nevertheless, future studies, including more comprehensive geographical, temporal, and genomic sampling, would be required to further clarify the observed population trends. For instance, a broader temporal sampling of the Alaskan wolf population would be required to tease apart the relative role of different factors that may have contributed to the distinct demographic history that this population has exhibited over the last centuries. Additionally, as SNP chip data could be affected by the ascertainment scheme used in variant discovery, we acknowledge that expanding the genetic information through whole‐genome data would allow further exploration of the patterns described here by implementing other methodological approaches such as model‐driven demographic inferences. Lastly, a detailed study of the SEAk population, including genetic and ecological data, would be required to test the hypothesis of inbreeding depression and its impact on the persistence of this ecologically and genetically distinct evolutionary lineage adapted to the unique coastal habitat of the Pacific Northwest.

AUTHOR CONTRIBUTIONS

Raquel Godinho and Astrid Vik Stronen conceived the idea; Carolina Pacheco carried out the data analysis; Bogumiła Jędrzejewska, Kamila Plis, Innokentiy M. Okhlopkov, Nikolay V. Mamaev and Sergei V. Drovetski contributed samples or data; Carolina Pacheco wrote the manuscript with input from Raquel Godinho and Astrid Vik Stronen and other coauthors reviewed the text.

CONFLICT OF INTEREST

The authors declare no competing interests.

OPEN RESEARCH BADGES

This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.17605/OSF.IO/3GTNM.

DATA AVAILABILITY AND BENEFIT‐SHARING STATEMENT

The wolf genotypes produced in this study can be accessed through The Open Science 515 Framework using the following link: https://doi.org/10.17605/OSF.IO/3GTNM.

Benefits Generated: A research collaboration was developed with the scientists who provided genetic samples, with all collaborators included as coauthors of this manuscript. Additionally, the results of this research have been shared with the provider and broader scientific community through public databases, as described above.

Data S1

Click here for additional data file.