Temporal and spatial instability in neutral and adaptive (MHC) genetic variation in marginal salmon populations

The role of marginal populations for the long-term maintenance of species’ genetic diversity and evolutionary potential is particularly timely in view of the range shifts caused by climate change. The Centre-Periphery hypothesis predicts that marginal populations should bear reduced genetic diversity and have low evolutionary potential. We analysed temporal stability at neutral microsatellite and adaptive MHC genetic variation over five decades in four marginal Atlantic salmon populations located at the southern limit of the species’ distribution with a complicated demographic history, which includes stocking with foreign and native salmon for at least 2 decades. We found a temporal increase in neutral genetic variation, as well as temporal instability in population structuring, highlighting the importance of temporal analyses in studies that examine the genetic diversity of peripheral populations at the margins of the species’ range, particularly in face of climate change.

The importance of marginal populations for the long-term maintenance of species' genetic diversity and evolutionary potential has long been discussed 1,2 but it is now particularly timely in view of range shifts caused by climate change 3,4 . According to the Centre-Periphery hypothesis, marginal populations inhabit unstable and poorly connected habitats, and may be expected to harbour less genetic variation and have lower evolutionary potential than those at the centre of the distribution 5,6 . However, although genetic diversity within populations seems to decline on average from the centre of the distribution to the periphery 6 , there is no conclusive evidence that geographically, historically or climatically marginal populations display lower average fitness 7 . Understanding the dynamics of populations at the species' range limits, including patterns of extinction and recolonization and their ability to adapt to environmental variation, is key to predict their responses to climate change 8 . Critically, studies looking at genetic diversity in central versus marginal populations have largely focused on contemporary patterns of genetic diversity, using primarily neutral markers 6 , which may not fully reflect the adaptive potential of populations 9 .
The genes of the major histocompatibility complex (MHC), which are some of the most studied in relation to adaptive genetic variation 10,11 , are useful markers to use in combination with neutral markers to reconstruct not only the genetic diversity but also the adaptive potential of marginal populations. MHC genes are central to immunity as they encode for proteins that present pathogen-derived antigens to T-cells, initiating the adaptive immune response 12 . MHC genes are amongst the most polymorphic genes in vertebrates and also some of the best studied 13 and the variation in the MHC residues that bind antigens from pathogens is thought to be maintained by balancing selection driven by pathogens 14,15 and also influenced by mate choice 16,17 . Evidence of selection on the MHC genes has been identified in many species as heterozygote advantage 18,19 , association of individual MHC alleles and/or genotypes with susceptibility to specific pathogens 20,21 , rare-allele advantage 22 and changes in allele frequencies under experimental infections 23 . Comparisons between MHC loci, or markers linked to them, and neutral markers can be used to infer differences in the relative levels of neutral and adaptive variation within and among populations, which can be variable not only across closely related species 24,25 but also within species depending on the spatial scale of analysis 26 .

Results
Microsatellite variability and population structuring. Individuals captured by anglers in the rivers Asón, Deva, Nansa and Pas (Northern Spain; Fig. 1) between 1948 and 2002 were genotyped at 13 putatively neutral microsatellite DNA markers. Deviations from HWE were only observed for loci CTAX (13 samples) and Sssp2210 (8 samples) following strict Bonferroni correction for multiple tests ( Table 1). The overall results of the analyses did not change by excluding these two microsatellites (data not shown) and we opted for including them. A total of 196 alleles were observed across neutral microsatellites for the whole sample. BOTTLENECK results indicated that allelic frequency distributions did not depart from the expected L-shaped distribution. The number of alleles ranged from five at locus SsaD486 to 26 at locus SsaD144b. Significant correlations in allele frequencies between adjacent temporal samples for all comparisons could suggest stability in allele frequencies, but none of the correlations was significant after applying strict Bonferroni correction for multiple tests (Supplementary table S1). Heterozygosity (Ho) increased significantly over time in the rivers Pas (Mann-Kendall trend test s = 8 P = 0.041), Nansa (s = 6 P = 0.042) and Pas (s = 10, P = 0.008) whilst no significant temporal change in Ho was observed in the Deva (s = 6 P = 0.117). Equally, allelic richness increased temporally in the rivers Asón (s = 8 P = 0.042), Nansa (s = 6 P = 0.042) and Pas (s = 8 P = 0.042), but not in the river Deva (s = 6 P = 0.117). Following population analysis by STRUCTURE, the estimated optimal number of genetic groups was K = 6 ( Fig. 2). The results indicated that the rivers Pas and Nansa, and to a lesser extent the river Ason, suffered a drastic change in population structure post-80 s, such that the genetic composition of these rivers in 2002 is rather different than that observed in 1950 and 1960. Similar results were obtained when the rivers were analysed individually (Supplementary figure S1). Genetic distance (D A ) between temporal samples of the same rivers (0.057 to 0.177) were of a similar order to genetic distances between river samples (0.060 to 0.215). The NJ-phenogram, although with low statistical support, suggested that samples from the rivers Ason and Deva tended to group by river and not by decade, whereas samples from the rivers Nansa and Pas were intermingled, with a tendency to associate by decade instead of river (Fig. 3). AMOVA results also indicated significant temporal heterogeneity within the Asón and Pas samples ( Table 2). Among river genetic variation was significant for each temporal sample, excluding the 1990 s samples, though there seemed to be a decrease in magnitude of F ST over time.  Considering  P-value = 0.531). Observed heterozygosities ranged from 0.17-1.00 across all samples (Table 1) Table S1). This indicates stability of allele frequencies over time for both Sasa-UBA-3′ UTR and Sasa-DAA-3′ UTR loci. Allele distributions mostly overlapped among the four populations, with some very low frequency alleles being only represented in one or two of the rivers (Supplementary Figure S2). AMOVA analysis indicated significant genetic structuring at the class I locus at all temporal periods tested, whereas for the class II locus, the 1980s and 1990s were significantly differentiated as well as samples from 2002 and the 1960s ( Table 3). The phylogenetic tree for the class I-linked marker indicated a relationship among rivers very similar to that for the neutral microsatellites, clustering the samples from rivers Asón and Deva according to river of origin, whereas the samples from rivers Nansa and Pas were largely intermingled within a cluster. In contrast, the class II marker, showed no structuring of samples based on river of origin in any of the rivers.    Discussion Peripheral (marginal) populations tend to be genetically and morphologically distinct as a consequence of their isolation and typically smaller size, and are considered particularly valuable because they can help preserve the evolutionary potential of the species 2 . Atlantic salmon populations in northern Spain represent peripheral populations at the southern limit of the species' range; these have been in decline since the 1960's and are now classified as endangered 51 . However, despite inhabiting the margins of the species' range and having small effective population sizes, these populations display levels of genetic diversity comparable to those reported for larger populations at the center of the distribution 50 . Northern Iberian rivers are thought to have been a refugium for Atlantic salmon during the last glacial maximum, and it is possible that this is the reason why these populations appear to harbour higher than expected ancestral mitochondrial DNA variation compared to more northerly European populations 28 . Additionally, stocking from different sources carried out in the 80s and in the 90s could have also contributed to the temporal differentiation of these populations 31 . Between 1972 (when stocking records start) and the 90s, these rivers (initially the Rivers Ason and Pas and then extending to the River Nansa and to lesser extent the River Deva) were stocked with high densities of eyed ova (200,000-300,000 annually) and fry (90,000-120,000 annually) mainly from Scotland and Iceland. Stoking from the 90s was carried out from native sources, and primarily from the river of origin during the last years 48,51 . Our results from neutral microsatellites indicate that there has been a temporal increase in genetic diversity (heterozygosity and allelic richness) in three of the four rivers over a 50 year period, but also some temporal maintenance of genetic identity in the river Deva. In contrast, increases in neutral genetic diversity in the rivers Ason, Nansa and Pas, coupled with the strong changes in their genetic background from the 80s and a temporal decrease in genetic structuring suggest that their genetic composition could have been affected by foreign stocking, as previously indicated using mtDNA 31 . We found no conclusive evidence of selection in the MHC-linked markers, apart from Sasa-UBA-3′ UTR, that together with three neutral microsatellites deviated from neutrality in samples from the 60s and 80s, suggesting the parallelism between the class I marker and the rest of the microsatellites. Results from neutral markers largely mirrored those of the class I MHC-linked marker (Sasa-UBA-3′ UTR) but not those of Sasa-DAA-3′ UTR (class II). This is perhaps not surprising given the differences in response to selection previously observed between both markers 46 . In this case, Sasa-DAA-3′ UTR did not indicate any clustering of samples, by river or decade. Such a pattern of variation could reflect adaptation to local conditions in these marginal populations, an scenario that might be expected given the homing behavior of Atlantic salmon and their tendency to form locally adapted populations 37 , but also genetic drift due to low effective population size. Recently introduced salmonid populations in Chile suggested that MHC class II functional diversity of invasive populations has decreased over time, in contrast to diversity at neutral markers which has remained very high 52,53 as a consequence of admixture 54 . Therefore, it is possible that, even if some neutral diversity has remained high in some Iberian salmon rivers, perhaps as consequence of foreign stocking, diversity at non-neutral markers may have been eroded over time due to geographical differences in selection 55 and to adaptation to local conditions (e.g. parasites) 56,57 . Our current results, in combination with previous studies on the same populations, indicate, despite their low effective population size 50 , these salmon harbour high neutral genetic diversity, atypical in marginal populations, highlighting the importance of the demographic history for the maintenance of the genetic diversity. This is particularly relevant in view of the predictions of the consequences of climate change for salmonids, i.e. movement of the thermal niche of salmon towards north as well as decreased production and population extinction in the southern range of species 58 . Our study highlights the importance of a adopting not only a spatial but also a temporal approach, considering both neutral as well as adaptive markers, in studies that examine changes in genetic diversity of peripheral populations at the margins of the species' range.

Methods
Origin of the samples and DNA extraction. Adipose fins from dead adult Atlantic salmon, captured by anglers in the rivers Asón, Deva, Nansa and Pas (Northern Spain; Fig. 1) were collected in 2002 and stored in 95% ethanol at 4 °C prior to genetic analysis. Dried scales from the same rivers collected from adult fish caught by anglers since 1948 were also included in the analyses. Due to their limited availability, historical scales were pooled across four decades following Ciborowski et al. 31 : 1948-1957, 1960-1963, 1980-1989, 1990-1996. No scale samples were available for any river from the 1970s, or for the river Nansa pre-1960. Therefore, 19 groups of samples, stratified by decade and river were generated for analysis (Table 1).
Total DNA was extracted using the Promega ™ Wizard SV 96 Genomic DNA Purification System.
Manufacturer's protocols were adhered to for modern adipose fin samples, but for historical scales we increased the incubation time during the elution steps to five minutes and decreased the elution volume to 80-100 μ l. Between one and three historical scale extractions were carried out in a dedicated ancient DNA laboratory, physically separated from PCR procedures. A blank control was extracted concurrently and subsequently amplified in PCR reactions. All eluted DNA was stored at − 20 °C.
PCR amplification and microsatellite genotyping. All  were run on a 3100 ABI Prism capillary sequencer using the Genescan-500 LIZ size standard. Alleles were scored using Genemapper V3.5 software (Applied Biosystems) and genotypes were manually checked.
Data analysis. Individuals with fewer than eight successfully genotyped loci were discarded from analysis (final sample sizes in Table 1). All loci were tested for conformity with Hardy-Weinberg equilibrium using the randomization test implemented in GENEPOP 65 and were also tested for neutrality using the Ewens-Watterson test 66 in ARLEQUIN v3 67 . Rates of allelic dropout (ADO) and false alleles (FA) were estimated according to Broquet and Petit (2004) 68 . Allelic richness (Ar) was calculated using FSTAT 69 . Statistical significance of temporal trends was tested using the Mann-Kendall trend test 70 74 was used to test how many genetic populations were represented by all individuals caught in each of the four rivers. We followed the methodology outlined in 74 . First, we constructed phylograms for all individuals from each river based on individual distance matrices calculated with the program POPULATIONS 75 using an allele sharing distance (ASD) method 76 to visualise whether there was any clustering of individuals into discrete population units. Following this, all individuals from each river were modelled in STRUCTURE. The program was run applying the admixture model, as this model was likely to be closer to the true nature of the history of these populations compared with a non-admixture model. The parameters of the simulations were burn-in length of 50,000 iterations; 100,000 MCMC repetitions; testing for K (the number of populations) between 2 and 8 over 10 repeated simulations. We estimated the correct value of K using the Evanno  Table 3. Temporal and spatial genetic structuring based on F ST in 4 marginal populations of Atlantic salmon (rivers Ason, Nansa, Pas and Deva in Northern Spain), estimated using MHC-linked markers. *p < 0.05 **p < 0.01 ***p < 0.001.  method 77 as implemented in STRUCTURE HARVESTER (http://users.soe.ucsc.edu/~dearl/software/struct_harvest/). We then used CLUMPP 78 and DISTRUCT 79 to summarise and represent the results. Spatial and temporal structuring was also analysed using AMOVA as implemented in ARLEQUIN and POPULATIONS was used to generate a consensus unrooted neighbour-joining tree (10,000 bootstrapped replicated) of the samples based Nei's D A distance 80 , which was visualized using TREEVIEW 81 . Sasa-UBA-3′ UTR and Sasa-DAA-3′ UTR data were analysed separately (as in 64 ). GENEPOP on the web 82 was used to estimate observed and expected heterozygosities (H o and H e ) for each locus in each sample. Allelic richness was calculated at each locus for each sample using FSTAT version 2.9.3 83 . Statistical significance of the temporal trends of genetic diversity (heterozygosity and allelic richness) was tested using the Mann-Kendall trend test 70 implemented in PAST 71 .
All markers, neutral and MHC-linked microsatellites, were tested for neutrality using Lositan 84,85 , under 50,000 simulations, estimated neutral mean F ST , infinite alleles mutation model, 99% confidence interval and false discovery rate of 0.1%. All populations were tested for recent bottlenecks using BOTTLENECK v.1.2.02 86 .