|Home | About | Journals | Submit | Contact Us | Français|
Habitat fragmentation is a ubiquitous by-product of human activities that can alter the genetic structure of natural populations, with potentially deleterious effects on population persistence and evolutionary potential. When habitat fragmentation results in the subdivision of a population, random genetic drift then leads to the erosion of genetic diversity from within the resulting subpopulations and greater genetic divergence among them. Theoretical and simulation analyses predict that these two main genetic effects of fragmentation, greater differentiation among resulting subpopulations and reduced genetic diversity within them, will proceed at very different rates. Despite important implications for the interpretation of genetic data from fragmented populations, empirical evidence for this phenomenon has been lacking. In this analysis, we carry out an empirical study in populations of an alpine meadow-dwelling butterfly, which have become fragmented by increasing forest cover over five decades. We show that genetic differentiation among subpopulations (GST) is most highly correlated with contemporary forest cover, while genetic diversity within subpopulations (expected heterozygosity) is better correlated with the spatial pattern of forest cover 40 years in the past. Thus, where habitat fragmentation has occurred in recent decades, genetic differentiation among subpopulations can be near equilibrium while contemporary measures of within subpopulation diversity may substantially overestimate the equilibrium values that will eventually be attained.
Habitat fragmentation is recognized as a leading cause of loss of biological diversity worldwide (Davies et al. 2001), and can affect the genetic structure of plant and animal populations (Templeton et al. 1990; Young et al. 1996; Hale et al. 2001), with potentially negative consequences for population persistence (Saccheri et al. 1998). When a population is subject to habitat fragmentation, the increased isolation of resulting subpopulations has two principal genetic effects: differentiation among the subpopulations increases while genetic diversity within them declines (Templeton et al. 1990; Young et al. 1996). Numerous studies have documented both these responses to decreased habitat patch connectivity (Soulé 1972; Hänfling & Brandl 1998; VanDongen et al. 1998; Williams et al. 2003; Andersen et al. 2004). Analytical and simulation studies predict that the two principal genetic responses to population subdivision will often proceed at dramatically different rates (Latter 1973; Crow & Aoki 1984; Varvio et al. 1986; Wang 1997). This result has important implications for the interpretation of population genetic data from recently fragmented populations (Varvio et al. 1986). However, despite significant interest in the genetic effects of habitat fragmentation and habitat patch isolation, this phenomenon has not previously been demonstrated in nature and has largely been ignored by empirical researchers.
Theoretical examinations of the dynamics of genetic variation in subdivided populations indicate that differentiation among subpopulations and diversity within subpopulations can approach equilibrium values along different trajectories and at very different rates (Latter 1973; Crow & Aoki 1984; Varvio et al. 1986; Wang 1997). When a formerly panmictic population becomes subdivided, random genetic drift will lead over time to both genetic differentiation among the resulting subpopulations and a loss of heterozygosity from within the resulting subpopulations (Hartl & Clark 1989). As genetic drift tends to drive allele frequencies within subpopulations towards fixation, genetic variability is removed from subpopulations and their average heterozygosity is reduced, while allele frequencies in different subpopulations tend to become, on average, more divergent. If some degree of gene flow among the subpopulations is maintained, it will act to counter the effects of drift and prevent complete fixation of alleles in subpopulations. In that case, both the differentiation among subpopulations and the heterozygosity within subpopulations will reach equilibrium values that are determined by the balance between genetic drift and gene flow.
Crow & Aoki (1984) modelled a subdivided population with migration, mutation, drift and a finite number of subpopulations, and examined the temporal dynamics of changes in subpopulation differentiation and local gene identity (homozygosity, which is inversely related to local heterozygosity). They showed both analytically and numerically that in such a subdivided population with migration among the local demes, the degree of differentiation among the subpopulations (GST) reaches close to the equilibrium value very rapidly, while homozygosity within subpopulations changes very slowly and takes many more generations to approach its equilibrium value. This result held for both island and stepping-stone models of population structure and with a wide range of parameter values. Varvio et al. (1986) performed numerical simulations using a finite island model of population structure and examined the temporal dynamics of the average gene diversity within subpopulations (HS) and differentiation among subpopulations (GST). Specifically, they examined the time required for a value differing by only 5% from the equilibrium value to be reached. They also demonstrated that GST approached equilibrium far more rapidly than did HS and concluded that in natural populations GST was the variable more likely to be close to its equilibrium value.
Here, we demonstrate different rates of response of the within- and among-subpopulation components of genetic variation to recent habitat fragmentation in an alpine butterfly. We have been studying populations of the meadow-dwelling butterfly, Parnassius smintheus Doubleday, in a series of alpine meadows atop three ridges in the eastern Canadian Rocky Mountains. These meadows are bordered by even-aged stands of lodgepole pine (Pinus contorta), subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii) and have been divided into 22 study sites representing habitat patches that are separated to varying degrees by intervening forest (figure 1). Because of a large fire that occurred in 1938 (Roland et al. 2000), the open, non-forested habitat on these ridges was much more extensive in the past. Since that time, re-establishment of forests and fragmentation of the ridge-top meadows has occurred (figure 2). Based on mark–recapture analysis of movements of adult P. smintheus indicating that forests act as strong barriers to dispersal, rates of movement along these ridges are predicted to have been substantially higher before re-establishment of forests occurred (Roland et al. 2000). Here, we present data that suggest that in response to the recent re-forestation of the ridge-tops, and ensuing reduction in movement rates of P. smintheus, the among- and within-subpopulation components of genetic variation in P. smintheus have changed at different rates. Specifically, we show that currently observed genetic differentiation among our study sites is best correlated with contemporary forest cover in the landscape, while currently observed genetic diversity within our study sites is better correlated with forest cover almost 40 years in the past.
Our study area comprised 22 study sites, each of which represents a patch of potential habitat for P. smintheus. These study sites were located along three ridges in the Kananaskis region of Alberta, Canada: Lusk Ridge (51°00′ N, 114°58′ W), Jumpingpound Ridge (50°57′ N, 114°53′ W) and Cox Hill (50°58′ N, 114°54′ W). Because movements of P. smintheus occur primarily along ridge-tops (Roland et al. 2000), and because the sites in our study area are distributed end-to-end along these ridges, with respect to butterfly movement, the sample sites are essentially distributed along a one-dimensional axis (figure 1). We have previously shown genetic differentiation among these sites and a significant pattern of isolation by distance along the ridge-top axis (Keyghobadi et al. 1999). Furthermore, the one-dimensional structure of the ridge-top habitats allows us to readily divide the distance between any two sample sites into components that are over meadow versus forest habitat. We have also previously examined the relative influence of each of these types of intervening habitat on the degree of genetic differentiation between pairs of study sites, and found that the amount of contemporary intervening forest explained a greater proportion of the variation in pairwise genetic distances (Keyghobadi et al. 1999). This result was in concordance with results from mark–recapture studies indicating that individuals are less likely to fly over forested habitat than over open, meadow habitat (Roland et al. 2000).
We were interested in examining whether the genetic diversity of P. smintheus subpopulations within habitat patches, measured as expected heterozygosity (He), was correlated with the contemporary degree of connectivity of each patch to all others in the study area, using a measure of connectivity that incorporates detailed information on both movement behaviour of the focal species, and landscape structure. For each study site (i.e. habitat patch) defined in our study, we measured its connectivity to the other patches using a modified version of the metapopulation measure of patch connectivity defined by Hanski (1994). Our measure indicates the potential number of immigrants, and therefore gene flow, into a habitat patch from all other possible patches in the study area. This measure of connectivity assumes that the rate of immigration into the target patch is a negative exponential function of the distance from a source patch, with different rates of exponential decline if the intervening distance is over forest versus meadow habitat, to account for the differential effects of these landscape elements on individual movement (Roland et al. 2000). The rate of immigration into a target patch from a particular source patch also depends on the area of the source patch, scaled to account for nonlinear effects of patch area on population size and emigration rate. Patch connectivity, as we have measured it, is the inverse of patch isolation, another commonly used concept in metapopulation and landscape ecology.
We measured the connectivity of each habitat patch to all other habitat patches in the study area as
where is the distance through forest between patches i and j, is the distance through meadow between patches i and j, and αf and αm are parameters that scale the rate of movement to distance for each of the two habitat types, respectively. Aj is the area of patch j. ξp is a parameter that scales local population size in a patch to the area of the patch, and ξem scales the probability of emigration to patch area. Summation is over all other patches in the study area.
Estimates of αf, αm and ξem were obtained using the virtual migration model (VMM; Hanski et al. 2000), which provides maximum-likelihood estimates of parameters of migration and survival from mark–recapture data. We used a modified version of the VMM that allows for estimation of more than one α value (Matter et al. 2004). Mark–recapture data from these same study sites were obtained in 1995 and 1996 for an investigation of the effects of landscape structure and population size on movement of individuals (Roland et al. 2000). Using the mark–recapture data from 1995, which fit the VMM better than did the data from 1996, the estimated parameter values were: αf=5.81, αm=1.34, and ξem=0.0 (Matter et al. 2004). Also from the mark–recapture data, indices of population size within each patch were obtained (Roland et al. 2000). The value of ξp was estimated to be 0.84 by fitting a power function to the relationship between patch area and the average of the population size indices from 1995 and 1996. Population sizes fluctuate yearly and we expected that the average of the 2 years would be a slightly better indicator of long-term population size.
Note that although similar in form, our measure of connectivity is different from the measure defined in the VMM (Hanski et al. 2000). That measure indicates the tendency for a particular habitat patch to send migrants out to all other patches in the study area, the opposite of what our measure indicates. Our measure of connectivity also ignores the possibility of gene flow from outside the study area. This assumption seems reasonable as the closest potential sources of immigrants are over 5km away, and including those sites in the measures of connectivity would have had little effect on the calculated values.
For a measure of connectivity that did not account for differential effects of intervening forest and meadow habitat, we used
where dij is the total distance between patches i and j and α=2.98. This α is a weighted average of the αf and αm values, weighted by the amount of forest and meadow habitat, respectively, separating all pairs of patches in the study area.
As an analogous measure of connectivity between each possible pair of patches in the study area we used
where Aij is the sum of the areas of patches i and j.
Using black and white aerial photographs (1:40000), the areas (ha) of all study patches were estimated as in Roland et al. (2000). Three distance measures (km) were also quantified for all pairs of patches: (i) total distance, which comprise (ii) the distance through forest and (iii) the distance through meadow. All distances were measured between the centroids of butterfly observations within each patch and along the ridges, since movements of P. smintheus are largely constrained along ridge-tops (Roland et al. 2000). We obtained landscape data from aerial photographs from both 1993 (only 2 or 3 years before samples for genetic analysis were collected) and from 1952 (the oldest available high-quality aerial photographs of the area).
Tissue samples from adult P. smintheus were collected from 17 of the 22 study sites (figure 1) in 1995 and 1996. Sample sizes are given in table 1. Samples were either small (approximately 0.15cm2) wing clippings or, for females and males respectively, approximately 25mg of thoracic or abdominal tissue. Genomic DNA was isolated using the QIAamp tissue extraction kit (QIAGEN). Each individual was typed at seven microsatellite loci (Ps50, Ps81, Ps85, Ps76, Ps163, Ps165 and Ps262) as described previously (Keyghobadi et al. 1999, 2002). Products of polymerase chain reaction amplification were electrophoresed and detected using an ABI 373A Automated Sequencer, and their sizes were determined and analysed using Genescan and Genotyper software (ABI).
For the samples analysed here, conformity of genotypes at the seven microsatellite loci to Hardy–Weinberg proportions have previously been tested (Keyghobadi et al. 1999, 2002). Because of excess homozygosity and evidence of null alleles at all loci (Keyghobadi et al. 1999, 2002), for each locus at each site, the frequency of the null allele was estimated and the frequencies of all other alleles were simultaneously re-estimated using the estimation–maximization algorithm, as in Keyghobadi et al. (1999). These adjusted allele frequencies were used to calculate unbiased estimates of expected heterozygosity for each locus within each patch (Nei & Roychoudhury 1974). Within each patch, expected heterozygosities were averaged over loci to give an overall indicator of genetic diversity, the average expected heterozygosity (He). Expected heterozygosity is a preferred measure of genetic diversity that is less sample-size dependent than alternative measures such as the observed number of alleles (Soulé & Yang 1973; Nei 1987). The significance of correlations between He and patch connectivity were assessed using Kendall's rank correlation (Sokal & Rohlf 1995).
Adjusted allele frequencies were also used to estimate genetic differentiation, GST (Nei 1973), for all pairs of patches. The significance of correlations between GST and pairwise patch connectivity were assessed using the Mantel test (Mantel 1967).
Genetic variation within patches was generally high but varied among patches, with the average number of alleles at seven microsatellite loci ranging from 6.7 to 11.3, and the average expected heterozygosity (He) ranging from 0.643 to 0.766 (table 1). The number of individuals sampled was positively correlated with the average number of alleles (r2=0.80, p=0.0001), but not with He (r2=0.22, p=0.93), confirming our decision to use He as a measure of within-patch genetic variability in our analyses and indicating that our results are not biased by differences in sample size among study sites.
Patch connectivity in 1993, using a measure that accounts for different rates of butterfly movement through forested and meadow habitats, was significantly positively correlated with within-patch genetic diversity (He; r2=0.29, p=0.04; figure 3a). Higher connectivity allows for a higher rate of incoming gene flow, which counters the effects of local genetic drift and has a positive effect on levels of genetic variation. This result confirms our concerns about the potential impact of increasing habitat patch isolation to reduce within-patch genetic variation in this species (Keyghobadi et al. 1999) and agrees with other studies that have found a significant, positive correlation between various measures of habitat patch connectivity and genetic variation (Soulé 1972; Saura et al. 1973; Descimon & Napolitano 1993; Schmitt et al. 1995; Hänfling & Brandl 1998; Shapcott 2000; Vucetich et al. 2001). Our result is unique however in demonstrating this correlation at a very fine spatial scale and in using a sophisticated measure of connectivity that incorporates information on both detailed landscape features and movement behaviour of the species. Our study area is dominated by the large, forested valley separating Lusk Ridge (sites C, D1, d2, and E) from Jumpingpound Ridge and Cox Hill (all other sites), and this may have an inordinate influence on our results. However, removal of the Lusk Ridge sites from the analysis, including from the calculations of connectivities of other sites, had little effect on the correlation between He and patch connectivity (r2=0.30, p=0.04). This is not surprising given that connectivity is largely determined by the immediately neighbouring patches, and patches located across the large valley have very little effect on each other's connectivity.
In a subdivided population, the genetic diversity within any subpopulation may be affected not only by the rate of gene immigration from other subpopulations but also by the size of the local subpopulation. We have, however, found no significant correlations between expected heterozygosity and Roland et al.'s (2000) indices of population size within patches from 1995 and 1996 (Electronic Appendix). Furthermore, we have not found significant correlations between expected heterozygosity and the harmonic mean of the population size indices for these years, nor with indirect indicators of population size such as habitat patch area or host plant abundance (Electronic Appendix). It may be that in a series of subpopulations that are connected by relatively high levels of gene flow, as in our study sites, habitat patch connectivity and the rate of gene immigration have a stronger influence on local genetic diversity than does local population size.
Quite surprisingly, we found that using only the total distance between patches, not accounting for reduced rates of movement through forest in our measure of patch connectivity for 1993, improved the observed correlation between He and connectivity (r2=0.39, p=0.008). The difference between this correlation coefficient and that obtained when reduced movement through forest is accounted for is not significant (ts=0.32, p=0.75; using the jackknife procedure as described in Sokal & Rohlf 1995), but is clearly opposite in direction to our expectation. We would expect a measure of connectivity that does not account for differential effects of intervening forest and meadow habitat on movement to be a poorer measure and, therefore, less strongly correlated with within-patch genetic diversity. Our result initially appears paradoxical given that both levels of genetic differentiation among patches and rates of movement are clearly related to the relative amounts of intervening forest versus meadow habitat (Roland et al. 2000; Keyghobadi et al. 1999).
This apparent paradox can, however, be explained by a different rate of response of among-patch genetic differentiation compared with within-patch genetic diversity to recent habitat fragmentation. As forest cover in our study area has increased since the 1938 fire, rates of movement, and therefore gene flow, among local populations are predicted to have declined significantly (Roland et al. 2000). However, as predicted by theory, genetic differentiation among patches would have responded more rapidly to this fragmentation than would the loss of diversity within patches. In our study sites, using population size estimates from mark–recapture data (Roland et al. 2000) and assuming that effective population size is one-tenth to one-half of census population size, the product of the number of occupied patches and the average effective population size of each patch is estimated to be of the order of hundreds or thousands. Also, observed rates of movements among patches (Roland et al. 2000) and a modest degree of genetic differentiation among patches (Keyghobadi et al. 1999) suggest that on average, the migration rate (Nm) is certainly greater than 0.1. Based on simulations of subdivided populations using a finite island model with migration (Varvio et al. 1986), for a group of populations of the size we have studied, and with similar migration rates, we would expect GST to approach its equilibrium value in tens of generations, while within-patch heterozygosity (or its complement, homozygosity) should approach its equilibrium in hundreds of generations, following fragmentation. Since P. smintheus is univoltine, these expected times to reach equilibrium translate directly to tens and hundreds of years, respectively. Such a difference in rate of response to recent habitat fragmentation would explain why genetic differentiation among patches reflects the contemporary patterns of forest cover in the landscape, while levels of genetic diversity within patches do not.
We were able to test this hypothesis by comparing the within- and among-patch components of genetic variation to both contemporary and past landscape structure using aerial photographs of the study area from 1993 (only 2 or 3 years before samples for genetic analysis were collected) and from 1952 (the oldest available high-quality aerial photographs of the area; figure 2).
When 1952 landscape data (intervening forest and meadow distances, and areas of patches) were used to calculate patch connectivity, the correlation with within-patch genetic diversity (He) was stronger than when we used 1993 landscape data (figure 3). Thus, He is better correlated with patch connectivity in the landscape 40 years ago than it is with patch connectivity today. In contrast, genetic differentiation between patches, measured as pairwise GST (table 2), was better correlated with an analogous measure of pairwise patch connectivity when 1993 landscape data were used than when 1952 landscape data were used (figure 4). These observed differences in correlation coefficients between the 1993 and 1952 data are not statistically significant (comparison of r values between 1952 and 1993 for the correlation between patch connectivity and He yields ts=0.13, p=0.89; for the correlation between pairwise connectivity and GST, ts=0.49, p=0.62). This is perhaps not surprising given the relatively small number of data points (17 sites) and the high error associated with estimates of genetic diversity and genetic differentiation, especially when they are based on a small number of independent loci. However, the trends we have observed in the comparisons of correlations between the 1993 and 1952 landscape data suggest that the among-patch component of genetic variation is better correlated with contemporary landscape structure, while the within-patch component is better correlated with the structure of a past landscape. The discordant responses of among- and within-patch genetic variation to contemporary landscape structure appear to be owing to a much more rapid response of among-patch differentiation to recent habitat fragmentation and population subdivision, as predicted by theoretical models.
In the scattergram showing correlation of GST and log pairwise connectivity from 1952 (figure 4b), there are two clouds of data points. The smaller cloud with log connectivity measures of approximately −12 to −8 represents pairs of patches from different sides of the forested valley that separates Lusk Ridge from Jumpingpound Ridge (figure 1). The larger cloud with log connectivity measures ranging from approximately −5 to 2 represents pairs of patches from the same side of the valley. In 1952 there was little forest cover on the ridges, so connectivity between patches on the same side of the valley was high, with a narrow range of values. The valley separating Lusk and Jumpingpound ridges was still forested, however, so connectivity between patches separated by the valley was much lower. Thus, pairwise connectivity values in 1952 were bimodal. In contrast, in 1993, when there was much more forest cover on the ridges, there was greater overlap between connectivity values for patches on the same versus opposite sides of the valley, and the pairwise connectivity values for 1993 are continuous (figure 4a). To remove any potential effect of the bimodal distribution of 1952 pairwise connectivity values on the analysis, we also examined correlations between log pairwise connectivity and GST for pairs of patches on Jumpingpound Ridge and Cox Hill only (i.e. we excluded all sites on Lusk Ridge from the analysis). The correlation between GST and log pairwise connectivity was again higher in 1993 (r2=0.15, p=0.01) than in 1952 (r2=0.12, p=0.03), confirming our previous results. As demonstrated previously, the analyses of correlations between heterozygosity and the connectivity of individual patches (figure 3) are not affected by the presence of the forested valley separating Lusk and Jumpingpound ridges because the connectivity of a given patch is mostly determined by the immediately surrounding patches.
Habitat fragmentation results in population subdivision and reduces connectivity among resulting subpopulations, but may also reduce habitat area and, hence, total census population size (Fahrig 2003). Theoretical models of subdivided populations account for the reduced effective population size and greater genetic drift that accompany the reduction of connectivity in a formerly panmictic population. However, any additional reduction in census population size owing to loss of habitat is generally not accounted for in the models, and could potentially confound our results by reducing even further the effective size of the population. The simulations of Varvio et al. (1986) suggest that a reduction in the size of local subpopulations (as might occur with habitat fragmentation) would reduce the time required for both He and GST to approach equilibrium. However, the number of generations required to reach equilibrium would remain an order of magnitude greater for He than for GST, even with a 1000 fold reduction in the size of local subpopulations. Therefore, any reductions in census population sizes that have accompanied the recent habitat fragmentation in our study area are only likely to have precipitated an even more rapid approach to equilibrium of among-patch differentiation (GST), with within-patch diversity (He) continuing to lag far behind.
The number of species experiencing habitat fragmentation and population isolation is rapidly increasing and the genetic consequences of habitat fragmentation remain a concern for conservation biologists (Templeton et al. 1990; Young et al. 1996). There exists an extensive literature on the genetic effects of habitat fragmentation, and numerous empirical studies have examined the genetic structure of subdivided populations within this context (Cunningham & Moritz 1998; VanDongen et al. 1998; Young et al. 1999; Mossman & Waser 2001; Vucetich et al. 2001; Williams et al. 2003; Wilson & Provan 2003; Andersen et al. 2004; Hänfling et al. 2004; Schad et al. 2004; Van Rossum et al. 2004). In genetic studies of such populations, the importance of considering historical factors and non-equilibrium dynamics in general has been stressed (Varvio et al. 1986; Cunningham & Moritz 1998; Kinnison et al. 2002; Schmitt & Hewitt 2004). However, the different rates of approach to equilibrium of the within- and among-patch components of genetic variation are an aspect of non-equilibrium dynamics that have largely been ignored. Our results indicate that this phenomenon can be detected in natural populations and should be considered when attempting to assess the genetic impacts of recent habitat fragmentation. The rapid approach of genetic divergence to equilibrium can be seen in cases where significant genetic differentiation among subpopulations is detected following very recent habitat fragmentation or other perturbations (Gerlach & Musolf 2000; Kinnison et al. 2002). Conversely, the slow approach of heterozygosity to equilibrium in subdivided populations may explain situations where heterozygosity is not correlated with subpopulation isolation (Young et al. 1999; Jäggi et al. 2000; Shapcott 2000; Mossman & Waser 2001; Hinten et al. 2003; Hänfling et al. 2004; Llorens et al. 2004), and suggests caution when interpreting such results. Lack of a contemporary correlation between patch isolation and heterozygosity may, in many cases, reflect a very slow approach to equilibrium of within-patch genetic diversity.
For fieldwork, we thank Norine Ambrose, John Brzustowski, Mark Caldwell, Sue Cotterill, Marian Forrester, Sherri Fownes, Fiona Johnson, Kris Sabourin, Chris Schmidt, Lynette Scott, Kirsty Ward, Arliss Winship and the Kananaskis Field Stations (University of Calgary). We thank Graeme Taylor, Dina Fonseca, James F. Crow and two anonymous reviewers for very helpful comments on the manuscript. Atte Moilanen kindly provided the programme to run the virtual migration model. Financial support was provided by NSERC, Parks Canada, Alberta Sports Recreation Parks and Wildlife Foundation, Friends of Banff National Park, NSF, and a Challenge Grant in Biodiversity (Alberta Conservation Association). This manuscript was partially prepared while N.K. was a fellow at the Smithsonian Institution.
As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.