Whether and how habitat fragmentation and population size jointly affect adaptive genetic variation and adaptive population differentiation are largely unexplored. Owing to pronounced genetic drift, small, fragmented populations are thought to exhibit reduced adaptive genetic variation relative to large populations. Yet fragmentation is known to increase variability within and among habitats as population size decreases. Such variability might instead favour the maintenance of adaptive polymorphisms and/or generate more variability in adaptive differentiation at smaller population size. We investigated these alternative hypotheses by analysing coding-gene, single-nucleotide polymorphisms associated with different biological functions in fragmented brook trout populations of variable sizes. Putative adaptive differentiation was greater between small and large populations or among small populations than among large populations. These trends were stronger for genetic population size measures than demographic ones and were present despite pronounced drift in small populations. Our results suggest that fragmentation affects natural selection and that the changes elicited in the adaptive genetic composition and differentiation of fragmented populations vary with population size. By generating more variable evolutionary responses, the alteration of selective pressures during habitat fragmentation may affect future population persistence independently of, and perhaps long before, the effects of demographic and genetic stochasticity are manifest.
The study of habitat fragmentation has primarily focused on its ecological effects and classic genetic impacts to populations, such as inbreeding and genetic drift [1–4]. That habitat fragmentation might affect natural selection by altering habitat characteristics as fragment and population size decrease has rarely been considered [5–7]. Consequently, we lack a clear conceptual and theoretical framework for predicting how habitat fragmentation affects the adaptive genetic composition and differentiation of fragmented populations of varying sizes, and how this might affect subsequent responses to environmental change.
One possibility, the ‘Directional Hypothesis’, proposes that habitat characteristics shift in a consistent manner during the habitat fragmentation process [6,7], resulting in directional relationships between these characteristics, population size and the extent of adaptive genetic variation and differentiation. For example, fragmentation might decrease population size while concurrently increasing isolation and reducing habitat quality, leading to greater organismal stress [8,9]. Hence, under this hypothesis, small, fragmented populations are predicted to consistently exhibit reduced adaptive variation/differentiation relative to large populations due to the resulting combined effects of genetic drift, restricted gene flow and/or inbreeding [2,4,10]. However, the Directional Hypothesis might not always typify the nature of adaptive variation. Some empirical work indicates that natural selection can overcome genetic drift at very small population sizes in the wild [11,12]. Moreover, environmental change associated with increased fragmentation and reduced habitat size can result in the increased maintenance of adaptive polymorphisms at small population size through, for example, balancing selection .
A possible alternative to the Directional Hypothesis is the ‘Variable Hypothesis’. It posits that habitat characteristics and resulting natural selection pressures become increasingly variable as habitat fragment size and population size decrease [6,7]. This stems from the idea that the evolutionary effects of habitat fragmentation are contingent upon the starting conditions within habitat fragments. For example, habitat complexity is usually greater at larger spatial scales , and the fragments that small populations occupy might simply be random samples of the fragments that large populations occupy (, references therein). Some of these small population fragments will thus typify the habitat heterogeneity of the larger population fragments while some will be more homogeneous. Therefore, under the Variable Hypothesis, the resulting extent of adaptive genetic variation and differentiation among populations is expected to become more variable with increasing fragmentation and decreasing population size.
These hypotheses provide a useful point of departure for an empirical exploration of the degree to which differences in the size of fragmented populations are associated with differences in adaptive variation and differentiation. Until recently, a lack of technological advancement and high-throughput capability in molecular ecology precluded such investigations [15–17]. Yet, as suggested by Wood et al. , if selection pressures within habitat fragments differ or vary more in intensity or direction among small than large populations, such assessments might improve the ability to predict population responses to future environmental change and the setting of species conservation priorities.
We examined patterns of putative adaptive variation and differentiation using detailed analyses of 164 single-nucleotide polymorphisms (SNPs) located in coding genes, with some also linked to QTLs for different phenotypic traits, in 14, varying-sized stream populations of brook trout (Salvelinus fontinalis). We specifically assessed the relationships between population size and: (i) genetic variation; (ii) extent of neutral population differentiation between population pairs and (iii) presence of signatures of diversifying or balancing selection between population pairs. We predicted that the Variable Hypothesis would better explain patterns of probable adaptive genetic variation and differentiation among these populations. Notably, previous work found no consistent difference in the characteristics of habitat fragments that small and large populations occupied, but rather more divergent habitat characteristics and greater spatial habitat variability among small populations . Importantly, both the adult census population size (N) and the effective number of breeders (Nb) were considered in our analyses, the latter being an analogue of effective population size (Ne) but for a single breeding event . Although some research has assumed a correspondence between N and Nb , Nb/N ratios can vary widely among populations of closely related species (e.g. ). It is also the genetic parameter Nb—not the ecological parameter N—that selection is associated with and that influences the magnitude of adaptive genetic variation and differentiation .
2. Material and methods
(a) Study site
Cape Race (CR) is a small, coastal barren region in southeastern Newfoundland, Canada, characterized by a parallel series of low-order streams harbouring pristine, resident populations of brook trout (electronic supplementary material, figure S1). CR streams are small (0.27–8.10 km range in axial length) and they can be comprehensively sampled, resulting in comparatively accurate and precise estimates of N and Nb. Fish residing in different streams form genetically distinct populations  that are estimated to have diverged from a common ancestor during the late-Wisconsinan glaciation (10–12 000 ybp; ). CR populations also exhibit considerable differences in life histories, probably because of population-specific differences in selective pressures following habitat fragmentation [22,23].
(b) Adult census population size (N) and effective number of breeders (Nb)
Estimates of population size for each study population were obtained from Wood et al.  with one exception (Hermitage Pond Brook), based on multi-year N and Nb estimates for most populations (electronic supplementary material, table S1). Either the Schnabel  or Petersen  method was applied to estimate annual N, whereas Nb was estimated using LDNe software  and based on data at 13 neutral microsatellite loci. (Additional details on calculating N and Nb for CR populations are available from .) Extremely low polymorphism was expressed in the small-population Hermitage Pond Brook at the 13 microsatellites, the neutral SNPs documented in this study, and at an additional 19 screened microsatellites (DJ Fraser 2014, unpublished results); this prevented us from providing a reliable estimate of Nb. An Nb estimate of 11 was therefore assumed for our analyses, based on a conversion from N of this population, and using the best-fit relationship between reliable estimates of Nb and N for all CR trout populations: inverse(Nb/N) = 5.9567 + 0.0115 × N (see Appendix I in Wood et al. ).
(c) Population genotyping
Tissue samples from a total of 446 individuals from 14 CR populations  were genotyped at coding gene SNPs previously developed for brook trout ([26,27]; see the electronic supplementary material, figure S1 and table S1, for names and locations of different streams, and numbers of individuals sampled per population). Briefly, these samples were obtained in the summer of 2011 (2012 for Hermitage Pond Brook) as adipose fin clips and stored in 95% ethanol until DNA was extracted, using a modified phenol–chloroform protocol. Samples within populations were obtained from individuals of the same age, unambiguously identifiable as being in their first year of life, based on body size [22,23]. Individuals were randomly sampled from 15 to 30 locations within each stream, using effort-standardized surveys (3-min electrofishing) conducted at each 50 or 100 m intervals originating at the stream mouth, depending on stream length; a subset of these samples was randomly selected for genotyping at 237 SNPs in this study. These SNPs have previously been positioned on a genetic map, tested for association with quantitative trait loci (QTL) at many physiological and growth traits, and annotated when feasible [26,27]. Details of SNP development, validation and sequencing at the Genome Quebec Innovation Center (McGill University, Montreal, QC, Canada) are provided in [26,27]. Three monomorphic SNPs were excluded from all analyses, as were four SNPs having heterozygosities exceeding 50% but for which the second homozygous genotype was missing. Four individuals with call rates of less than 85% across SNPs were also removed. Of 237 SNPs screened in our 446 individuals, 164 SNPs were polymorphic in 442 individuals and retained for analyses after the application of our quality criteria.
(d) Relationship between population size and genetic variation
We firstly assessed whether a relationship existed between heterozygosity, the simplest measure of within-population genetic variation, and population size. The analysis accounted for potential non-independence of loci among populations and within individuals , which is often not undertaken.
We applied generalized linear mixed models (GLMMs), using heterozygosity as a binary variable (homozygous = 0, heterozygous = 1), for 71 406 combinations of locus-individual combinations as the response (1.5% of combinations were missing data), the fixed continuous covariates of N or Nb (each mean-centred) and random intercept effects for the 14 populations, 164 loci and 442 individuals (nested within populations). The overall relationship between heterozygosity and population size was estimated across 164 random slope effects of loci (i.e. loci-by-population size effects) that we included in the models. We also estimated the covariance between loci intercepts and loci slope effects. GLMMs were fit with a binomial residual distribution and a logit link function under Laplace approximation, using the ‘glmer’ function of the R-package LME4 , and executed in R v. 2.15.3 . Significances of fixed and random terms were tested using deviance-based likelihood ratio tests (LRT) between nested models with and without a respective term.
(e) Adaptive genetic differentiation among populations
We conducted a genome scan for each population pair to determine the probable extent of adaptive genetic differentiation among populations (i.e. total genome scans = 91), based on the numbers of SNPs exhibiting signatures of diversifying or balancing selection between each population pair (see below). Genome scans were conducted by implementing the commonly used software LOSITAN . We estimated loci that deviated significantly higher (diversifying selection) or lower (balancing selection) from neutral genetic differentiation (FST) as predicted by the average relationship of heterozygosity and FST. We first calculated an average ‘neutral’ FST and forced this average FST by repeating runs with each of 50 000 simulations. The expected number of populations in these simulations was 14. The simulation p-values were corrected for the false discovery rate (FDR) after . Only loci with q-values of less than 0.001 were regarded as showing significant signatures of selection. This approach corresponded to an overall FDR of 0.002 as we estimated q-values separately at both ends of the probability distribution.
(f) Relationship between population size and neutral population differentiation
We analysed the relationship between population size and neutral genetic differentiation as estimated from pairwise genome scans. This provided ‘neutral’ FST estimates between each population pair, once SNPs exhibiting signatures of selection in each pairwise comparison were removed from the dataset. Our pairwise estimation process assumed that each pairwise FST estimate might be affected by one or two population sizes, their interaction and/or by the geographical distance separating the population pair (the shortest distance between each pair's stream mouths). Preliminary analyses indicated a positive but non-significant relationship between pairwise FST (range: 0.018–0.828) and geographical distance (range: 0.03–15.29 km). However, to remove any geographical influence from our data, we firstly obtained model residuals from the FST–distance relationship. We then added to these residuals the FST value corresponding to the average pairwise geographical distance, to obtain estimates that were standardized with the average geographical distance. We then fit linear mixed models (LMMs) for the response of ‘remaining variation in FST’ with the two fixed continuous covariates of population size (either N or Nb) for both populations contributing to each respective FST value, plus the interaction between both population sizes. As all data obtained from pairwise comparisons were correlated by two different populations simultaneously, we fit random intercept effects for the two populations in all models. All LMMs were fit with a Gaussian error distribution using the ‘lmer’ function of the LME4 R-package. Significances of the fixed terms were tested using F-tests with denominator degrees of freedom (ddf) adjusted according to .
(g) Relationship between population size and adaptive genetic differentiation
We also analysed the relationship between population size and putative adaptive genetic differentiation as obtained from pairwise genome scans using GLMMs. We regarded adaptive genetic differentiation as the detection probability of an investigated locus to be under selection (either balancing or diversifying selection) out of the total number of loci analysed for a given population pair. Accordingly, for each pairwise genome scan, one population size, two population sizes, their interaction and/or the geographical distance separating the population pair (the shortest distance between each pair's stream mouths) might influence the probability to detect loci under selection. Not all of our loci are independent as some are linked [26,27] and this might inflate or deflate detection probabilities, although effects of linkage on genome scans remain to be investigated . To be able to account for linkage group identity of loci in our follow up GLMM analyses, we coded our results as a binary response of each locus in each pairwise comparison (either balancing or diversifying selection; no selection detected = 0, selection detected = 1).
We analysed the binary responses with the three fixed continuous predictors of population sizes (N, Nb) for population pairs, their interaction and geographical distance, fitted with a binomial error distribution and a logit link function. As all binary data were obtained from pairwise comparisons, they were correlated by two different populations simultaneously. To account for this, random effects for each of the two populations and for population pairs were fitted in all models (population pairs effects accounted for non-independence across loci within a given comparison). Lastly, we included linkage group identity random effects to account for non-independence among loci within and across comparisons. Random locus effects were confounded with linkage group effects for many loci, and therefore were not fitted. This was because no linkage was detected in previous studies for 44 out of 164 loci [26,27], and thus we treated these loci as single locus linkage groups. Significance testing of model terms was conducted as for heterozygosity models.
(a) Relationship between population size and genetic variation
After accounting for variation caused by non-independence in the data, we detected a significant, positive relationship between multi-locus genetic diversity and both population size measures (N or Nb) (table 1 and figure 1). The largest variance for heterozygosity was present among loci (61% and 62% of the total variance for Nb and N, respectively). The estimated variance among loci for slopes of heterozygosity with population size was also significant (table 1). Furthermore, a negative correlation was exhibited between random locus effects for intercepts and slopes (ρ = −0.42, −0.33 for Nb, N; table 1). That is, loci with low heterozygosity at overall average population size exhibited a larger change in heterozygosity with population size change relative to loci with high heterozygosity. Less but significant variance was also present among populations (12% and 14% of the total variance for heterozygosity for Nb and N, respectively) and among individuals within populations (both models; 0.5% of the total variance for heterozygosity).
(b) Adaptive genetic differentiation among populations
The percentage of outlier SNPs differentiating population pairs ranged between 0% and 8%, with a total of 31 and 145 loci showing significant signatures of balancing and diversifying selection, respectively, in at least one pairwise genome scan. Of these, 18 (balancing) and 13 (diversifying) were uniquely exhibited in a single pairwise genome scan. The highest numbers of loci under selection in any pairwise genome scans were 4 (balancing) or 18 (diversifying).
(c) Relationship between population size and neutral genetic differentiation
Depending on the population pair, an average of 113.5 ± 21.3 (1 s.d.) SNPs were used to estimate ‘neutral’ FST after removing outliers (and ‘monomorphic’ SNPs corresponding to that pair). Only non-significant relationships were obtained between variation in pairwise FST based on ‘neutral’ SNPs (and corrected for pairwise geographical distance) and both population size measures (electronic supplementary material, table S2). However, a trend of decreasing FST with increasing population size was observed when using Nb (F1/12 = 3.0, p = 0.111; electronic supplementary material, figure S2) but less so when using N (F1/12 = 1.8, p = 0.207; electronic supplementary material, figure S2), as the rank of N value for a population did not necessarily concur with the same rank for the Nb value (electronic supplementary material, figure S3).
(d) Relationship between population size and adaptive genetic differentiation
We found significant relationships between both population size measures (N, Nb) and the proportion of loci under balancing selection, as well as under diversifying selection for Nb, whereas the geographical distance separating population pairs had only a non-significant effect on the extent of adaptive genetic differentiation (table 2). Accounting for the potential non-independence of loci on the same linkage groups, these significant relationships (or the non-significant trend between N and diversifying selection) were characterized by an influence of the interaction of population sizes from both populations on the probability to detect loci under selection (table 2). Relationships were also consistently stronger for estimates of Nb than N. For balancing selection, the highest probability was modelled for combinations of smallest with largest population sizes, an intermediate probability for small with small populations and the lowest probability for large with large populations (figure 2). For diversifying selection using Nb, the relationship (or trend using N) was similar as for balancing selection, but the differences in probabilities between small with small populations and small with large populations were nearly absent here and effects were weaker (figure 2).
Our study revealed evidence of putative adaptive variation in brook trout over a broad range of population sizes for vertebrate species (N = 73–8416; Nb = 11–249) and at a fine geographical scale (0.03–15.29 km). Probable adaptive genetic differentiation was greater between small and large populations or among small populations than among large populations, based on signatures of balancing and/or diversifying selection at SNPs linked to phenotypic traits. These patterns were stronger for estimates of Nb than N, an important finding given that theoretically it is through Nb that selection should influence the magnitude of adaptive genetic differentiation [20,35]. What we did not find was evidence of reduced signatures of selection in small populations relative to large populations, despite that overall genetic variation (multi-locus heterozygosity) and neutral differentiation (pairwise neutral FST) were consistent with pronounced genetic drift in small populations.
Our results strengthen insights from recent work which indicate that selective pressures change during the habitat fragmentation process as fragment and population size decrease [5–7], even once population sizes become very small. These changes lead to (i) more variability (selective pressures vary more among small populations than among large populations); (ii) more ‘buffering’ (selective pressures favour genetic polymorphism underlying some traits in at least some small populations, as shown by evidence for balancing selection in comparisons primarily involving small populations); and to a lesser extent (iii) more diversification (selective pressures in small and large populations can be systematically different).
Our results therefore have elements supporting both the Directional and Variable hypotheses outlined above, so they do not entirely mirror what is known of the physical habitat of our study populations that can affect fitness in brook trout (, references therein). For example, Wood et al.  found that small trout populations were more often associated with more spatially divergent physical habitats than large populations, and we found that this does in fact result in more divergent selective pressures among small than large populations. Contrastingly, we also found some evidence that small- and large-population selective pressures can be consistently different from each other, suggesting that other underlying environmental selective pressures might be involved in driving the observed adaptive divergence. Many SNPs used in this study have been linked to growth, reproductive and stress response traits in brook trout [26,27]. In fact, of the 17 SNPs linked to QTL that were also polymorphic in CR populations, 16 were outliers in our study in two to 12 pairwise population comparisons. Additional and/or cryptic selective pressures might therefore relate to: (i) temporal variability in the physical habitats among small and large populations that affect growth or survival, specifically greater temporal variability within small populations , which might explain the observed prominence of balancing selection in small populations; (ii) diverse, seasonal aspects of intraspecific competition for accessing breeding areas or overwintering habitat, both of which are limited in some CR streams and may be associated with body size [22,23,36]; and/or (iii) interspecific interactions with pathogens or parasites, and less so with competitors or predators (brook trout are the only fish species in most CR streams and avian/mammalian predation is virtually non-existent; ). In sum, interactions between all biotic/abiotic factors likely drive the composition and extent of adaptive differentiation among natural, fragmented populations of varying size.
We further suggest that the nature of adaptive variation and differentiation in CR brook trout, a product of natural habitat fragmentation since the Last Pleistocene deglaciation, might not substantially differ from the human-induced fragmentation that many species are being subjected to today. The primary difference between natural and human-induced habitat fragmentation perhaps lies more in the rate at which environmental selective pressures within habitat fragments are altered, than in how these pressures are altered . Indeed, conditions likely vary within landscapes with the end result depending on the initial conditions, regardless of whether fragmentation is natural or human-induced. Isolated populations reduced in size by the relatively slower, incremental environmental changes of natural fragmentation might be more likely to adapt and persist in the long term, but the collective phylogenetic and geological evidence suggests that fragmentation for many CR trout populations arose rapidly (; T Burdon and DJ Fraser 2014, unpublished data).
(a) Possible caveats and assumptions
We used the proportion of SNPs exhibiting signatures of selection to characterize the probable extent of adaptive genetic differentiation between populations of varying size. Although the brook trout genetic map into which this study's SNPs were embedded was sufficiently dense to cover almost all chromosomes [26,27], our metric should be considered heuristic for it assumed that all SNPs were ‘equal’ and did not covary in shaping adaptive differentiation . The use of population pairwise genome scans permitted us to investigate the extent to which, within a ‘matrix’ of populations of varying size, signatures of adaptive genetic variation were most evident. As these pairwise comparisons were therefore not strictly independent, two population random effects and linkage group identity effects were fit in all statistical models, and each population was represented proportionately equally across the 91 genome scans.
The overall percentage of SNPs exhibiting signatures of selection in our study was high, most likely due to the number of populations compared pairwise. Indeed, the range within population pairs (0–8%) was otherwise consistent with that documented for other species  and for brook trout, based on the same SNPs [38,39]. Furthermore, a high level of selection was not unexpected with these SNPs given (i) their association with traits potentially underlying fitness in brook trout [26,27] and (ii) that CR trout populations are known to differ extensively in physical environmental characteristics  and at several phenotypic traits that have been shown to proffer adaptive value [22,23]. Moreover, it is unlikely that hierarchical population structure might have increased the rate of false-positive signatures of selection given the recent common ancestry of our study populations, and that population isolation appears to have happened rather uniformly (T Burdon and DJ Fraser 2014, unpublished data).
Finally, although we adopted a stringent FDR for detecting outlier SNPs, we cannot rule out the possibility that a few (especially small) study populations might have had a history of severe bottlenecks, which can increase the generation of false-positive signatures of selection in genome scans, especially for diversifying versus balancing selection . On the other hand, more recent simulations suggest that high population differentiation can result in a higher rate of false-negatives , yet our study still found greater evidence for selection signatures in pairwise comparisons involving small populations.
(b) Evolutionary and conservation implications
Comprehensive studies on the adaptive genetic composition and differentiation of fragmented populations of known, varying sizes are rare. In our naturally fragmented trout populations originating from a common ancestor, signatures of putative adaptive genetic differentiation were most evident among small and large populations or among small populations. Such findings lend support to the hypothesis that conditions do exist in nature wherein the evolutionary trajectories of very small populations are still very much affected by natural selection, in addition to drift [11,12,42]. They also support the proposition that a fuller understanding of adaptation will come from considering how intra- or interspecific interactions might be altered by habitat changes that are directly or indirectly linked to population size, or how selective pressures might shift temporally as populations expand or decline in size. Similarly, our results suggest that balancing selection at some traits can act in nature at small population size and hence potentially buffer, to some degree, against the loss of genetic diversity at loci of high evolutionary importance [10,13].
The observed differences in selective pressures between small and large populations and the increase in their variability at smaller population size also have potentially important conservation ramifications in the face of growing, worldwide habitat fragmentation of natural populations. These results suggest that habitat fragmentation generates both directional and increasingly diversified changes to selective pressures as population size is reduced. Thus, (i) large and small populations can represent distinct entities, and (ii) collectively, different small populations may in fact harbour unique variation that is adaptive in a wide range of circumstances. Such knowledge is beneficial for a more informed approach to setting biodiversity conservation programmes and priorities.
A last but perhaps most important conservation implication of our work stems from the evidence of more varying selective pressures in small than large populations. This finding supports a previous contention that responses to environmental change, and ultimately the probability of persistence, may become more variable as habitat is fragmented and population sizes are reduced . Indeed, such ‘evolutionary stochasticity’ (to distinguish it from demographic, environmental or genetic stochasticity) might in fact occur well before populations become sufficiently small to succumb to classic genetic effects such as genetic drift and inbreeding depression. Collectively, our work reiterates the need for further integration of natural selection into conservation biology theory [5,20,43].
Research undertaken in this study complies with the requirements of the Canadian Council on Animal Care (CCAC).
Funding was provided by NSERC (Natural Sciences and Engineering Research Council) Discovery grants to J.A.H. and D.J.F.
We thank M. Yates, J. Wood, S. Belmar-Lucero, L. Weir, A. Harbicht, M. Bonamy and A. Meli for assistance with field sampling, R. Carson for assistance with microsatellite genotyping and A. Belisle from Genome Quebec for assistance with SNP genotyping.
- Received February 12, 2014.
- Accepted July 2, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.