Severe population bottleneck and cranial morphology change in the Mednyi Island subspecies of Arctic fox Vulpes lagopus (Carnivora: Canidae)

Introduction

Arctic foxes, Vulpes lagopus L., 1758, living on Mednyi Island suffered a drastic decline in population size in the late 1970s due to an outbreak of mange epizootic (Goltsman et al. 1996, 2005, Ploshnitsa et al. 2013). For several years, cub mortality was at 90% (Goltsman et al. 1996, Goltsman and Kruchenkova 2001). This constituted a population bottleneck, where variability was lost randomly in the population (Nei et al. 1975, Tajima 1989, 1996, Bouzat 2010). After the bottleneck the population stabilized (Goltsman et al. 2011), although numbers remained much lower than before, with only about 90 individuals in 2005 (Goltsman et al. 2005), sharply contrasting with up to 1000 individuals before the bottleneck (Geptner and Naumov 1967), rendering this subspecies endangered (Goltsman et al. 1996). Demographic bottlenecks result in a loss of genetic variation and increased inbreeding (Wright 1969, Theodorou and Couvet 2006). Genetic studies have shown that the current Mednyi population displays low variability (Dzhykiya et al. 2007, Geffen et al. 2007, Ploshnitsa et al. 2012), a probable consequence of the bottleneck in the 1970s (Ploshnitsa et al. 2012). The disease also occurred on the nearby Bering Island, but it had no appreciable effect on the Arctic fox population size, which remained stable at about 600 adult animals (Ryazanov 2002), compared to 2000–4000 Arctic foxes in the mid-twentieth century (Geptner and Naumov 1967). Bering Island, at an area of 1667 km², is the bigger of the two Commander Islands supporting Arctic foxes, and can presumably accommodate a larger Arctic fox population than Mednyi Island, at only 186 km² (Ploshnitsa et al. 2013).

Both Mednyi and Bering Islands Arctic fox populations have been isolated from the mainland population of Arctic foxes (and most likely from each other) since the ice cover retreated after the Last Glacial Period, 10,000 years ago (Goltsman et al. 1996, 2005, Geffen et al. 2007, Dzhykiya 2008). It has been shown (Nanova 2009) that the cranial morphology of Mednyi Arctic foxes differs not only from the mainland population, but also from the neighboring population on Bering Island. Although overall cranial length is similar for both Bering Island and Mednyi Island Arctic foxes, the rostrum of Mednyi Arctic foxes is relatively shorter in comparison to Bering Island (Nanova 2009, Nanova and Prôa 2017). This relative shortening of the rostrum in Mednyi foxes has been previously interpreted as directly related to the biomechanical requirements of hunting prey (Nanova and Prôa 2017). Yet, after passing through the population bottleneck the diet of Mednyi Arctic foxes has changed considerably: Northern fulmars (Fulmarus glacialis L., 1761) and storm petrels (Oceanodroma furcata Gmelin, 1789 and Oceanodroma leucorhoa Vieillot, 1818) became the main food source, and the consumption of other bird species, marine invertebrates and otarid products, foodstuffs equally available as before, substantially decreased (Goltsman et al. 2011, Bocharova et al. 2013). These changes are expected only to occur over longer periods of time, and it is possible that passing through a bottleneck has randomly affected the genetic structure of the population and consequently altered the expressed phenotypic traits. Indeed, genetic drift may play a role in peak shifts on phenotypic adaptive landscapes, but the extent to which it contributes to phenotypic differentiation in adaptive radiation is not known (Schluter 2000, p. 119). Arctic foxes in the Scandinavian Peninsula have also suffered a recent bottleneck (Dalén et al. 2006, Nyström et al. 2006) which created population differentiation within Scandinavia that did not exist prior to the bottleneck.

Bottleneck detection is, therefore, critical for the interpretation of the demographic history of endangered populations, with consequences for conservation management (Ploshnitsa et al. 2013). Methods for detecting bottlenecks were devised for use with molecular genetic markers (Harpending and Jenkins 1973, Cornuet and Luikart 1996), but it is possible to study the pattern of phenotypic divergence in the same way as neutral genetic data (Spurgin et al. 2014). The possibility of using morphometric data means that material from past populations (where no DNA can be retrieved) can be studied and compared with material from current populations. Studies using morphometric traits as proxies for genetic data (Relethford and Lees 1982, Williams-Blangero and Blangero 1989, Williams-Blangero et al. 1990, Konigsberg and Blangero 1993, Relethford 1994, Relethford and Crawford 1995, Wescott and Jantz 1999, Tatarek and Sciulli 2000) use the R-matrix method of Harpending and Jenkins (1973), which treats phenotypic data as genetic markers to provide estimates of genetic distances between populations, as well as estimates of F_ST, or subpopulation heterozygosity, in relation to total population heterozygosity.

In this paper, we investigated whether differences in cranial morphology between Mednyi Island Arctic foxes and Bering Island Arctic foxes could be a result of the severe population bottleneck suffered by the Mednyi population in the 1970s. Firstly, linear measurements were taken in Arctic fox crania collected in the wild before and after the bottleneck, on both Mednyi and Bering Islands, and differences among them were assessed with pairwise analysis of variance (ANOVA). Secondly, we test the effects of the bottleneck on the variability of the population by detecting deviations from expectations under mutation-drift equilibrium, following Relethford and Blangero (1990). Thirdly, we applied a model to study whether or not the observed divergence between populations could be explained by random genetic drift alone.

Materials and methods

A sample of 117 dry crania of adult Arctic foxes were used in this study. Crania are housed in the Zoological Museum of M.V. Lomonosov Moscow State University, and were collected in the wild on Mednyi Island before the bottleneck event (29 specimens, including 16 males, 10 females and three of unknown sex) and after it (24 specimens, including nine males, six females and nine of unknown sex), and on Bering Island before (32 specimens, including 17 males and 15 females) and after (32 specimens, including 15 males and 17 females). Seven linear measurements of the cranium (Figure 1; raw data provided as supplementary information, Table S1) were taken with a Sylvac digital calliper (Sylvac, Crissier, Switzerland) to an accuracy of 0.1 mm. The measurements are defined as follows: overall condylobasal length (CBL), defined as the distance between the base of incisors at the midline and the posterior curve of the condyles; absolute rostrum length (RL), defined as the distance between base of incisors at the midline and ectorbitale; braincase length (BL), defined as the ectorbitale-akrokranion distance; nasal bone length (NASL), defined as the maximum length of the nasal bone along the lateral edge; maximum braincase width (BRCW), defined as the maximum width of the cranium at the parietal bones; postorbital width (PORW), defined as the width of the cranium at the postorbital constriction; braincase height (BRCH), defined as the distance between the maximum bend of parietal bones (excluding the sagittal crest) and the base of the cranium between the condyles. These measurements were chosen because they have been used to describe cranial morphology accurately, and are also likely to be integrated with one another (Nanova 2009).

Figure 1:

Superior and lateral views of an Arctic fox cranium showing the measurements taken for this study.

CBL, overall condylobasal length; RL, absolute rostrum length; BL, braincase length; NASL, nasal bone length; BRCW, maximum braincase width; PORW, postorbital width; BRCH, braincase height.

The Commander Islands Arctic foxes can be considered as a single population with random mating, subdivided into two subpopulations, Mednyi and Bering. For the purpose of this study, each island subpopulation is further divided into before and after the bottleneck, giving a total of four subpopulations with size estimates taken from the literature: Mednyi population before the bottleneck, 1000 individuals (Geptner and Naumov 1967); Mednyi population after the bottleneck, 90 individuals (Goltsman et al. 2005); Bering before the bottleneck, 3000 individuals (Geptner and Naumov 1967); and Bering after the bottleneck, 600 individuals (Ryazanov 2002). The pairs of subpopulations compared are designated throughout this paper as “Mednyi Before vs. Mednyi After”, “Bering Before vs. Bering After”, “Mednyi Before vs. Bering Before” and “Mednyi After vs. Bering After”.

An ANOVA was performed for each pair of subpopulations, for each measurement. Since previous research (Nanova and Prôa 2017) had found that the differences between males and females are consistent across Arctic fox populations of Mednyi and Bering Islands, no sex distinction was made during analyses; sex information was nevertheless retrieved and is available with the raw data as supplementary information, Table S1.

Originally devised for genetic data (Harpending and Jenkins 1973), the F_ST bottleneck detection method was developed to accommodate phenotypic data (Relethford and Blangero 1990, Relethford 1991, Relethford and Crawford 1995, 2013, Relethford et al. 1997). The method treats morphometric data as genetic markers and computes an R-matrix of genetic distances between subpopulations from the pooled within-group phenotypic variance-covariance matrices (the within-group diversity, average of the individual diversities for individual traits), scaled on the heritability of a trait (h²), and outputting an F_ST, which is a measure of deviation from the Hardy-Weinberg equilibrium (heterozygosity) for subdivided populations with random mating (Gillespie 2004, p. 119). Whether with genetic or phenotypic data, in a population with random mating, F_ST converges to 1 when neither mutations nor gene flow are present (genetic drift removes the within-group variance, leaving only the between-group variance); when mutations are absent but gene flow is present, F_ST will approach 0 (the between-group variance will be removed by gene flow, leaving only the within-group variance); in the presence of both mutations and gene flow, F_ST takes values between 0 and 1 (neither within-group nor between-group variance are removed) (Relethford et al. 1997, Lynch and Walsh 1998). F_ST can therefore be used to gage the effect of a population bottleneck: basically, the higher the F_ST, the greater the effect of random genetic drift in a population.

In the absence of heritability estimates for Arctic foxes, the phenotypic covariance matrix was assumed to be proportional to the additive genetic covariance matrix (Cheverud 1988) and was scaled to h²=1. The genetic distances are often computed with heritability estimates of h²=1 in order to produce minimum estimates of F_ST (González-José et al. 2001, 2007, Betti et al. 2009, 2010, Delgado 2017). Nevertheless, exploratory runs of other values of h² were tested (not shown) and they did not affect the results. Computation of distances was performed in the software RMET, version 5.0 (Relethford 2003), freely available at http://employees.oneonta.edu/relethjh/programs/. In total, 65 Relethford-Blangero analyses were ran; all the results are in Table 2. A summary table of each analysis, from which the F_ST is computed [using unbiased F_ST, recommended by Relethford (Relethford 1991)], is available as supplementary information (Tables S3–S67). Based on the F_ST properties, we developed general predictions:

1. Prediction 1: Mednyi Before and Mednyi After are the same population, so we expect F_ST to be very low when only these two subpopulations are considered, and for all the traits. This is because if they are the same populations they are related by genealogy (effectively gene flow).

2. Prediction 2: Same for Bering Before vs. Bering After as in Prediction 1.

3. Prediction 3: Comparing Mednyi Before vs. Bering Before should yield much higher F_ST, because they have been drifting apart from a common ancestor, and because we know there is no gene flow between them (the two islands are isolated).

4. Prediction 4: Same for Mednyi After vs. Bering After as in Prediction 3.

5. Prediction 5: If the bottleneck did have an effect on the cranial traits (i.e. that drift randomly removed within-group variation), then we expect F_ST in Predictions 1 to be higher than F_ST in Predictions 2.

6. Prediction 6: If the bottleneck did have an effect on the cranial traits (i.e. that drift randomly removed within-group variation), then we expect F_ST in Predictions 4 to be higher than F_ST in Predictions 3.

Bottleneck tests were run on all measurements, but more specific predictions were made, concerning three traits, RL, PORW and BRCH, following noteworthy ANOVA results:

1. Prediction 7: RL and PORW are statistically significantly different in Mednyi Before vs. Mednyi After, therefore, if they are different due to the effect of random drift, we expect the F_ST to be higher than for other traits, meaning that these traits were more affected by the bottleneck event.

2. Prediction 8: BRCH, which is not statistically significantly different in any subpopulation comparison, is expected to have a low F_ST in all cases, compared to other traits, meaning that it is a trait not affected by the bottleneck event.

3. Prediction 9: Removing RL and PORW from the analysis should decrease F_ST.

4. Prediction 10: Removing BRCH from the analysis should increase F_ST.

Finally, to complement the bottleneck analyses, the within-group and between-group variance-covariance matrices were compared to determine whether or not the observed diversity could be explained by random genetic drift alone. Following the claim (Cheverud 1988) that, in contemporary populations, the within-group variance-covariance matrix is often proportional to the genetic variance-covariance matrix, the latter can be substituted by the former. Comparing the between-group variance-covariance matrix and the within-group variance-covariance matrix (as a surrogate of the average genetic variance-covariance matrix) was accomplished by using the method of Ackermann and Cheverud (Ackermann and Cheverud 2002). The null hypothesis of divergence by random genetic drift alone is rejected if the slope of the regression (β) deviates significantly from 1. When using a significance level of α=0.05, it is expected that a true null hypothesis has a 5% chance of being rejected (a type I error). This test was proved to be robust in falsifying the underlying assumptions (Prôa et al. 2013), and has been used consistently (Ackermann and Cheverud 2004; Prôa 2016, Prôa and Matos 2017). The analyses were run in R (R Development Core Team 2018), using code available in the literature (Prôa et al. 2013) and modified.

Results

The results of ANOVA (Table 1) show that differences in RL (F=4.63, p=0.03618) and PORW (F=6.75, p=0.01223) are statistically significant between Mednyi Before and Mednyi After. No statistically significant differences are found between Bering Before and Bering After in any of the traits, with the exception of PORW (F=15.33, p=0.00023). Differences in CBL and PORW are not statistically significant between Mednyi Before and Bering Before, but are statistically significant between Mednyi After and Bering After (CBL, F=9.17, p=0.00377; PORW, F=44.22, p≤0.0001); differences in RL, BL, and NASL are statistically significant both between Mednyi Before and Bering Before (RL, F=14.38, p=0.00035; BL, F=7.26, p=0.00917; NASL, F=86.16, p≤0.0001), and between Mednyi After and Bering After (RL, F=10.46, p=0.00208; BL, F=13.56, p=0.00054; NASL, F=181.17, p≤0.0001). Differences in BRCW are statistically significant between Mednyi Before and Bering Before (F=6.87, p=0.01113), but are not statistically significant between Mednyi After and Bering After. The only trait that showed no statistical significance in any of the comparisons was BRCH.

Values of computed F_ST with standard errors are presented in Table 2. F_ST of Mednyi Before vs. Mednyi After and Bering Before vs. Bering After are very low for all traits, as predicted in Predictions 1 and 2. F_ST in Bering Before vs. Bering After is so low for some traits (CBL, BL, NASL, BRCW, BRCH) that it is close to zero. When considering Predictions 3 and 4, F_ST is relatively high for some traits, but not for all: NASL and PORW yield an F_ST higher than other traits, but it is only when some traits are removed from the analysis that the value rises; even so, it still tends to remain below 0.4, far from converging to 1. As predicted by Prediction 5, F_ST of Mednyi Before vs. Mednyi After is higher than F_ST of Bering Before vs. Bering After for every trait, except PORW. As predicted by Prediction 6, F_ST of Mednyi After vs. Bering After is higher than F_ST of Mednyi Before vs. Bering Before for every trait, except RL and BRCW. The four specific Predictions 7–10 were correct: RL and PORW yielded F_ST higher than for other traits in Mednyi Before vs. Mednyi After; BRCH yielded F_ST low in all cases, compared to other traits; removing RL and PORW from the analysis decreased F_ST; removing BRCH from the analysis increased F_ST.

The result of the drift test was that divergence by genetic drift alone could not be excluded as an explanation for the differences between subpopulations. In none of the analysis the slope of regression deviated statistically significantly from 1 (Mednyi Before vs. Mednyi After, β=1.3600, p=0.3151; Bering Before vs. Bering After, β=0.2902, p=0.1885; Mednyi Before vs. Bering Before, β=1.0416, p=0.4695; Mednyi After vs. Bering After, β=0.4712, p=0.1926).

Discussion

The results showed that differences in cranial morphology between the Mednyi population and the Bering population after the 1970s could be attributed to the population bottleneck caused by an outbreak of mange epizootic.

As predicted, comparing each island subpopulation before the bottleneck with the respective subpopulation after the bottleneck (Predictions 1 and 2) yielded a very low F_ST, demonstrating low between-group variance. Dividing each island population into two subpopulations, one before and one after the bottleneck, allowed us to estimate the effect of random drift in each island population, caused by the bottleneck, which seems to have been relatively low in Bering Island Arctic foxes, in spite of the population size reduction from 2000 to 600 individuals. Five out of the seven traits measured showed no effect of drift at all in Bering Arctic foxes, while in the Mednyi Island Arctic foxes only two out of seven traits showed no effect of drift. We therefore conclude that differences in cranial morphology on Mednyi Island Arctic foxes could be attributed to the severe population bottleneck in the 1970s.

Comparing the Mednyi Island population before the bottleneck with the Bering Island population before the bottleneck (Prediction 3) was expected to return high F_ST for all traits, due to lack of Arctic fox migration between the two islands, and because they have been diverging from a common ancestral for a considerable amount of time. Yet, though relatively high for NASL and PORW, F_ST remained substantially low, and that may be because, even though there is no gene flow between the two populations, they still belong to the same species, and 10,000 years of isolation may not have been enough divergence time for a signature to be found in the cranial morphology.

We interpret the higher F_ST in Mednyi Before vs. Mednyi After (Prediction 1), compared to Bering Before vs. Bering After (Prediction 2), as a bottleneck signature. A considerably larger effect of random drift took place in the Mednyi population, sufficient for it to be detected by this method. This was not true, however, for one trait, PORW, which suggests this trait was much less subject to the effect of drift in Mednyi Island Arctic foxes, or that it was in fact more affected by it in Bering Island Arctic foxes. Either way, it seems to be a trait that has consistently changed in both the Mednyi population and the Bering population through time (as seen by the ANOVA results), and we interpret that as a consequence of going through the bottleneck.

That F_ST is higher when comparing Mednyi After vs. Bering After (Prediction 4), than in Mednyi Before vs. Bering Before (Prediction 3), indicates a pattern of divergence consistent with random genetic drift. Because the population was larger, Bering Island Arctic foxes lost less genes at random (20% of individuals survived). Mednyi Island Arctic foxes lost many more genes at random (only 9% of individuals survived), including probably many genes that it shared with Bering, which means they now “drifted apart” even more than they had before the bottleneck event. RL and BRCW escape this pattern, and this may be due to their being subjected to stronger selective pressure than other traits: RL contraction in Mednyi Arctic foxes was interpreted as necessary to keep the bite force large at a larger gape angle to catch large prey (Nanova and Prôa 2017). Indeed, we started the study presented here precisely because we had found differences in RL between Mednyi Before and Mednyi After which could be attributed to the bottleneck effect (Nanova and Prôa 2017). Extending the question to other traits which are likely to have been affected by the bottleneck if RL was (the cranium is an integrated whole), in this study we found that Bering Island Arctic foxes were also affected by the bottleneck (there are differences not in RL, but in PORW in Bering Before vs. Bering After), but in a different way from Mednyi Island Arctic foxes. Traits that show differences before and after the bottleneck event are not the same in both populations, which shows the randomness of the event that contributed to the recent divergence in morphology in these populations. RL, the trait that varied most in previous studies, is in fact the one least affected by the bottleneck.

Predictions 7–10, concerning individual traits, were correct. Statistically significant differences in RL and PORW in Mednyi Before vs. Mednyi After can be attributed to the effects of the bottleneck event. Differences in BRCH were not found statistically significant in Mednyi Before vs. Mednyi After, meaning a higher degree of similarity between populations unlikely to be due to a strong effect of random drift, shown by the low F_ST. Removing from the analysis the traits likely affected by drift, i.e. RL and PORW, resulting in a lower F_ST, means that the ensemble of all other traits suffered less the effects of the bottleneck. Likewise, removing BRCH, resulting in a higher F_ST, means that the ensemble of all other traits were more affected by the bottleneck. No predictions were made on traits whose differences were statistically significant between Mednyi and Bering Island Arctic foxes, both before and after the bottleneck, because it would be difficult to detect any effect of drift on traits which had already diverged between island populations.

A null hypothesis of divergence by genetic drift alone could not be excluded, meaning a strong effect of random drift in these populations, though these results could also be explained by the fact that all the subpopulations considered here are too closely related for non-random factors to be detected by the method of Ackermann and Cheverud (Ackermann and Cheverud 2002).

Application of bottleneck tests requires that population declines have a high probability of being detected and that bottlenecks are not regularly inferred for stable populations (Peery et al. 2012). Indeed, bottleneck tests have failed to detect well-known population collapses in Scandinavian lynx [Lynx lynx L., 1758; (Spong and Hellborg 2002)], California sea otters [Enhydra lutris nereis Merriam, 1904; (Aguilar et al. 2008)] and Amur tigers [Panthera tigris altaica Temminck, 1844; (Henry et al. 2009)]. In our study the population decline was observed in the field and is firmly established in the literature. Our findings of strong effect of random genetic drift within the same population, as measured from the cranial morphology, can be a consequence of the bottleneck effect. Bottleneck detection is critical for the interpretation of the demographic history of the endangered Mednyi Island Arctic fox, and the next step could be a Bayesian approach to the demographic history of these Arctic foxes.