A Benchmark for Atlantic Salmon Conservation: Genetic Diversity and Structure in a Southern European Glacial Refuge before the Climate Changed

Abstract

Atlantic salmon Salmo salar supports highly valuable commercial and recreational fisheries in Europe, but its stocks are currently overexploited and threatened by climate change. Its southernmost populations (in northern Spain) play a key role in conserving the species’ original genetic diversity, which is endangered due to decades-long (1970s to 1990s) massive stocking with non-native stocks. Their decline is well documented, but the effect of stock transfer and conservation efforts is unclear. Nine microsatellite loci were amplified from archival samples (scales from 1958–1959) from eight Spanish rivers to analyse the species’ natural genetic dynamics before its decline started. Allelic richness was high in the historical populations (the 1950s) and above most contemporary estimates. Private alleles were found in most rivers, indicating high local uniqueness and relative isolation among river basins. Some alleles are regional markers since they are rare or absent from contemporary northern European populations. Effective population size suggested good conservation status, with higher values than those estimated for contemporary populations. Strong population structure and genetic differentiation between rivers were found, with limited gene flow, restricted to geographically close populations. Our estimates of historical genetic diversity and structure from southernmost salmon populations are a powerful benchmark to guide conservation programs.

1. Introduction

Preservation of genetic diversity not only is paramount for the long-term conservation of species [1,2] but also enhances the capacity of ecosystems to withstand, recover, and adapt to environmental perturbations [3,4,5]. However, anthropogenic pressures on natural populations have long modified the genetic structure and gene flow patterns of fish species worldwide. These pressures hinder current attempts to identify which gene pools/stocks are genetically differentiated and which must be prioritised for conservation.

Atlantic salmon (Salmo salar) is an anadromous species with a very high socioeconomic status in Europe as it supports valuable commercial and recreational fisheries. However, since the second half of the 20th century, the species experienced alarming reductions in population size and genetic diversity due to overexploitation, anthropogenic impacts in fresh waters (i.e., water and habitat quality degradation, connectivity loss, and population admixing with non-native domesticated strains), and possible ecosystem change with reduced feeding opportunities for Atlantic salmon in the North Atlantic Ocean occurring in 2005. Production was reduced after that with smaller, slower growing, and later maturing salmon [6,7]. Consequently, the species was assessed as vulnerable (VU A2ace) in Europe [8]. These impacts have been more severe in southern European populations, where effective population sizes have been reported to be rather low [9] and negatively affected by current climate change [10,11,12]. Apart from the climate-driven ecosystem regime shifts leading to poorer trophic conditions that have occurred in their marine feeding habitats [11,12], southern populations have experienced the most serious impacts upon their local environmental conditions (mainly temperature and river flow) of its entire distribution [7,11,13]. These southern European salmon populations also have the singularity of following a westerly migration route towards Greenland instead of following the easterly branch of the North Atlantic current into the Norwegian Sea [14,15]. Better knowledge of these marine migration routes can help understand the causes of the ongoing salmon population decline.

Spanish populations have been profusely studied in terms of their genetic variability, population structure, and response to different management and conservation measures [16]. These studies were typically based on samples collected from the 1980s onwards, which corresponds to a period of drastic reduction in population sizes. This collapse was driven by the concurrence of an abrupt acceleration in the anthropogenic warming trend and the warm phase of the Atlantic Multidecadal Oscillation, which caused a major regime shift in biophysical conditions throughout the North Atlantic salmon feeding grounds [11]. Massive stocking of non-native salmon was performed in Spanish rivers between the 1970s and 1990s following the populations’ collapse in an attempt to aid population recovery [17]. Salmon were mainly imported from Scotland and released (between 20,000 and 100,000) annually at similar densities in all rivers [12]. However, information regarding the genetic diversity and structure of Spanish populations before this time is scarce [12,18,19]. In addition, conclusions from previous studies are strikingly disparate. Some results indicate a strong loss of diversity and population structure, significant levels of introgression, and a general population homogenization across river basins [12,20,21]. In contrast, others show a low impact of admixture with non-native strains [22] and a reduction in the effective population size of native stocks, while maintaining relatively high levels of genetic variation [18]. Finally, others reveal different effects between western and eastern Spanish rivers [23].

In order to carefully assess the current conservation status of Spanish Atlantic salmon populations and to understand the effects of anthropogenic disturbances, climate-driven environmental changes, and conservation efforts, it would be desirable to have an overview of historical genetic diversity (prior to stocking and the severe impacts of climate change). This can provide a benchmark for restoration attempts [24]. To this end, salmon scales collected in 1958–1959 in eight river basins that cover most of its Spanish distribution were analysed for variation at nine microsatellite loci. Microsatellites can be superior markers to SNPs for studying genetic diversity in natural populations, so microsatellite DNA polymorphism analysis can help support the success of conservation and restoration projects as an economically advantageous research technique [25,26]. Microsatellites were chosen in the present study because they also allow comparison with the significant amount of data already available in the literature and databases.

Our main objectives were to (a) study the genetic diversity of each population; (b) infer connectivity, migration rates, and population structure across river basins; (c) obtain indicators of the historical conservation status of each population, such as effective population size and evidence of bottleneck events. This information will provide insight into the changes experienced by Spanish salmon populations in the last 70 years and a benchmark for comparisons with current populations in the context of ongoing climate change.

2. Materials and Methods

2.1. Study Area

A total of 374 archival samples of dried scales of Atlantic salmon from 8 northern Spanish rivers (Navia, Narcea, Sella, Deva-Cares, Nansa, Pas, Asón, and Bidasoa; Figure 1) were analysed. Archival samples were collected in 1958–1959 from returning salmon by the National Inland Fishing and Hunting Service of the Spanish Ministry of Agriculture (Madrid). Scales of salmon returning after two winters at sea were selected.

The Spanish rivers constitute the southernmost part of the Atlantic salmon distribution in Europe. In this area, the geomorphology and hydrogeology of the rivers are determined by the proximity of the Cantabrian Mountains to the sea, which produces very steep, V-shaped valleys that widen in their downstream sections. The main courses of the river basins are narrow and short, less than 100 km long, except for the River Navia, with a length of 159 km. In general, there is a dominance of fast-water habitats, and all are very well-preserved streams free from human activities causing water or habitat quality degradation.

However, the local hydroclimatic conditions in the study area have markedly changed over the last decades. Almodóvar et al. [11] detected a significant upward trend in regional temperature anomalies during 1950–2011, with an abrupt increase in the late 1980s. During this period, local climate warming concurred with a significant continuous increase in the frequency of low-flow events and a decrease in the minimum flow over seven or thirty consecutive days in the main rivers draining into the Cantabrian Sea and the Bay of Biscay. Such hydroclimatic changes were significantly linked to the abrupt decline observed in salmon captures in the study area since the 1970s, which was strongly accelerated from the late 1980s onwards [11]. In addition to the archival samples, reference samples of Atlantic salmon microsatellites for the eastern Atlantic (SALSEA) [27] were used for comparison.

2.2. Trends in Population Abundance

Atlantic salmon catch per unit effort (CPUE) data collected by the management authorities of Asturias, Cantabria, and Navarra in the study rivers for the last 73 years (1950–2022) were analysed. River catches were used as a proxy for abundance because this is the only measure of Atlantic salmon abundance consistently recorded in Spanish rivers over time. Until the early 2000s, the local authorities pooled salmon catches from the three rivers of Cantabria (Nansa, Pas, and Asón); so for consistency, we pooled their catches for the whole 73-year series. The salmon fishing season lasts from March to July, and the fishing regulations, number of fishing licences, and, thus, fishing effort have not changed significantly over the study period [11]. The nonparametric Mann–Kendall test, as modified by Yue and Wang [28] to account for serial correlation, was employed to detect significant annual upward or downward trends of salmon abundance in all the study rivers using the modifiedmk version 1.6 R package [29]. To detect regime shifts in the total salmon catch time series, a split moving-window boundary analysis [30] was used. This test computes the statistical contrast of the means of two halves of a window (width of 20 years here) as it moves along the time series through ANOVA tests (α = 0.01).

2.3. DNA Extraction and Microsatellite Genotyping

The QIAamp DNA Mini Kit (QIAGEN, IZASA, Madrid, Spain) was used to extract genomic DNA from archival scales. The quality and concentration of DNA were determined using spectrophotometry and were verified by 0.8% agarose gel electrophoresis. The following nine microsatellite loci were analysed: Ssa197, Ssa85 [31], SSOSL417, SSOSL85, SSOLS311 [32], SSOSL438 [33], SsaF43 [34], SSspG7, and SSsp2210 [35]. For each locus, a polymerase chain reaction (PCR) was performed in a volume of 25 μL. PCR reactions contained 1X Mg-free PCR buffer, 1.5–2.5 mM MgCl₂, 0.16 mM of each primer, 200 μM of dNTPs, 1–1.25 U of Biotools HotSplit DNA polymerase (Biotools, Madrid, Spain), and 50–100 ng genomic DNA template. Specific annealing temperatures for the loci were 50 °C (SSOSL85, SSOSL311), 53 °C (SSOSL417, SSOSL438), 55 °C (SsaF43, SSspG7), or 58 °C (Ssa197, SSsp2210, Ssa85).

The amplifications were performed in a PCR machine GeneAmp^® PCR System 9700 (Applied Biosystems, Waltham, MA, USA) using the following conditions, with occasional modifications to adapt to specific primers and/or samples: 95 °C 5 min, 30–40 cycles of 95 °C 20–40 s, 50 °C up to 58 °C, 20–40 s, and 72 °C 20–40 s, with a final extension at 72 °C 10 min. An ABI PRISM 3730 sequencer (Applied Biosystems, Waltham, MA, USA) was used to visualise the PCR products, and allele scoring was performed manually with the Peak Scanner^TM Software v1.0 (Applied Biosystems, Waltham, MA, USA). To avoid cross-contamination while working with the archival samples, the entire process of buffer preparation, DNA extraction, and PCR amplification were performed in a laminar flow hood. The working area and equipment were cleaned with UV light followed by swabbing with ethanol before each new step in the protocol, and strict cleaning procedures were respected by all staff. Negative and positive controls were applied in all steps of the extraction and genotyping process.

2.4. Genetic Diversity

MICRO-CHECKER v2.2.3 (University of Hull, Hull, UK) [36] was used to assess the frequency of null alleles and scoring errors due to stuttering or large allelic drop-out. The combined use of two or three methods has been suggested as the best strategy for minimizing the false-positive and false-negative rates [37]. For this reason, two different analysis methods were additionally used to test for the presence of null alleles: CERVUS v3.0.3 (Field Genetics, London, UK) [38] and ML-NullFreq (Montana State University, Bozeman, MT, USA) [39].

GENEPOP v4.1 (Université de Montpellier, Montpellier, France) [40] was used to perform tests for departure from Hardy–Weinberg equilibrium and linkage disequilibrium for each locus and sample. The statistical significance was evaluated using the Bonferroni corrections. GENETIX v4.05.2 (Université de Montpellier, Montpellier, France) [41] was used to estimate Wright’s fixation index (F_IS) for samples’ deviation from Hardy–Weinberg expectations for heterozygote disequilibrium following Weir and Cockerham [42]. Allele frequencies, number of alleles, and observed (H_o) and expected (H_e) heterozygosities were calculated at the population level with GENETIX v4.05.2 [41].

2.5. Connectivity and Population Structure

FSTAT v2.9.3 (Université de Lausanne, Lausanne, Switzerland) [43] was used to compute allelic richness (A_R) and the genetic differentiation (F_ST) between population pairs using sequential Bonferroni-corrected p-values (10,000 permutations). GENEALEX v6.501 (Australian National University, Acton, Australia) [44] was used to visualise linearised F_ST values F_ST/(1 − F_ST) among locations via principal coordinates analysis (PCoA).

The Bayesian clustering method implemented in STRUCTURE v2.3.4 (Stanford University, Stanford, CA, USA) [45] was applied to further explore population structure. Structure analyses were performed for 1 to 10 clusters (ΔK) with 10 replicates for each simulated cluster. Analyses were run using an admixture model with correlated allele frequencies, 1,000,000 MCMC generations, and a burn-in period of 250,000 steps. Optimal ΔK was determined using STRUCTURE HARVESTER v0.6.94 (University of California, Los Angeles, Los Angeles, CA, USA) [46,47]. Replicates were aggregated using CLUMP v1.1.2 (Stanford University, Stanford, CA, USA) [48] and graphically displayed using DISTRUCT v1.1 (Stanford University, Stanford, CA, USA) [49]. Since STRUCTURE can capture major structures in a dataset but overlook finer scale structure [46], a hierarchical approach was followed. STRUCTURE was further applied independently to the identified genetic clusters until either each river represented its own cluster or optimal ΔK for a group of populations approached 1.

Relationships between genetic differentiation and landscape characteristics were examined using two approaches. First, the Mantel test was used to assess the significance of the correlation between linearised F_ST [50] and pairwise geographic distance between locations (measured as the shortest sea distance between river mouths [51]). Mantel tests with 9999 permutations were conducted in the R v3.3.3. package ade4 v1.7-11 [52,53]. Distance-based canonical redundancy analysis (dbRDA) of pairwise differentiation, implemented in the R package vegan v2.5-2 [54], was used as an additional method to study the correlation between F_STs and geographic distance.

Evidence for demographic bottlenecks was examined following two approaches. First, BOTTLENECK v1.2.02 (INRA-URLB, Montpellier, France) [55,56] was used to test for heterozygosity excess, assuming a two-phase mutation model (TPM) with 80% stepwise mutations (SMM) [57] and 10,000 iterations. The statistical significance of heterozygous excess was tested by one-tailed Wilcoxon’s signed-rank test. In addition, we estimated the ratio of the number of alleles to the range of allele size (M-ratio) following Garza and Williamson [58]. The M-ratio was most likely to correctly detect a population size reduction if the bottleneck was more ancient, prolonged, and had a large Θ value (Θ = 4N_eµ) after the initial population decline [59]. The M_P_VAL software (NOAA Fisheries Santa Cruz Laboratory, University of California, Santa Cruz, CA, USA) [58] was used to estimate the M-ratio, which was compared with a critical value of M (Mc) from a theoretical population in mutation-drift equilibrium, implemented in the CRITICAL_M software (NOAA Fisheries Santa Cruz Laboratory, University of California, Santa Cruz, CA, USA) [58], assuming pre-bottleneck effective population size of 50, 100, 500 or 1000 and a mutation rate (µ) of 5 × 10⁻⁴. According to the recommendations of Perry et al. [60], the proportion of one-step mutations (p_g) was set to 0.22, and the mean size of non-one-step mutations (Δ_g) to 3.1.

BAYESASS v3.0.3 (University of California, Davis, CA, USA) [61] was used to estimate recent migration rates (m) among rivers, using 2,000,000 burn-in and 20,000,000 iterations. Delta values for allele frequencies, inbreeding coefficients, and migration rates were set to 0.6, 0.7, and 1, respectively; these values provide adequate mixing within the ideal range of 20% and 60% [62].

The coalescent method implemented in MIGRATE-N v3.2.7 (Evolution and Genomics, Ballwin, MO, USA) [63,64] was chosen to explore migration rates between and within rivers. Estimations of mutation-scaled migration rates M (M = m/µ) and Θ (Θ = 4N_eµ) were calculated using a Bayesian search strategy and a Brownian motion microsatellite model. Parameter space was searched using 10 short chains and 1 long chain with 3 replicates for 20,000,000 generations, an increment step of 20, and a burn-in of 250,000. Likewise, the parameter space was explored using 4 chains with an adaptive heating scheme (temperatures: 1.0, 1.2, 1.5, 3.0) to ensure that run results did not reflect local likelihood peaks. Finally, the Bayesian assignment procedure of Rannala and Mountain [65] implemented in GENECLASS v2 (INRA-URLB, Montpellier, France) [66] was used to identify putative first-generation migrants in Spanish rivers. The Paetkau et al. [67] resampling method was used to estimate the probability of each individual being a migrant based on 10,000 simulated individuals and a type 1 error threshold significance of 0.01.

2.6. Demographic Parameters

Estimates of census population size (N_c) were determined by the harmonic mean of ten years of annual river catches of returning adults collected by local management authorities. The effective number of breeders within a reproductive year (N_b) was calculated using the linkage disequilibrium (LD) method implemented in NeESTIMATOR v2.1 (Department of Agriculture and Fisheries, Queensland Government, Brisbane, Australia) [68]. A minimum allele frequency cutoff value of 0.02 was employed and 95% confidence intervals were obtained using the jack-knife method. This approach is based on LD between alleles from unlinked neutral loci within populations with random mating and implements a bias correction under a wide range of sample sizes.

Several studies [69,70,71,72,73] have shown that N_b (LD) estimates are robust to predict gene flow when migration rates are low and flow occurs between populations with weak genetic differentiation, even with small effective population sizes and missing data adjustment. This method appears to be the most suitable to estimate effective population size, showing consistent values across different demographic scenarios. In cases with overlapping generations, N_b estimates can be biased. For this reason, we applied the method developed by Waples et al. [74], in which two simple life-history traits are used to adjust genetic estimates of N_b to correct biases due to age structure. N_b was adjusted using the ratio between adult life span (AL) and age-at-maturity (α), following the equation:

N_b(adj) = N_b/(1.103 − 0.245Log (AL/α))

N_b is more easily quantifiable and thus can be a useful tool for managers, but remains less used than N_e [9]. Adjusted age-at-spawning was estimated between 2.65 to 3.04 years in Spanish rivers [18]; hence, an average age-at-maturity of 2.8 years was assumed for all rivers. The AL value was calculated as described by Waples et al. [74], using a maximum breeding age (ω) of 5 years for Atlantic salmon. The N_e(adj) was calculated using the equation proposed by Waples et al. [74]:

N_e(adj) = N_b(adj)/(0.485 + 0.758Log (AL/α))

The effective size ratios N_b(adj)/N_c and N_e(adj)/N_c also were calculated following Perrier et al. [9].

3. Results

3.1. Long-Term Trends in Population Abundance

Total salmon abundance in the study area decreased significantly over the 1950–2022 period (Mann–Kendall test; τ = −0.493, adjusted Kendall p < 0.001, Sen’s slope = −13% change per decade). The same downward trend was observed in the individual rivers (Navia, τ = −0.611, p < 0.001, Sen’s slope = −3.8%; Narcea τ = −0.316, p < 0.001, Sen’s slope = −9%; Sella, τ = −0.494, p < 0.001, Sen’s slope = −11%; Deva-Cares, τ = −0.674, p < 0.001, Sen’s slope = −15.1%; Cantabria’s rivers, τ = −0.708, p < 0.001, Sen’s slope = −15%; Bidasoa, τ = −0.210, p < 0.05, Sen’s slope = −6.1%) (Figure 2). There was a regime shift in CPUE in the early 1970s followed by an abrupt decline from 1988–1989 and a more recent shift in 2008–2009 (Figure 2). Total salmon abundance dropped by 45% from the 1950–1972 period (mean = 5519.4 ± 1770.9) to the 1973–1988 period (mean = 3020.6 ± 1207.9), decreasing from this period to 1989–2008 (mean = 1728.9 ± 563.1) by 43%. Finally, salmon populations were further reduced by 52% from 2009 onwards (mean = 820.6 ± 343.7). On the whole, total salmon abundance was reduced by 85% from 1950–1972 to 2009–2022.

3.2. Genetic Diversity

Variations at nine microsatellite loci were examined in archival samples (the 1950s) from eight Spanish rivers. The locus Ssa85 was removed from the data set because it was difficult to amplify in some individuals and because of the occurrence of some alleles in the Navia, Narcea, and Asón rivers using MICRO-CHECKER (at the Bonferroni confidence level), CERVUS (F_(null) ≥ 0.20) and ML-NullFreq (F_(null) ≥ 0.11, p < 0.05) analyses.

A total of 110 alleles were found among all archival salmon samples. The mean number of alleles per locus ranged from 6.13 for the River Pas to 8.75 for the River Bidasoa (

Table 1. Genetic diversity indices in archival samples (the 1950s) of Atlantic salmon from Spain: sample size (N), observed number of alleles (A), allelic richness (A_R), expected (H_e) and observed heterozygosity (H_o), F_IS values, and deviations from Hardy-Weinberg equilibrium (HWE). No deviation from HWE expectation was statistically significant after Bonferroni corrections.

River		Total	Locus
			Ssa197	SSOSL85	SSOSL311	SSOSL417	Sssp2210	SsspG7	SsaF43	SSOSL438
Navia	N	48	48	48	48	48	48	48	48	47
	A	7.875	8	9	11	8	8	9	6	4
	AR	7.191	7.706	8.178	9.756	7.098	7.458	8.454	5.387	3.489
	Ho	0.623	0.708	0.625	0.646	0.563	0.792	0.813	0.479	0.362
	He	0.643	0.808	0.659	0.717	0.635	0.806	0.786	0.412	0.317
	FIS	0.040	0.133	0.062	0.11	0.125	0.029	−0.023	−0.154	−0.13
	HWE	0.081	0.047	0.268	0.132	0.112	0.391	0.448	0.115	0.209
Narcea	N	48	48	48	48	48	48	48	48	48
	A	7.875	10	7	9	9	9	10	5	4
	AR	7.515	9.433	6.454	8.918	8.701	8.837	9.116	4.729	3.929
	Ho	0.677	0.729	0.708	0.708	0.688	0.792	0.813	0.604	0.375
	He	0.687	0.8	0.661	0.766	0.707	0.788	0.799	0.632	0.34
	FIS	0.024	0.099	−0.061	0.086	0.038	0.005	−0.006	0.054	−0.092
	HWE	0.192	0.095	0.292	0.156	0.848	0.515	0.538	0.342	0.278
Sella	N	49	49	49	49	49	49	49	49	49
	A	8.75	12	11	9	10	9	11	5	3
	AR	8.217	11.727	10.304	8.405	9.3	8.121	10.183	4.693	3
	Ho	0.722	0.776	0.776	0.776	0.837	0.878	0.694	0.571	0.469
	He	0.731	0.853	0.814	0.758	0.761	0.794	0.779	0.544	0.545
	FIS	0.023	0.101	0.057	−0.012	−0.09	−0.095	0.119	−0.041	0.149
	HWE	0.194	0.061	0.228	0.541	0.141	0.111	0.048	0.416	0.149
Deva-Cares	N	46	46	46	46	46	46	46	46	46
	A	8.5	10	8	12	12	9	9	4	4
	AR	7.811	9.886	6.988	10.749	10.51	8.865	7.792	3.988	3.706
	Ho	0.745	0.804	0.717	0.935	0.739	0.826	0.674	0.717	0.544
	He	0.697	0.844	0.688	0.822	0.729	0.762	0.642	0.604	0.488
	FIS	−0.057	0.058	−0.032	−0.127	−0.004	−0.073	−0.039	−0.177	−0.102
	HWE	0.031	0.213	0.435	0.029	0.572	0.195	0.412	0.063	0.283
Nansa	N	48	44	47	48	47	46	48	48	47
	A	7	10	11	11	6	6	5	4	3
	AR	6.702	9.748	10.347	10.368	5.984	5.509	4.724	4	2.937
	Ho	0.678	0.841	0.745	0.813	0.809	0.565	0.646	0.729	0.277
	He	0.667	0.866	0.77	0.855	0.715	0.599	0.584	0.673	0.274
	FIS	−0.006	0.041	0.043	0.06	−0.12	0.067	−0.096	−0.074	0.002
	HWE	0.472	0.297	0.31	0.218	0.105	0.288	0.238	0.272	0.638
Pas	N	47	46	47	47	45	47	47	47	46
	A	6.125	11	7	8	6	5	5	5	2
	AR	5.96	10.633	6.615	7.741	5.998	4.984	4.969	4.744	1.997
	Ho	0.583	0.783	0.617	0.766	0.733	0.532	0.553	0.596	0.087
	He	0.619	0.839	0.694	0.825	0.745	0.614	0.564	0.585	0.083
	FIS	0.068	0.078	0.121	0.082	0.027	0.144	0.03	−0.007	−0.034
	HWE	0.017	0.14	0.102	0.154	0.412	0.089	0.452	0.555	1
Asón	N	48	48	48	48	48	48	48	48	48
	A	7.625	9	7	10	7	9	9	5	5
	AR	7.191	8.378	6.71	9.116	6.706	8.417	8.688	4.928	4.587
	Ho	0.703	0.792	0.813	0.667	0.854	0.688	0.75	0.521	0.542
	He	0.683	0.774	0.724	0.736	0.708	0.725	0.815	0.469	0.516
	FIS	−0.019	−0.012	−0.112	0.104	−0.196	0.062	0.09	−0.101	−0.04
	HWE	0.289	0.539	0.12	0.131	0.013	0.297	0.131	0.225	0.441
Bidasoa	N	40	38	37	40	35	40	40	40	36
	A	8.75	8	9	12	10	8	11	7	5
	AR	8.67	7.989	8.89	11.845	10	7.861	10.971	6.861	4.944
	Ho	0.716	0.868	0.73	0.825	0.714	0.7	0.9	0.575	0.417
	He	0.728	0.79	0.721	0.873	0.709	0.787	0.875	0.63	0.437
	FIS	0.029	−0.086	0.002	0.067	0.008	0.123	−0.016	0.099	0.061
	HWE	0.138	0.194	0.574	0.173	0.556	0.087	0.541	0.172	0.391

). The number of alleles varied from two at locus SSOSL438 to 12 at locus SSOSL311. Allelic richness showed significant differences among rivers (ANOVA, p < 0.05). Higher values were observed in the rivers Sella, Deva-Cares, and Bidasoa (A_R mean = 8.2, range: 7.8–8.7), i.e., those with higher flows. Lower values were found in the rivers Navia, Narcea, Nansa, Pas, and Asón (A_R mean = 7.1, range: 6.0–7.5), i.e., those with lower flows. Observed heterozygosity ranged from 0.58 (the River Pas) to 0.75 (the River Deva-Cares). No statistically significant deviation from Hardy–Weinberg equilibrium or any evidence of linkage between pairs of loci was detected.

Results from the comparison of the genetic diversity parameters (A_R, H_e, and H_o) and effective population size N_e between the archival (the 1950s) and contemporary (the 1990s and 2000s) samples of Spanish Atlantic salmon are shown in

Table 2. Comparison of average allelic richness (A_R), and expected heterozygosity (H_e), observed heterozygosity (H_o), and effective population size (N_e) between archival (the 1950s) and contemporary (the 1990s and 2000s) samples of Atlantic salmon from Spain.

River	Period	Av. AR	Av. He	Av. Ho	Av. Ne	Sources
Navia	1950s	7.19	0.643	0.623	667	Present study
	1990s	-	-	-	-
Narcea	1950s	6.68	0.617	0.607	325	Present study
	1990s	6.32	0.536	0.487	114	[75]
Sella	1950s	8.22	0.731	0.722	1088	Present study
	1990s	6.65	0.532	0.487	107	[75]
Deva-Cares	1950s	7.81	0.697	0.745	729	Present study
	1990s	5.97	0.689	0.564	-	[18,19,75,76]
	2000s	4.68	0.660	0.690	-	[19]
Nansa	1950s	6.70	0.667	0.678	197	Present study
	1990s	5.34	0.772	0.620	68	[18,19]
	2000s	4.88	0.680	0.710	-	[19]
Pas	1950s	5.96	0.619	0.583	180	Present study
	1990s	5.15	0.769	0.650	-	[18,19]
	2000s	4.48	0.670	0.670	-	[19]
Asón	1950s	7.19	0.683	0.703	1051	Present study
	1990s	5.20	0.773	0.650	42	[18,19]
	2000s	4.99	0.700	0.740	-	[19]
Bidasoa	1950s	8.67	0.728	0.716	689	Present study
	2000s	-	0.850	0.750	-	[77]

. Average allelic richness A_R for the historical (the 1950s) samples was higher than contemporary averages. However, the average expected He and observed Ho heterozygosity in archival samples were not apparently different from those of contemporary populations, with the exception of the lower 1990s values for the rivers Narcea and Sella. Finally, the average effective population size N_e for the archival 1950s samples was always higher than contemporary values, showing an evident decline.

All study rivers contained private alleles, ranging in number from 1 to 8 (for the River Bidasoa) (

Table 3. Private alleles in archival samples (the 1950s) of Atlantic salmon from Spain after comparison with the SALSEA baseline. For each microsatellite locus, the first column indicates the allele and the second its frequency (in percentage) in the river. The number of private alleles in each river and locus (#) are also included.

River	# Private Alleles per River	Ssa197		SSOSL85		SSOSL311		SSOSL417		Sssp2210		SSspG7		SsaF43		SSOSL438
Navia	3	-	-	*214	1.0%	*173	1.0%	*209	1.0%	-	-	-	-	-	-	-	-
Narcea	1	-	-	-	-	-	-	*185	1.0%	-	-	-	-	-	-	-	-
Sella	5	*151	2.0%	-	-	-	-	*159	1.0%	-	-	*106	2.0%	-	-	-	-
												*110	7.1%
												*138	1.0%
Deva-Cares	6	-	-	-	-	*123	1.1%	*189	2.2%	*158	3.3%	-	-	-	-	*128	1.1%
						*171	1.1%	*199	1.1%			-	-	-	-
Nansa	1	-	-	*216	3.2%	-	-	-	-	-	-	-	-	-	-	-	-
Pas	1	-	-	-	-	-	-	-	-	-	-	-	-	*128	6.4%	-	-
Asón	3	-	-	*204	1.0%	*127	1.0%	-	-	*138	3.1%	-	-	-	-	-	-
Bidasoa	8	-	-	-	-	*137	2.5%	*175	2.9%	*170	1.3%	-	-	*112	1.3%	*142	1.4%
						*149	6.3%	*183	5.7%
						*163	3.8%
	# Private alleles per locus	1		3		7		7		3		3		2		2
	Total alleles	14		15		19		18		13		15		9		7

). All the microsatellite loci showed private alleles, with SSOSL311 and SSOSL417 showing the highest number of allelic variants across rivers (7 in both cases). Frequencies for these alleles ranged from 1.04% to 7.10%, the highest value corresponding to SSspG7*110 from the River Sella.

After comparison with the SALSEA baseline for the coincident microsatellite loci (Ssa197*, SSsp2210*, SspG7*, and SsaF43*), a few of these private alleles were found exclusively in Spanish rivers and were absent or extremely rare in northern European samples. For example, Ssa197*151 and SSspG7*110 were only found in the historical samples from the River Sella and the 2000s samples from the rivers Ulla and Eo; SSspG7*106 and SsaF43*128 were found in historical samples from the rivers Sella and Pas, and only (at a very low frequency) in three isolated populations from Ireland and Norway and Iceland, respectively.

3.3. Connectivity and Population Structure

All pairwise F_ST values were statistically significant after the Bonferroni corrections (

Table 4. Pairwise F_ST values in archival samples (the 1950s) of Atlantic salmon from Spain. All values were significant after the Bonferroni correction.

	Navia	Narcea	Sella	Deva-Cares	Nansa	Pas	Asón
Narcea	0.019
Sella	0.054	0.039
Deva-Cares	0.059	0.040	0.008
Nansa	0.064	0.049	0.049	0.035
Pas	0.085	0.079	0.098	0.082	0.083
Asón	0.052	0.048	0.036	0.044	0.052	0.085
Bidasoa	0.033	0.037	0.043	0.052	0.071	0.067	0.052

). They ranged from 0.008 to 0.098 (average F_ST = 0.054). The lowest level of genetic differentiation was observed between the rivers Sella and Deva-Cares (F_ST = 0.008, p < 0.010) and rivers Navia and Narcea (F_ST = 0.019, p < 0.01), whereas the highest level of differentiation was observed between the rivers Sella and Pas (F_ST = 0.098, p < 0.01).

Principal coordinate analysis (PCoA) based on linearised F_STs showed the separation of three groups in the studied rivers (Figure 3). The first axis, which accounted for 54.5% of molecular variance, separated the River Pas population from the rest of the rivers. The second axis, which explained 25.6% of the variance, separated the rivers Narcea, Navia, and Bidasoa from the Sella, Deva-Cares, Nansa, and Asón group, while the River Pas occupied an intermediate position.

The results of the Evanno test suggested that ΔK = 3 was the optimal cluster number for the first hierarchical level for the STRUCTURE analysis. Furthermore, the most significant increase in the LnP(D) values was detected at ΔK = 3. This Bayesian approach confirmed the same three genetic units found in the PCoA analysis (Figure 4). The River Asón population was assigned to the Sella, Deva-Cares, and Nansa clusters, to which it had the greatest percentage of membership (Q value = 62.0%). At the second and third hierarchical levels, all rivers were genetically different, except for Deva-Cares and Sella populations, which were indistinguishable from each other and presented a very low genetic divergence (F_ST = 0.008), so they were considered together for further analyses.

Mantel and dbRDA failed to detect isolation by distance among the analysed populations and significant correlations between genetic and geographic distances (Mantel: r = −0.083, p > 0.05; dbRDA: r = 0.801, p > 0.05). The results of the bottleneck tests did not show any such effect on the studied populations (

Table 5. BOTTLENECK tests of archival samples (the 1950s) of Atlantic salmon from Spain. The expected recent bottleneck is presented as p-values from Wilcoxon’s signed-rank test, assuming a two-phase mutation model (TPM) with 80% stepwise mutations (SMM). M-ratio and critical value of M (Mc) assuming pre-bottleneck effective population size (Ne) of 50, 100, 500, and 1000, and a mutation rate (µ) of 5 × 10⁻⁴.

River	SMM	M-ratio	Mc
			Ne = 50	Ne = 100	Ne = 500	Ne = 1000
Navia	p = 0.981	0.771	0.754	0.744	0.701	0.685
Narcea	p = 0.680	0.781	0.754	0.744	0.701	0.685
Sella	p = 0.422	0.777	0.758	0.746	0.701	0.683
Deva-Cares	p = 0.875	0.782	0.756	0.748	0.703	0.680
Nansa	p = 0.273	0.761	0.754	0.744	0.701	0.685
Pas	p = 0.273	0.802	0.758	0.749	0.703	0.799
Asón	p = 0.809	0.776	0.754	0.744	0.701	0.685
Bidasoa	p = 0.727	0.765	0.756	0.749	0.700	0.677

). The BOTTLENECK test did not detect significant excess heterozygosity for the TPM model with 80% stepwise mutations (SMM) (p > 0.05, Table 5). Similarly, M-ratio values were not significantly lower than the simulated critical value of M (Mc) under mutation-drift equilibrium for the lower and upper values of pre-bottleneck effective population size values (p > 0.05, Table 5). Both tests did not show old or recent bottlenecks in the archival samples of Spanish salmon populations.

Recent migration rates between most isolated populations were low and asymmetric (m = 0.020–0.100) or very low (m < 0.020) (Figure 5). Higher migration rates were detected between the Navia (source) and Narcea (destination) populations (m = 0.119, range 0.099–0.141); Narcea and Navia (m = 0.245, range 0.181–0.309), and the Nansa and Sella + Deva-Cares populations (m = 0.169, range 0.074–0.264). Migration rates estimated by MIGRATE were very low (0.002–0.022) and showed no clear patterns. In addition, several of the chains did not reach convergence.

Finally, GENECLASS analysis identified nine individuals as migrants (p < 0.01) (

Table 6. Estimates of the census size (N_c), the effective number of breeders (N_b, [CI 95%]), the adjusted number of breeders (N_b(adj)), the adjusted effective size (N_e(adj)), and the effective size ratios N_b(adj)/N_c, and N_e(adj)/N_c in archival samples (the 1950s) of Atlantic salmon from Spain. Populations in the Deva-Cares and Sella were pooled together after considering the 1st-level STRUCTURE clusters.

River	Nc	Migrants Removed	Nb	Nb(adj)	Ne(adj)	Nb(adj)/Nc	Ne(adj)/Nc
			[CI 95%]	[CI 95%]	[CI 95%]
Navia	701	1	320 [162–∞]	285 [145–∞]	667 [338–∞]	0.407	0.951
Narcea	1298	2	156 [80–∞]	139 [71–∞]	325 [166–∞]	0.107	0.250
Deva-Cares + Sella	2238	2	1167 [161–∞]	1040 [144–∞]	2431 [336–∞]	0.465	1.086
Nansa	124	0	95 [54–275]	85 [48–245]	197 [111–572]	0.681	1.592
Pas	136	1	86 [48–259]	77 [43–230]	180 [100–538]	0.566	1.323
Asón	1051	2	645 [197–∞]	575 [175–∞]	1344 [410–∞]	0.547	1.278
Bidasoa	102	1	331 [83–∞]	295 [74–∞]	689 [173–∞]	2.890	6.753

). The migrants’ populations of origin were in accordance with their percentage of membership estimated from STRUCTURE analysis and their population migration rates estimated by BAYESASS, except for one migrant found in the River Pas. This individual was identified as being from the River Navia but showed a high percentage of membership (Q value = 93.5%) to the Pas cluster, and overall, there was a very low migration rate between rivers (average m = 0.013).

3.4. Demographic Parameters

Estimates of the census size, the effective number of breeders, the adjusted number of breeders, the adjusted effective size, and the effective size ratios are shown in Table 6. N_b(adj) was lowest for the rivers Pas and Nansa (77 and 85 respectively) and highest for the Deva-Cares and Sella populations (1040), with the remaining values being in the 139–575 range. For N_e(adj), a similar trend was observed, with the lowest values found for rivers Pas and Nansa (180, 197) and the highest values for the Deva-Cares and Sella populations (2431), while the rest of the populations varied between 325 and 1344. The River Narcea had the lowest N_b(adj)/N_c value (0.107), while the rest of the populations showed higher values (0.497–0.681). The River Bidasoa presented an unusually high value (2.890). The N_e(adj)/N_c values ranged between 0.250 (Narcea) and 1.592 (Nansa), with the Bidasoa again showing a very high value (6.753).

4. Discussion

4.1. Genetic Diversity

The historical (the 1950s) allelic richness of Spanish salmon populations was comparable to those found in some northern European populations [78,79] and was markedly higher than those observed in the 1990s Spanish samples [18]. These findings indicate a healthy level of genetic diversity prior to a severe loss in the following four decades. This genetic diversity loss contrasts with the increase of allelic richness observed in other salmonid populations after stocking [80] but could be explained by the loss of private alleles, an increase in the fishing effort [76], and the population collapse caused by climate-driven environmental changes in their breeding and feeding areas [10,11], leading to genetic drift in the resulting small populations. The increase in allelic richness in the Spanish 2000s samples [76], surpassing historical pre-stocking values, could be explained as a result of the genetic introgression caused by the stock transfer from northern Europe [12].

The absence of significant deviations from Hardy–Weinberg equilibrium (HWE) is consistent with other wild populations presenting low levels of disturbance (i.e., Icelandic populations [81]). In contrast, Consuegra et al. [18] found significant deviations from HWE (heterozygote deficit) in the rivers Nansa (samples from 1996 and 1997) and Deva (samples from 1993 and 1996), while Griffiths et al. [82] found that genotype frequencies in samples from the rivers Narcea, Asón, and Sella from 2004 significantly deviated from HWE. According to Griffiths et al. [82], these deviations from HWE expectations may be attributed to the decline in salmon population sizes in the Spanish river basins, combined with the effects of historical stocking and supportive breeding measures.

The number of private alleles found in historical samples and their absence in current populations offers further evidence of the loss of local genetic diversity. The existence of this native genetic diversity provides evidence of the relative historical isolation of salmon populations in these rivers before climate change and the effects of stocking on affected populations.

4.2. Connectivity and Population Structure

Historical Spanish populations showed a strong and significant genetic structuring very similar to that found in the River Teno complex (Norway and Finland) salmon population (F_ST = 0.065) [83]. Likewise, a pattern of spatial genetic differentiation was also described in the Spanish rivers, with an overall F_ST value of 0.019 [75]. Separating the pairwise comparisons by year (1993 vs. 1999), an increase in genetic differentiation (from 0.016 to 0.023) was found [75]. Even then, these values are far from those obtained here for the historical (the 1950s) populations (0.054). An overall F_ST value of 0.053 for Spanish rivers sampled in 1988, which dropped to 0.017 in 1996 and recovered to 0.033 in 2007, was reported [21]. These comparisons support the view of original marked genetic differentiation between populations, which was disrupted by human-mediated pressures, such as stock transfer and possible alterations of homing behaviour due to climate change, and a slow later recovery.

The results from the STRUCTURE and PCoA analyses confirmed a significant level of genetic structure in the region, which would correspond to the situation existing before the effects of climate change and foreign stocking intensified in the late 1980s. Indeed, right after the end of stock transfer practices (the 1990s), some Spanish populations (from which the rivers Sella, Narcea, and Cares overlap with our study) could be considered a single genetic unit [21]. In 2007, after several years of supportive breeding, the genetic structure changed into three separate units. This suggests that the removal of human-mediated transfer practices allowed the recovery of the natural genetic structure of salmon populations from Spanish rivers, even if the allocation to genetic units did not exactly match the clusters found in the 1950s.

While most of the allocations to genetic clusters are consistent with geographic location (with neighbouring rivers being the most closely related), other allocations are surprising and warrant more complex explanations. The River Pas population seems to be very isolated genetically, but it is geographically located between the other river basins. On the other hand, the population from the River Bidasoa (the easternmost river basin) was included in the same cluster as the rivers Narcea and Navia (the westernmost river basins). This cluster could be explained by the existence of undocumented pre-1960s human-mediated transfers, which could have homogenised these distant populations. Another compatible scenario is the existence in the past of a widespread lineage, linked by frequent migration, corresponding to the Narcea-Navia-Bidasoa cluster, from which smaller populations (Bidasoa) would have eventually become isolated, diverging into distinct genetic units. Therefore, larger populations (Narcea-Navia) could have retained more of the original signal of panmixia. However, the separation of these populations into different clusters in the lower hierarchical level appears to be more consistent with the geographical location of the populations and should reflect their connectivity and spatial genetic structure more accurately.

A STRUCTURE analysis including samples from the rivers Deva, Asón, Pas, and Nansa, and samples pooled by decades since the 1950s was recently carried out by Ciborowski et al. [19]. As temporal variation added heterogeneity to the dataset and the microsatellites used were not all the same as in the present study, it is difficult to compare the results of both analyses. There were some similarities and some discrepancies. As in the second hierarchical level, each river was assigned to its own cluster. However, populations from the rivers Nansa and Pas could be assigned to the same cluster considering the 1960s data in isolation, while in our study, the Pas population is clearly isolated from the rest at the first hierarchical level.

Significant migration only occurred between contiguous rivers, which were recovered in the same PCoA and STRUCTURE (1st-level) clusters. Estimates obtained by MIGRATE were generally very low and differed from the BAYESASS migration rates in magnitude and direction. Estimates from MIGRATE tend to be comparatively low and do not necessarily reflect connectivity inferred by BAYESASS (commonly interpreted as contemporary connectivity) [24]. In this light, we decided to focus on BAYESASS migration rates for interpretation.

Estimated migration rates for Atlantic salmon have been rarely reported for the Spanish rivers. An earlier study in the 1990s [18], using microtagging-recapture and Bayesian assignment tests based on microsatellites to estimate the proportion of salmon migrants among the rivers Asón, Pas, and Nansa, obtained a range from 0.07 to 0.31. Migration rates obtained for most population pairs in our study are within or below this range (even though some of them are far from each other), indicating infrequent departures from homing behaviour among historical populations.

More studies on Atlantic salmon migration rates have been carried out in the rest of its geographical distribution. Several tagging studies have shown that straying in salmonids is rare, typically ranging between 1 and 10% in geographically close rivers [24]. In the River Teno (Norway and Finland, samples from 1979 to 2001), recent migration rates (using BAYESASS) varied between 0.003–0.260 [84]. In Canadian rivers from Newfoundland and Labrador, a similar range was obtained (0.011–0.303) using contemporary samples and the same methodology [24]. Both studies considered such migration rates as low. Once again, estimates from our historical populations seem to overlap with those previous rates, but with lower maximum values. It is remarkable that migration from larger to smaller populations and vice versa were detected. Some authors [18] suggested that larger, more stable populations should act as a sink attracting strays from neighbouring rivers, while others [85,86,87] found that larger populations act as a source. Empirical studies show that larger populations can behave as sinks and that patterns can switch over evolutionary time scales [24]. In any case, we found that gene flow from smaller to larger populations was quantitatively prevailing. On the other hand, warmer and drier local environmental conditions seem to cause increased straying between neighbouring rivers in southern populations of Atlantic salmon, resulting in a higher gene flow between close populations [77]. Thus, ongoing climate change might increase straying rates in the future. Values obtained in current populations from Spain could be compared with migration rates from historical populations to extrapolate and test this hypothesis across the whole Cantabrian watershed.

The absence of isolation by distance is not uncommon in Atlantic salmon populations at large spatial scales [24,88]. Considering the strong homing behaviour of the species, it seems intuitive that genetic distances between populations from different river basins have an important geographical component. However, it seems more likely that the spatial genetic structure in the Spanish populations may be determined by the geological and ecological characteristics of the rivers (isolation by the environment [89]). Thus, the studied rivers can be classified into Atlantic siliceous rivers (Navia, Narcea), Atlantic calcareous rivers (Sella, Deva-Cares, Nansa, Pas, Asón), and Pyrenean rivers (Bidasoa). Salmon are more likely to stray into rivers that are ecologically similar to their river of origin [90], which is consistent with the correspondence between the 1st-level genetic clusters and the grouping of the rivers on the former classification. Other environmental features of the rivers (temperature, water discharge, etc.) could explain further subdivision in 2nd-level clusters and isolation of the River Pas, the closer relationship between the Sella and Deva-Cares samples, or the exclusion of Nansa (geographically closer to Sella and Deva-Cares).

4.3. Demographic Parameters

The population abundance trends estimated in this study are consistent with previous studies in the area [11], showing a widespread and strong population decline linked to climate change from the 1980s. Consequently, the estimates of the census size, the effective and adjusted number of breeders, the adjusted effective size, and the effective size ratios in historical samples were generally closer to values found in well-preserved, high-latitude salmon populations than to current Spanish populations, which experienced a strong decline linked to climate change from the 1980s [10,11]. For example, the average N_b(adj) from our study was 356, the average N_b(adj) estimate from Canadian (Quebec) populations was 211 [9], and the average N_b estimated for some Spanish rivers was 107.5 [75]. The same applies to N_e(adj): our estimates from historical populations ranged from 180 to 2431 (average 833), while estimates for the 1990s and 2000s populations ranged from 37 to 96 (average 60 [76]), 64–260 [17,75] or 38–175 [18]. These results are consistent with the general declines in Ne observed in salmon populations worldwide [91].

As evidence of the reasonably good conservation status of the historical populations, all had N_e above the minimum threshold (N_e = 95) for retaining 90% of genetic diversity over 100 years [76], and several showed N_e values close to or well over the threshold (N_e = 1000), below which long-term maintenance of evolutionary potential is uncertain [92]. Considering other N_e estimates from the study area [76], which included three of the studied rivers (Cares, Sella, and Narcea), it is established that the decline in effective population size between the 1990s and 2007–2008 already started between the 1950s and the 1990s, probably due to marine overfishing and climate-driven environmental changes in both freshwater and marine habitats. Contrary to Consuegra et al. [18], who found evidence of bottlenecks in the more recent (1993–1998) samples of the populations from the rivers Deva, Nansa, Pas, and Asón, there was no signal of bottlenecking in the populations sampled in the 1950s. This is other evidence of the good conservation status of the historical populations that highlights the importance of their decline in the last decades of the 20th century.

5. Conclusions

Historical (the 1950s) Atlantic salmon populations from Spain showed a remarkable native allelic uniqueness, different indicators of good conservation status (not dissimilar from current Canadian and Scandinavian populations), and a marked population structure with limited migration between close river basins. Since then, study populations have experienced a continuous decline that was accelerated in the 1980s, which could be attributed to intense climate-driven biophysical changes in their marine and freshwater habitats. This decline contributed to a reduction in genetic diversity since the 1990s, which was reversed due to massive stocking with non-native stocks, a practice that could not prevent a loss of native alleles and the population collapse. Furthermore, the synergistic effects of a warmer climate, drier rivers, and the homogenizing effect of human-mediated introgression have driven the genetic erosion of native population structure.

Despite efforts to recover targeted stocks, the current situation is alarming and merits urgent measures to prevent local extinction. Stocking practices have been highly detrimental and their effects, combined with those of climate change, can be considered the main anthropogenic impact affecting salmon survival and the integrity of the species’ genome regionally. Comparisons of current populations with the benchmark provided in our study could help identify which factors have affected salmon populations the most in this area in the last 70 years, in order to guide conservation programs and evaluate the adaptive potential of the species under ongoing climate change. Efforts should target recovering the historical structure of populations and preserving their unique gene pools.