Phylogeography and Conservation Genetics of the Ibero-Balearic Three-Spined Stickleback (Gasterosteus aculeatus)

RESEARCH ARTICLE Phylogeography and Conservation Genetics of the Ibero-Balearic Three-Spined Stickleback (Gasterosteus aculeatus) Marta Vila 1 *, Miguel Hermida 2 , Carlos Ferna ´ndez 2 , Silvia Perea 3 , Ignacio Doadrio 3 , Rafaela Amaro 2 , Eduardo San Miguel 2 1Universidade da Coruña, Evolutionary Biology Group (GIBE), Facultade de Ciencias, Campus da Zapateira, A Coruña, Spain, 2Universidade de Santiago de Compostela, Departamento de Xene ´tica, Facultade de Veterinaria, Avenida Carballo Calero s/n, Lugo, Spain, 3Museo Nacional de Ciencias Naturales, Departamento de Biodiversidad y Biologı ´a Evolutiva, CSIC, Jose ´Gutie ´rrez Abascal 2, Madrid, Spain *[email protected] Abstract Genetic isolation and drift may imperil peripheral populations of wide-ranging species more than central ones. Therefore, information about species genetic variability and population structure is invaluable for conservation managers. The Iberian populations of three-spined stickleback lie at the southwestern periphery of the European distribution of Gasterosteus aculeatus. This teleost is a protected species in Portugal and Spain and local extinctions have been reported in both countries during the last decades. Our objectives were (i) to determine whether the Iberian populations of G.aculeatus are unique or composed of any of the major evolutionary lineages previously identified and (ii) to assess the evolutionary potential of these peripheral populations. We genotyped 478 individuals from 17 sites at 10 polymorphic microsatellite loci to evaluate the genetic variability and differentiation of the Ibero-Balearic populations. We also sequenced 1,165 bp of the mitochondrial genome in 331 of those individuals in order to complement the estimates of genetic diversity in the Ibero-Balearic region. We predicted the evolutionary potential of the different sites analysed based on the contribution of each of them to total allelic/mitochondrial diversity. An intraspecific phylogeny at European level was reconstructed using our data from the mitochondrial cytochrome b gene (755 bp) and published sequences. The so-called Transatlantic,European and Mediterranean mitochondrial lineages were found to be present in the Ibero-Balearic region. Their phylogeography suggests a history of multiple colonisations. The nuclear results show, however, a strong correlation between population structure and drainage system. The following basins should be prioritised by conservation policies in order to preserve those populations with the highest evolutionary potential: the Portuguese Vouga and Tagus as well as the Spanish Majorca and Limia. Maintenance of their connectivity, control of exotic species and monitoring of habitat properties are strongly recommended in those areas. Genetic variation alone cannot, however, ensure the persistence of these peripheral southern populations of G.aculeatus. On the one hand, the analysis of a historical sample from Eastern Spain (Penyscola) revealed no genetic erosion, which suggests a fairly sudden extinction of that population. On the other hand, the reintroduction program implemented in PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 1 / 29 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Vila M, Hermida M, Ferna ´ndez C, Perea S, Doadrio I, Amaro R, et al. (2017) Phylogeography and Conservation Genetics of the Ibero-Balearic Three-Spined Stickleback (Gasterosteus aculeatus). PLoS ONE 12(1): e0170685. doi:10.1371/journal.pone.0170685 Editor: Tzen-Yuh Chiang, National Cheng Kung University, TAIWAN Received: August 17, 2016 Accepted: January 9, 2017 Published: January 24, 2017 Copyright: ©2017 Vila et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: The mitochondrial haplotypes supporting the results of this article have been deposited on GenBank under accession numbers KX910738–KX910784, whereas the microsatellite datasets are available in the FigShare repository (entry DOI: 10.6084/m9.figshare. 4013079 and DOI: 10.6084/m9.figshare.4018740). A fasta file containing the concatenated 48 mitochondrial haplotypes can be found at FigShare repository (entry DOI: 10.6084/m9.figshare. 4028796). the Valencian Community has mostly failed despite our finding of similar level of genetic diversity between the wild source and the captive-bred released individuals. Introduction The combined effects of isolation and drift may cause peripheral populations of wide-ranging species to be more imperilled than central ones [1,2]. This is because low genetic variation is expected to decrease their potential for continuous adaptation [3]. Therefore, knowledge about their genetic diversity and population structure is invaluable for conservation policy and management. However, policy-makers and managers usually need a more practical approach to apply genetic data to conservation biology [4]. One such example is predicting and ranking the evolutionary potential of different breeds [5] or wild populations [6] by assessing their contribution to global diversity [7]. The decline experienced by peripheral populations is well exemplified by the three-spined stickleback (Gasterosteus aculeatus Linnaeus, 1758), a mainly circumboreal-north-temperate fish, widely known as model species in evolutionary biology [8]. The seminal paper by Foster et al. [9] reviewed the situation of this species at the southern edges of its distribution, highlighting the case of the Spanish populations. Indeed, this teleost has suffered a northward trend of local extinctions in the Iberian Peninsula [10–13]. G.aculeatus is classified as Endangered in Portugal [13], but it is currently absent from the Spanish Catalogue of Endangered Species [14]. However, the three-spined stickleback is included in all sub-national red lists of Spanish regions where the species occurs. Such a regional classification triggered the undertaking of conservation work. For instance, a captive breeding and reintroduction program started in 2002 in the Valencian Community (Eastern Spain) [15]. At that time, conservation decisions based on genetic evidence were hampered by the lack of information about most of the Iberian populations of this species. To date, only Araguas et al. [16] and Sanz et al. [17] have studied the phylogeography and conservation genetics of this species in Catalonia (Northeastern Spain). The former defined up to four conservation units based on the population structuring revealed by nuclear microsatellite markers, whereas the latter grouped those clusters into two Evolutionarily Significant Units (ESUs). In addition, the insular population of G.aculeatus from Majorca was recently defined as a different ESU [18]. By contrast, the knowledge about the evolutionary history and genetic diversity of Iberian G.aculeatus from the Atlantic basins is scarce (but see [17–21]). For further information on the concepts of ESU and conservation unit, readers are referred to Funk et al. [22]. Bearing in mind the noticeable divergence of the Portuguese individuals analysed by Sanz et al. [17] and the fact that the Iberian Peninsula is one of the European freshwater fish biodiversity hotspots, both at inter and intraspecific level [23], we addressed the following questions. First, to determine whether the Iberian populations of G.aculeatus were part of any of the major evolutionary lineages identified by prior literature [20,24] or were genetically unique. Second, to investigate the genetic diversity and population structure of this iconic species along the southwestern edge of its European range. These results were used to predict the evolutionary potential of these peripheral populations as well as to discuss their implications for conservation. Materials and Methods Ethics statement Ethics approval of all procedures involving vertebrate animals is legally required under the Spanish legislation (Royal Decree 1201/2005 and Law 32/2007, on the protection of animals Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 2 / 29 Funding: This work was supported by grants from Xunta de Galicia (PGIDIT06RF026101PR, GRC2014/050) and Ministerio de Ciencia y Tecnologı ´a (CGL2010-15231/BOS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. used for experimentation and other scientific purposes), which is a transposition of the European Directive 86/609/EEC. In agreement with article 18 and annexes VII and XI of the aforementioned Royal Decree, all animal procedures performed as part of the experimental work described in this paper have received prior and explicit approval from the competent authorities, defined in article 3e of the Law, and substantiated in the corresponding regulations of the Spanish autonomous communities. Thus, permissions for fieldwork and the concomitant experimental procedures (nonlethal sampling) were issued by the Xunta de Galicia (permit #35/2007), Diputacio ´n Foral de Bizkaia (reference number 9621) and Parc Natural de S’Albufera (permit issued by the director, Mr. Rabassa, on April 18 th , 2007 to S Perea and I Doadrio) in Spain, and, in application of an analogous transposition of the European Directive, by the Instituto da Conservac¸ão da Natureza (Licenc¸a N˚259/2007/CAPT, Licenc¸a N˚260/2007/ CAPT, Licenc¸a N˚261/2007/CAPT) and Direcc¸ão Geral dos Recursos Florestais (Credencial de pesca N˚ 88/2007, Credencial de pesca N˚ 89/2007, Credencial de pesca N˚ 90/2007) in Portugal. Samples provided by other researchers (S7: Ha¨nnu Ma¨kinen but collected by Ignacio Doadrio, S14: Centro de Investigacio ´n Piscı ´cola El Palmar, S16: Francisco Go ´mez Caruana and S17: Ju¨rgen Geist) were collected following procedures reviewed by their institutions and approved by the competent authorities (CSIC’s Ethics Committee, Generalitat Valenciana, Centro de Acuicultura Experimental and Technische Universita¨t Mu¨nchen, respectively). Sampling Altogether, 478 sticklebacks from 17 freshwater sites were sampled for the present study: 389 of them were obtained at 14 locations to cover most of the hydrogeographical areas were the species is naturally present in Spain and Portugal (Fig 1). Sampling at S1-S13 and S15 took place between June and August of 2007 and 2008, but for S7 (Antela) sampled in 2004. Adult fish were collected using fishing nets. Animals were quickly (<5 min) processed by clipping a piece of the caudal fin with sterile scissors and immediately released at the capture site afterwards. Tissue was preserved in 95% ethanol. We additionally analysed 30 individuals (coded as S14, Valencia) from a translocated population as well as 30 specimens from the extinct population of Penyscola (coded as S16), kept at the National Museum of Natural Sciences (Madrid). We also surveyed 29 individuals (coded as S17, Gu¨nz, a tributary to the Danube) from Southeast Germany in order to have an external reference for the microsatellite results. More details about the samples can be found at S1 Appendix. DNA extraction and scoring of molecular markers Genomic DNA was extracted from fin tissue using DNeasy1Blood & Tissue Kit (QIAGEN). For comparison purposes, we aimed at using the same set of markers as Ma¨kinen et al. [19] and Ma¨kinen & Merila¨[20]. We initially genotyped the 478 sampled individuals (up to 30 per locality) using the same 18 nuclear microsatellites used by Ma¨kinen et al. [19], except for locus stn122, (replaced by stn82 in the present work), which failed to amplify for individuals from Guisande (S4) and Rato (S5) during preliminary tests. Microsatellites were amplified as specified in Pe ´rez-Figueroa et al. [21]. PCR products were separated in an ABI Prism 3730xl Analyzer (Applied Biosystems, ABI). Alleles were scored using GENEMAPPER 4.0 (ABI). We sequenced a total of 1,165 bp of mtDNA in 17–22 individuals chosen at random within each sampling locality. We used the primers described by Ma¨kinen & Merila¨[20] to obtain 755 bp from the cytochrome b (cytb) and 409 bp from the control region (cr). PCR amplifications for cytbwere performed in a final volume of 20 μL containing 100ng DNA template, 1 X Buffer, 2.5mM MgCl 2 , 0.2 mM dNTPs, 0.5 μM each primer, 1 U GoTaq Flexi DNA Polymerase Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 3 / 29 (Promega). Thermocycler profile was as follows: initial denaturation (94˚C, 2’); 35 cycles of denaturation (94˚C, 1’), annealing (56˚C, 1’) and extension (72˚C, 2’); final extension (72˚C, 10’). PCR amplifications for the cr were similar to those of cytb. The only difference was the annealing temperature: 55˚C. PCR products were electrophoresed on 2% agarose gels and visualised under UV light after ethidium bromide staining. Products were purified with MultiScreen kit (Millipore) and used as template for direct sequencing on an ABI Prism 3730xl. DNA sequences were inspected and aligned using SEQSCAPE 2.5 (ABI). Alignments were straightforward, as only a mononucleotide indel was found in the cr fragment. Measures of genetic diversity and paring down of molecular markers Microsatellites: Genetic diversity, HW and gametic equilibrium per sampling locality. The number of alleles per locus ranged between six (loci stn38 and stn46) and 54 (locus 1125pbbe) [28]. Basic descriptors of genetic diversity per locus can be found at S1 Table. The initial inspection of the 18 loci dataset revealed markers stn3, stn12 and stn174 to show intermediate alleles (1 bp difference). Therefore, these three loci were excluded from further analyses. The frequency of null alleles was then calculated using Oosterhout’s estimator as implemented in MICRO-CHECKER 2.2.3 [29]. Eight markers showed evidence for null alleles (S2 Table), so we ended up with a final robust dataset of ten loci [30]: 1125pbbe [31], stn19, stn21, stn38, stn46, stn57, stn79, stn110, stn163 and stn195 [32]. Fig 1. Approximate distribution area of G.aculeatus across Portugal and Spain (redrawn from Ribeiro et al. [13] and Doadrio et al. [23]). Readers are referred elsewhere [16,25–27] for a more detailed geographic distribution of the species in Northern Spain. Boundaries of icthyographic provinces are marked by dotted lines. Sampling locations surveyed for the present study (red squares) are coded as in Table 1. Yellow stars correspond to localities where the species has been reintroduced. The asterisk (S14*) indicates that the sample we analysed (Valencia, yellow star) comes from an ex-situ breeding program, using individuals from Catalonia (red circle). The cross indicates that site Penyscola (S16) went extinct in the early 1990’s. doi:10.1371/journal.pone.0170685.g001 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 4 / 29 Number of alleles observed and averaged over loci, observed and unbiased expected heterozygosities, as well as analyses of Hardy-Weinberg (HW) segregations were computed with GENODIVE 2.0b25 [33]. Allelic richness and allelic private richness were obtained using HP-RARE [34]. Tests for gametic phase disequilibrium were computed using the web version of GENEPOP 4.2 [35]. For this, we applied the test for each pair of loci in each population and default parameters as well as Fisher’s method to test the null hypothesis of random association for each locus pair across all populations. Newly obtained mitochondrial haplotypes. Haplotypes defined by the cytb and cr fragments, as well as and their frequencies and other standard indices of genetic variation such as the number of segregating sites (S), nucleotide diversity per gene (π), haplotype diversity (Hd) and the average number of nucleotide differences (k) were calculated in DNAsp 5.1 [36]. New haplotypes were deposited in GenBank (S3 Table). Identification of nuclear gene pools As isolation-by-distance (IBD) scenarios are not suitable for model-based Bayesian clustering of genotypes, we tested the correlation between genetic (F ST /(1F ST )) and geographical distances (raw and log 10 transformed). We calculated the matrix of shortest pairwise water distances among sampling sites using the Ruler Path option in Google Earth. For comparative purposes, linear distances were calculated for each pair of sampling locations with Geographic Distance Matrix Generator 1.2.3 [37] (S4 Table). Mantel tests and their significance after 1000 permutations were run in GENODIVE. After ruling out IBD, we used STRUCTURE 2.3.4 [38] and BAPS 6.0 [39] to unravel the nuclear genetic structure. These two programs implement model-based Bayesian clustering algorithms that minimise Hardy-Weinberg and linkage disequilibrium. We ran STRUCTURE to cluster individuals without prior information on sample origin; simulations were run assuming the admixture ancestry model and correlated allele frequencies. After preliminary runs (data not shown) aiming at evaluating the Markov chain Monte Carlo (MCMC) length needed for the summary statistics to converge, we set up a burn-in of 200,000 iterations followed by 500,000 iterations for parameter estimation. Each simulation was run 20 times, exploring values for K(the total number of clusters to be constructed in a given simulation) ranging from one to 18. We preliminary determined the number of clusters by applying both the methods of Evanno [40] and Pritchard [38] as implemented in CLUMPAK [41]. We used CLUMPP [42] to permute the admixture coefficients for the 20 independent runs resulting for each K-value. Then, we ran DISTRUCT [43] to visualise the output from CLUMPP. Following Meirmans [44], we discussed the clustering results based on their biological relevance. With regard to BAPS, we firstly applied a non-spatial genetic mixture analysis [45] clustering both groups of individuals (sampled localities) and just individuals. Briefly, this MCMC-based algorithm groups samples into variable user-defined numbers Kof clusters. Then, the best partition of data into Kclusters is identified as the one with the highest marginal log-likelihood. We performed ten independent simulations for each value of Kfrom 1 to 18. Lastly, we performed an admixture analysis based on mixture clustering (of the resulting 17 groups of individuals). Differentiation between sampling sites In the light of the results provided by the clustering algorithms, pairwise differentiation among sampling localities was calculated using the nuclear dataset. Firstly, we calculated the values of the IAM-based θunbiased estimator of F ST [46] and their significance after 10,000 permutations. Then, we obtained the harmonic mean of D[47] across loci. We used GENODIVE for these calculations. Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 5 / 29 Mitochondrial phylogeny The phylogeny of the mitochondrial haplotypes resulting from concatenation of cytband cr was inferred calculating a 95% statistical parsimony network using TCS 1.2 [48]. The resulting haplotype network revealed so many alternative mutational connections (loops) that it was impossible to displayed clearly. Hence, the network for each mitochondrial fragment was calculated separately. Both the lack of phylogenetic resolution as from the cytb+crconcatenated dataset (preliminary Bayesian and Maximum Likelihood trees not shown) and the fact that Malhi et al. [49] only published haplotypes of cytbon the Scottish three-spined sticklebacks led us to the use of the cytbfor phylogenetic analyses at European level. For this, we gathered a 755 bp cytbdataset of 172 haplotypes by aligning and collapsing our Iberian cytb sequences and the homologous fragments reported in prior literature [17,20,24,49–54] (S2 Appendix). We again calculated a 95% statistical parsimony network using TCS. The earliest diverging haplotype/s were inferred using the homologous fragment of a Japanese specimen of G.aculeatus (Accession number AB094627). Due to the difficulty of representing in such a network the actual sample size of some haplotypes, readers are referred to S1 Fig and S2 Appendix for further details on the frequency and geographic distribution of those mitochondrial variants. Population prioritisation for conservation The evolutionary potential of different populations can be predicted from the contribution of each of them to the overall species level allelic diversity [55,56]. Such a contribution was calculated for both the nuclear and mitochondrial datasets using METAPOP 2.0.a1 [57]. We ran the Population Analysis implemented in METAPOP applying rarefaction to correct for sample size and using Nto determine the weight given to each subpopulation when calculating averages. This software was also used to rank populations according to their relative contribution (c GDpool ) to produce a single (hypothetical) pool of maximal gene diversity. For this, we set λ= 1 (equal weights to withinand between-population diversity), 1000 individuals and 2000 steps to perform the simulated annealing algorithm. The rationale of this calculation is that maximisation of genetic diversity is equivalent to maximisation of effective population size. In addition, maximising genetic diversity is expected to lead to maximum allelic richness in the long-term [58]. Results Genetic diversity: microsatellites The final dataset of ten loci contained an overall proportion of missing genotypes of 0.46% (ranging from zero at locus stn57 to 0.8% at 1125pbbe). We scored a total number of 251 alleles. Each individual resulted in a different multilocus genotype, except for two individuals from S13. These two specimens differed, however, in their genotypes for excluded loci 7033pbbe,stn174 and stn132. The total number of alleles per locus ranged between six (stn46) and 54 (1125pbbe) (S1 Table). At population level, the lowest allelic richness was obtained at S1 (Txingudi) and S8 (Salas), whereas the highest was found at S9 (Vouga), a Portuguese locality. Focusing on gene diversity, the lowest values were found at S4 (Guisande) and S12 (Sado), and the highest ones at S3 (Gobelas), S10 (Vouga) and S14 (Valencia). Majorca (S15), Vouga and Tagus (S11) showed the highest private allelic richness (Table 1). We observed two significant deviations from gametic equilibrium at locality level. The two loci involved in each case (stn21xstn195 at S16 (Penyscola) and stn19xstn163 at S11 (Tagus), adjusted P-value for 5% nominal level = 0.0014 and 0.0011, respectively) corresponded to different linkage groups according to [32] (S1 Table). Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 6 / 29 No significant deviations from Hardy-Weinberg (HW) equilibrium were detected at the population level, after correcting for multiple tests (adjusted p-value for 5% nominal level = 0.003). Genetic diversity: mtDNA The concatenated (cytb +cr) data matrix contained 331 sequences and defined 48 haplotypes (47 if only nucleotide substitutions were considered, see below) deposited in FigShare [59] (S3 Table). Overall, nucleotide diversity (±SD) was π= 0.0059 ±0.0032, haplotype diversity Hd = 0.953 ±0.004 and average number of nucleotide differences k= 6.775 ±3.201. The most diverse localities were the Portuguese S9 and S10 (Mondego, Vouga), whereas only one haplotype was found in the following localities: the three Basque sites (S1-S3), S8 (Salas) and the German S17 (Gu¨nz) (Table 2). Variation in the cytochrome b sequences was due to 28 segregating sites, which defined 28 haplotypes. Twenty-four of them were new and deposited in Genbank (Accession numbers KX-910738-KX910761). Twenty-one substitutions were parsimony informative and most changes (82.14%) corresponded to third codon positions and implied no amino acid replacement. Three Northwestern Spanish localities (S4 Guisande, S5 Rato and S7 Antela) were the most diverse as from cytbsequences, whereas eight others resulted monomorphic (the Basque S1-S3, S8 Salas, S13 Mira, S14 Valencia, S16 Penyscola and S17 Gu¨nz) (Fig 2). The haplotype we found at the translocated site of Valencia (S14) did not match the one reported at river Orlina [17], the source of individuals according to the environmental Table 1. Indices of nuclear genetic diversity calculated for the 17 localities where G.aculeatus was sampled. Code Population River Basin N Missing data Num AR PAR H O H E G IS S1 Txingudi Bidasoa 17 0 2.9 2.88 0.27 0.388 0.434 0.106 S2 Castaños Nervio ´n 30 0 4.1 3.9 0.01 0.560 0.547 -0.024 S3 Gobelas Nervio ´n 27 0 6.6 5.89 0.28 0.689 0.655 -0.052 S4 Guisande Miño 30 0 6.7 5.33 0.2 0.380 0.396 0.041 S5 Rato Miño 30 1(stn21) 5.1 4.59 0.16 0.543 0.548 0.009 S6 Asma Miño 30 1(stn19, stn110, stn46, stn163, stn79 and stn195) 6.3 5.43 0.21 0.560 0.584 0.040 S7 Antela Limia 30 3(1125pbbe), 1 (stn21) 9.2 7.24 0.44 0.549 0.570 0.036 S8 Salas Limia 17 1(1125pbbe), 1(stn46), 1(stn79) 2.9 2.89 0.02 0.469 0.467 -0.004 S9 Mondego Mondego 28 0 6.5 5.71 0.53 0.579 0.565 -0.024 S10 Vouga Vouga 30 0 10.0 7.98 0.88 0.683 0.643 -0.062 S11 Tagus Tagus 30 1(stn195), 1(stn21, stn110, stn38) 6.5 5.62 0.71 0.594 0.598 0.007 S12 Sado Sado 30 1(stn163, stn79) 4.4 3.77 0.01 0.331 0.356 0.069 S13 Mira Mira 30 0 3.9 3.49 0.01 0.443 0.419 -0.059 S14 Valencia Orlina 30 0 5.2 4.76 0.04 0.637 0.633 -0.006 S15 Majorca 30 0 9.8 7.42 1.83 0.560 0.584 0.041 S16 Penyscola Ju ´car 30 2(stn38) 5.8 5.16 0.15 0.599 0.584 -0.026 S17 Gu¨nz Danube 29 0 4.0 3.55 0.41 0.538 0.526 -0.023 Calculations based on the 10 microsatellite loci unaffected by null alleles and complying with expected allele sizes. Missing data are displayed as number of individuals failing at a given locus (in parentheses). Num = number of alleles observed and averaged over loci, both AR (allelic richness) and PAR (private allele richness) were calculated using a sample size of 16 diploid individuals and averaged over loci, H O = observed heterozygosity, H E = unbiased expected heterozygosity. None of the localities sampled deviated from Hardy-Weinberg expectations (no value of the inbreeding coefficient (G IS ) significantly differed from zero after Bonferroni correction, test based on 10,000 permutations, adjusted nominal level (5%) = 0.003). The three largest values for H E ,AR and PAR are marked in bold, whereas the three lowest ones appear underlined. doi:10.1371/journal.pone.0170685.t001 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 7 / 29 authorities [15]. Rather, that haplotype matched the one obtained from individuals from Ullals lagoon [17] (one of the two stocks maintained in the Delta de l’Ebre Natural Park, [16]), but for a mononucleotide indel in the control region (S3 Table). The variation found in the control region was defined by 23 substitutions and a mononucleotide indel; all of them were parsimony informative. We found 26 haplotypes defined by substitutions. The mononucleotide insertion defined one more haplotype, shared by some Portuguese samples. The 23 new haplotypes were deposited in Genbank (Accession numbers KX910762-KX910784). The most diverse localities as from cr were the northern Portuguese S9 Mondego and S10 Vouga, followed by S5 Rato and S6 Asma, both at Miño river basin. Again, several localities resulted monomorphic for this marker (the Basque S1-S3, S8 Salas and S17 Gu¨nz) (Fig 2). Identification of nuclear gene pools Genetic and geographic distances were not correlated (p>0.05) regardless of using Mantel’s or Spearman’s rstatistics and raw or log 10 transformed geographic distances. Therefore, Bayesian clustering methods were justified. The surveyed localities showed a high degree of nuclear structuring. Actually, there were different levels of organisation present in the genetic structure. Firstly, all but two sampling sites (Valencia and Penyscola) were differentiated according to the partition of data (K= 17) supported by the highest Pr(K=k) in STRUCTURE (Fig 3A). The proportion of individual membership to the 17 th cluster was very low (only >5% at samples from Txingudi, six individuals from Antela, two from Guisande, one from Asma and another one from Vouga) and therefore difficult to notice (sea green colour) at Fig 3A. Secondly, the K= 12 clustering pattern was highlighted as relevant by both Evanno’s method and the plateau of L(K) (Fig 3A,S2 Fig). It is worth noting that STRUCTURE produced four Table 2. Measurements of mitochondrial genetic diversity for G.aculeatus (concatenated dataset: cytb+cr). Code Locality N S h πHd k S1 Txingudi 17 0 1 0 0 0 S2 Castaños 20 0 1 0 0 0 S3 Gobelas 21 0 1 0 0 0 S4 Guisande 21 6 5 0.0023 0.714 2.695 S5 Rato 20 840.0026 0.658 3.005 S6 Asma 19 5 5 0.0015 0.778 1.789 S7 Antela 20 6 80.0014 0.732 1.579 S8 Salas 17 0 1 0 0 0 S9 Mondego 17 10+1 30.0042 0.662 4.853 (5.382) S10 Vouga 20 12+1 8 0.0024 0.816 2.8 (3.068) S11 Tagus 18 7 3 0.001 0.542 1.137 S12 Sado 21 2 3 0.0009 0.643 1.029 S13 Mira 20 1 2 0.0004 0.505 0.505 S14 Valencia 22 3 3 0.0008 0.628 0.978 S15 Majorca 20 5 60.0005 0.516 0.589 S16 Penyscola 19 3 3 0.0013 0.608 1.509 S17 Gu¨nz 19 0 1 0 0 0 N= number of sequenced individuals; S= number of segregating sites (indel after symbol +); h= number of haplotypes; Hd = haplotype diversity; π= nucleotide diversity; k= average number of nucleotide differences (recalculated using indel as fifth character in parentheses). The three largest values of each measurement are shown in bold; a p-value = 0.052 doi:10.1371/journal.pone.0170685.t002 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 8 / 29 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 9 / 29 result of a second wave of Atlantic migrants expanding southwards (Fig 5). Our interpretation about the number of colonisation events should be used with caution, as it is based on the assumption that each of them was due to a group of migrants carrying little ancestral polymorphism. Confirming whether Western Iberian G.aculeatus experienced admixture of lineages due to different migration events or if they are simply carrying part of the polymorphism of a single ancestral population is beyond the scope of this work. For the time being, we advocate that several waves of North Atlantic migrants colonised Western Iberia. This idea is based on the 25 alternations between full glacial (stadial) and relatively mild conditions occurred during the Last Glacial period (119–8 kya) [62]. Such fluctuations likely made northern G.aculeatus to reach the Western Iberian Margin in several occasions. Indeed, dramatic hydrological changes took place along the Western Iberian Margin during the last 61 kya [63,64]. For instance, the Oceanic Polar Front reached latitude 42˚N four times within the 40–10 kya interval [65]. Bearing in mind the high marine dispersal ability of larvae and juveniles of this species (up to 110 Km in the Northeast Atlantic) [66], the arrival of two waves of migrants to River Tagus remains plausible. Fig 6. Geographic distribution of the three cytb haplotypes of G.aculeatus classified as regional as well as their phylogenetic position. Further information about the geographic occurrence of each haplotype can be found in S2 Appendix. Map reprinted from Presentationmaps.com under a CC BY license, with permission from Michael Roscoe, original copyright 2010. doi:10.1371/journal.pone.0170685.g006 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 16 / 29 The star-like phylogeny centred in the disjunct haplotype H19 indicated a population expansion in Northwestern Spain (Fig 5), as eight private haplotypes unambiguously derived from it (Fig 8). H19 was also present in two Irish individuals (S2 Appendix). At present, we cannot determine whether this finding is a hint of the northwards expansion of this haplotype from Northwestern Spain or if it actually dispersed southwards from northern latitudes and subsequently diversified in Northwestern Spain. Both possibilities are plausible as other intermittent poleward flows have been recorded along the Eastern Boundary of the North Atlantic [67]. Genetic data support such a northward dispersal from Western Iberia to the British Isles for some coastal species [68], but not for others [69]. The star-like phylogeny around the disjunct H4 haplotype (present in Portugal, Black Sea, British Isles and USA) indicates transoceanic dispersal. Its derived H29 was private to the American locality of Maine, whereas haplotype H7 was only found in Scotland and the H127-H129 variants were endemic to two geographically close sites in Southwestern Ireland. Transoceanic dispersal has been previously suggested for G.aculeatus [70,71]. The arrival of North American sticklebacks to the Western Iberian Peninsula is plausible if considering that Fig 7. Geographic distribution of the 23 cytb haplotypes of G.aculeatus classified as local as well as their phylogenetic position. Further information about the geographic occurrence of each haplotype can be found in S2 Appendix. Map reprinted from Presentationmaps.com under a CC BY license, with permission from Michael Roscoe, original copyright 2010. doi:10.1371/journal.pone.0170685.g007 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 17 / 29 the North Atlantic Current was displaced south of 42˚N and this allowed icebergs to drift and melt at Iberian latitudes [72]. The question remains whether the Transatlantic dispersal of G. aculeatus was only eastwards or if any of the haplotypes present in Europe, e.g. H4, dispersed into North America at some point. The phylogeography of the nine-spined stickleback Pungitius pungitius does not support such a westward dispersal [73], but further sampling and markers will be needed to test that hypothesis, proved in other species [74]. The three Basque sites shared haplotype H35 with the Western French sites of Borgneuf and Oleron [20], in agreement with the Iberian Poleward Current flow. The populations surrounding the Bay of Biscay most likely derived from the expansion of the widespread haplotypes H1 and H41. Indeed, the low allelic richness but intermediate gene diversity at Castaños may be indicative of a recolonisation [75]. Such a pattern of fairly low allelic richness but high gene diversity was also found in Valencia (S14), where the translocation from Catalonia somehow mimicked a bottleneck due to a long-distance founding event. Fig 8. Phylogenetic position and geographic origin of the 134 private cytb haplotypes of G.aculeatus.Further information about the geographic occurrence of each haplotype can be found in S2 Appendix. doi:10.1371/journal.pone.0170685.g008 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 18 / 29 Fig 9. Proportional contribution to allelic (A) and mitochondrial diversity (B) of each of the 17 sampling sites using data from nuclear microsatellites and mitochondrial haplotypes (cytb+cr), respectively. Positive values show a loss of diversity when the locality is removed from the pool. Triangles indicate the contribution to overall allelic/mitochondrial diversity. The contribution to withinand betweensubpopulation diversity is represented by the white and grey bars, respectively. Vertical axis: proportional contribution expressed as a percentage. doi:10.1371/journal.pone.0170685.g009 Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 19 / 29 Conservation Genetics The threats reported for G.aculeatus in Portugal and Spain are common to distant areas such as North America: introduced species and human-derived changes on water quantity and quality [9,11,76]. Indeed, the abundance of the Iberian three-spined stickleback is negatively correlated to invasive introduced fish (Gambusia holbrooki,Lepomis gibbosus,Micropterus salmoides) and crayfish (Procambarus clarkii) [25,77]. Less is known, however, about the influence of climate change on the extinction of populations of G.aculeatus at low latitudes. The three-spined stickleback is likely displacing northwards because of global warming [78]. Thus, physiochemical constraints such as a rapid increase of water temperature and low dissolved oxygen may be a threat in peripheral isolated populations of southern latitudes [79], particularly in isolated ponds [80]. Prioritising efforts. Majorca (S15), Vouga (S10) and Antela (S7) were unambiguously indicated by both the nuclear allelic and mitochondrial diversity analyses as the Ibero-Balearic populations of G.aculeatus on whose habitats conservation efforts should be focused, as their loss would dramatically impact on the overall allelic and mitochondrial diversity. Prioritisation of highly diverse and effectively large populations (rather than unique but small and depleted ones) is a better approach to maintain the evolutionary potential of a species. This is because allelic-diversity measures are more correlated with long-term and total response to selection, an observation valid to unlinked neutral markers [55,56]. The fact that we did not include data by Araguas et al. [16] in our analyses (see S1 Appendix) prevents us from drawing robust conclusions at national level. However, our results will be most valuable not only for the Portuguese but also for the Spanish authorities. This is because conservation policies are mainly dependent on the Spanish regional governments. We found that G.aculeatus from Majorca contributed the most to the overall nuclear and mitochondrial diversity of this species in the Ibero-Balearic region. This result reinforces the interest in preserving this isolated peripheral population from the loss of genetic diversity due to drift. This insular population was recently postulated as an Evolutionarily Significant Unit (ESU) [18]: not only the three-spined sticklebacks from Majorca were genetically distinct from all mainland populations, but also body size was reduced in adult individuals and resembled sub-adult phenotypes from Northern Europe. Water temperature is likely playing a role in this case, as the breeding period was either extended from spring to November or shifted to autumn [18]. Vouga (S10) and Antela (S7) were the most interesting Iberian sites in terms of prioritising conservation efforts. The loss of both of them would diminish nuclear and mitochondrial diversity in the Ibero-Balearic region about 7 and 20%, respectively. Therefore, fine scale monitoring of G.aculeatus and in these two basins is also needed to evaluate their effective population size. River Vouga was free of dams until 2015 when the Ribeiradio—Ermida hydroelectric project (40˚44’32’’N, 8˚19’10’’W) was fully operative. These two dams are in the middle reach of river Vouga, where G.aculeatus is likely absent from. However, such physical barriers are known to cause hydrological alterations (e.g. flow regime, sedimentation) leading to severe changes in river function (e.g. hypersalinisation of estuaries, temperature fluctuations) [81] that might threat the three-spined stickleback at the lower reach of Vouga. Besides, exotic species such as M.salmoides,L.gibbosus,G.holbrooki and P.clarkii are present in the Vouga and Mondego river basins [82,83]. Our admixture results revealed that Vouga exchanged individuals with other Portuguese basins. On the one hand, one of the four individuals (out of 478) showing significant nuclear admixture was collected at Vouga (Vou06: 61% Vouga, 24% Mira). On the other hand, Vouga and Mondego shared the widespread H1 haplotype and the local H67 variant. These results suggest some degree of dispersal along the Portuguese coast. This is plausible from both biological [84] and oceanographic perspectives. Besides the Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 20 / 29 aforementioned Iberian Poleward Current (currently flowing northwards during winter), the Portugal Coastal Current transports water southwards during summer [64]. Marine dispersal of G.aculeatus was inferred in California between rivers Ventura and Santa Clara (whose mouths are 10 km apart) either because of drainage conjunction during a Holocene lowering of the sea level or because of contemporary and occasional gene flow via ocean dispersal during winter storms [85]. The mouths of Vouga and Mondego are approximately 60 km apart, but the intermediate coastal area is scattered by lagoons and streams that likely shorten the effective distance for dispersal between these two river basins. Indeed, exchanges between the Vouga and Mondego drainages have been reported for other freshwater fish species [86]. Genetically diverse G.aculeatus from Antela (S7) occur in a habitat that suffered drastic human-mediated changes. Not only were two dams built downstream of Antela (see below), but also the original lake was desiccated, artificial channels were constructed, water quality was impaired by pollution and the exotic L.gibbosus and P.clarkii occur in that area. The permanent water body of the Antela Lake is now reduced to a 26 km x 17 m channel, inhabited by the three-spined stickleback, although seasonal floods occur in the surroundings [87]. In the light of the nuclear and mitochondrial diversity obtained (Table 1), we expect the effective population size (N e ) of Antela to be intermediate (contemporary N e ~ 50, long-term N e ~ 300) between the ones obtained for the upper reach of Miño and two of its tributaries (Rato and Guisande) [21]. A similar conclusion of sufficient size and genetic diversity to maintain evolutionary potential needs to be confirmed for Antela, especially in the light of the density fluctuations reported for this species at the Limia basin in 2009 and 2010 [23]. The three-spined stickleback populations inhabiting the Limia river basin were fairly continuous, rather than localised as it happened in other basins such as Miño (Doadrio personal observation). The lack of gene flow between Antela and Salas since 1949 (building of As Conchas dam, downstream of Antela) likely contributed to the differentiation obtained between these two sites: their differentiation was lower than the Basque Castaños-Gobelas intrabasin comparison, but similar to the one between Rato and Asma, within the Miño river basin (Table 3). In addition, Salas is basically a cul-de-sac, as the Salas reservoir (built in 1971) isolates this population from the rest of the Limia river basin. If the Salas low effective population size is confirmed, the only solution to prevent its extinction may be augmentation through a translocation from the same basin. However, such an action makes, a priori, no sense unless the causes for the population decline of these peripheral sites are determined and solved. We postulate river Tagus to be next to Majorca, Vouga/Mondego and Antela in the priority list for conservation managers. We acknowledge that losing Tagus would not substantially affect allelic and mitochondrial diversity. Indeed, we are aware that decisions on the conservation of populations should be taken based on allelic diversity, rather than heterozygosity, especially if few markers are available [56]. However, the contribution of Tagus to a hypothetical pool with maximal gene diversity (GDpool) was the highest (S4 Appendix). This is important because expected heterozygosity may be correlated with short-term response to selection [55], but see [56]. Sample S11 was actually taken in a tributary of Tagus: River Almansor, a Mediterranean-type river. This kind of watershed experiences a predictable annual cycle of autumnwinter floods and extended summer droughts. These changes result in a series of disconnected ponds during the dry season and affect biological community traits [88,89]. This is illustrated by the earlier breeding season of G.aculeatus reported at river Almansor [77], which resembles the shift reported in Majorca [18]. It is therefore tempting to postulate high temperature as a major threat for G.aculeatus in the peripheral populations of the Ibero-Balearic region. Water temperature has been positively correlated with parasite abundance in freshwater G.aculeatus from Iceland [90] and specific parasites have also been suggested as partially responsible for the decline of the Spanish populations of this species [23]. However, the relative importance of Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 21 / 29 water temperature versus other environmental factors for the conservation of G.aculeatus in these peripheral Mediterranean-type rivers is yet to be determined. Abundance of the threespined stickleback was found to increase at river Torgal, a Mediterranean-type tributary of river Mira, during a sequence of dry years to subsequently decrease when wetter conditions resumed [91]. Our sample S13 was actually taken at river Torgal. In addition, a certain level of coexistence between G.aculeatus and P.clarkii has been reported in intermittent streams from Catalonia [25]. Stable versus ephemeral reaches of rivers Ventura and Santa Clara were associated with the main genetic breaks of G.aculeatus from Southern California [85]. It is therefore reasonable to think that such a particular hydrological regime may have influenced the evolutionary history of G.aculeatus from Central/Southern Portugal. About the extinct and translocated sites. The proportional contribution of the population from Penyscola to the Ibero-Balearic allelic and mitochondrial extinct was fairly low. However, its extinction substantially increased the contribution of the other Mediterranean populations, namely the translocated site of Valencia, to a synthetic pool of maximal genetic diversity (S4 Appendix). We failed to identify any of the three main genetic factors to extinction risk reviewed by Frankham [92] at Penyscola. Its intermediate level of genetic diversity (Tables 1and 2) did not match the reduced heterozygosity and allelic richness expected in declining populations. There was no indication of substantial inbreeding either (Table 1). Lastly, mutation accumulation is likely a negligible factor in this case, both in the light of the large number of generations in isolation required to be effective and the fairly low nuclear private allelic richness (Table 1). Therefore, our preliminary conclusion is that the extinction of Penyscola was too abrupt to leave any genetic footprint. This scenario somehow resembles the case of the high nuclear (and mitochondrial) diversities critically reported for the endangered Iberian fish Anaecypris hispanica and attributed to recent fragmentation and decline of its populations [89,93]. Overall, the loss of Penyscola did not substantially affect to the nuclear and mitochondrial diversity present in our 17 sampled sites. The annual release between 2002 and 2014 of almost 108,000 captive bred individuals descending from 77 translocated individuals mostly failed to repopulate two of the Mediterranean sites where G.aculateus became extinct in the 1980/90’s [15]. The Catalonian River Orlina (source of that translocation and analysed by Araguas et al. [16]) and sample S14 showed similar levels of genetic diversity. Our slightly higher values of allelic richness and expected heterozygosity may result from the use of ten markers, i.e. only stn46 and stn195 are in common with Araguas et al. [16], as stn3,stn174 and stn132 had to be excluded from our analyses because of null alleles. The project of reintroducing G.aculeatus at Algemesı ´(39˚14’28.65’’N 0˚22’55.66’’W, Natural Park of Albufera) has come to an end because of the abundance of exotic species. However, the reintroductions at River Bullent (38˚52’49.31’’N 0˚5’6.12’’W, Pego-Oliva Marsh Nature Reserve) pose a more promising scenario, i.e. a lower impact of exotic species and a better status of other threatened species as Valencia hispanica and Salaria fluviatilis. We recommend a genetic monitoring of the individuals to be released in order to ensure the same level of genetic variability found at the source. Translocations of G.aculeatus proved successful elsewhere [50,85,94]. Future studies will have to determine whether the low recapture at the stocking area within the Pego-Oliva Marsh Nature Reserve is due to dispersal through the surrounding wetland [15], to physiochemical constraints or even the use of captive bred individuals [94]. Conclusions The Transatlantic,European and Mediterranean mitochondrial lineages recognised and defined for G.aculeatus by prior literature were present in the Ibero-Balearic region. By contrast, the Black Sea,Adriatic and Irish lineages were absent. Our results suggested that Atlantic Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 22 / 29 three-spined sticklebacks colonised the Iberian Peninsula several times during the end of the Last Glacial period. G.aculeatus showed a strong nuclear population structure highly concordant with the Iberian hydrological pattern. The Portuguese populations of Vouga and Tagus, as well as the Spanish Antela (Northwest Spain) and Majorca should be prioritised by conservation policies. Their loss would severely erode the genetic diversity of G.aculeatus in the Ibero-Balearic region, so maintenance of their connectivity with nearby populations, control of exotic species, and monitoring of habitat properties are recommended in those areas. We would like to highlight the conservation value of the Central and Southern Portuguese populations, not simply because of their genetic and habitat distinctiveness, but rather because of the indication of sporadic coastal dispersal involving Mira and Sado. The site of Penyscola showed an intermediate level of genetic diversity, so that its extinction was likely too abrupt to leave any genetic footprint. Lastly, the ex-situ conservation program implemented in the Valencian Community mostly failed despite ensuring a similar level of genetic diversity between the source and the captive-bred individuals. Environmental changes caused by human activity have been frequently argued to explain local declines and extinctions reported during recent decades, particularly with regard to the influence of exotic species and modification of drainage patterns caused by land-use changes. However, the persistence of peripheral populations at southern latitudes may be also conditioned by the tolerance of G.aculeatus to physiochemical changes related to global warming. Our allelic-diversity results are expected to contain information regarding the evolutionary potential for adaptation to environmental changes. Nevertheless, a deeper knowledge about the relative contribution of extrinsic and intrinsic factors to fitness is needed in order to properly undertake conservation actions that will preserve the three-spined stickleback at low latitudes. Supporting Information S1 Fig. 95% statistical parsimony network obtained for 172 mitochondrial haplotypes (755 bp) of G.aculeatus.Our cytbdata were aligned and collapsed with the homologous fragments reported in Europe by prior literature [17,20,24,49–54]. The basal haplotype/s were inferred using the homologous fragment of a Japanese specimen of G.aculeatus (Accession number AB094627). Readers are referred to S2 Appendix for further details on the frequency and geographic distribution of those mitochondrial variants. (TIF) S2 Fig. Identification of the most likely number of G.aculeatus populations by the analysis of microsatellite data with STRUCTURE 2.3.4. Up: estimated log probability of data for the different number of inferred clusters (K); bars correspond to standard deviation, after 20 independent runs. Down: rate of change in the log probability of data between successive Kvalues (k). It is worth noting that Kmay be underestimated if there is hierarchical structure. (PDF) S1 Table. Basic descriptors of the 18 polymorphic microsatellite loci selected to assess the genetic diversity and differentiation of G.aculeatus.k, number of alleles; AR, allelic richness overall populations (based on sample size of 16 diploid individuals); H O , observed heterozygosity; H E , expected heterozygosity. Significance of F IS was obtained after 1,000 randomisations of alleles;  p0.001. Standard deviations are showed in parentheses. The three highest values of AR and H E are marked in bold, whereas the lowest appear underlined. (DOCX) S2 Table. Proportion of null alleles in loci initially used to genotype the 17 sites where G. aculeatus was sampled for the present study. Loci stn3, stn12 and stn174 excluded from this Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 23 / 29 calculation due to presence of intermediate (1 bp difference) alleles. (DOCX) S3 Table. Geographic and frequency distribution of the 48 mitochondrial haplotypes (cytochrome b (cytb) and control region (cr)) of G.aculeatus (rows) obtained by the present study. (DOCX) S4 Table. Geographic distances (km) between sampling localities. Above diagonal: straight line distances. Below diagonal: shortest pairwise water distances among sampling sites using Google Earth. Note that samples from S14 (Valencia) were geographically coded as located at river Orlina (42˚23’6.09”N, 3˚1’59.29”E), the source of individuals for the translocated specimens founding S14. (DOCX) S1 Appendix. Details about the samples used in the present study. (PDF) S2 Appendix. Original code, frequency and geographic origin of the 172 mitochondrial haplotypes (cytb) of G.aculeatus used in the present work. (XLSX) S3 Appendix. Output of CLUMPAK [41]. (PDF) S4 Appendix. Proportional contribution to allelic and mitochondrial diversity of each of the 15 extant sites using data from nuclear microsatellites and haplotype sequences (cytb +cr), respectively. (PDF) Acknowledgments We thank Sergio de Juan, Patricia Ornelas and Begoña Valca ´rcel for their help during fieldwork. Armando Caballero, Neus Marı ´-Mena, Alejandro Martı ´nez-Abraı ´n, Andre ´s Pe ´rezFigueroa, Manuel Pimentel and three anonymous reviewers improved earlier versions of the manuscript with their comments. Hannu Ma¨kinen, Alison Bell, Jacquelin DeFaveri and Nu ´ria Sanz generously shared some datasets with us. Dr. Ma¨kinen also provided us with DNA samples from individuals previously sampled by Ignacio Doadrio at Antela (S7). Francisco Go ´mez Caruana, Ju¨rgen Geist and Centro de Investigacio ´n Piscı ´cola El Palmar collected the individuals from Penyscola, Gu¨nz and Valencia, respectively. We also thank the forest guards (Xunta de Galicia) for their assistance during fieldwork. David Romero assisted us with artwork. We appreciate the support provided by the Genetics Department (Faculty of Veterinary, USC), Parc Natural de s’Albufera, Fluviatilis Lta. and IZATE Asociacio ´n Medioambiental. We are deeply indebted to David C Lees for his thorough revision of the final version of this work. Author Contributions Conceptualization: MV ID RA ES. Formal analysis: MV MH CF SP. Funding acquisition: MV ID ES. Investigation: MH CF SP ES. Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 24 / 29 Methodology: MV MH CF SP. Project administration: ES. Resources: ES RA ID. Visualization: MV. Writing – original draft: MV. Writing – review & editing: MV MH SP ID ES. References 1. Leppig G, White J. Conservation of peripheral plant populations in California. Madroño. 2006; 53 (3):264–74. 2. Haak A, Williams J, Neville H, Dauwalter D, Colyer W. Conserving Peripheral Trout Populations: the Values and Risks of Life on the Edge. Fisheries. 2010; 35(11):530–49. 3. Allendorf F, Luikart G, Aitken S. Units of conservation. Conservation and the Genetics of Populations. 2nd ed: Wiley-Blackwell; 2013. p. 316–51. 4. Hoban S, Hauffe H, Pe ´rez-Espona S, Arntzen J, Bertorelle G, Bryja J, et al. Bringing genetic diversity to the forefront of conservation policy and management. Conservation Genetic Resources. 2013; 5:593– 8. 5. Medugorac I, Veit-Kensch CE, Ramljak J, Brka M, Markovic B, Stojanovic S, et al. Conservation priorities of genetic diversity in domesticated metapopulations: a study in taurine cattle breeds. Ecology and Evolution. 2011; 1(3):408–20. doi: 10.1002/ece3.39 PMID: 22393510 6. Coleman RA, Weeks AR, Hoffmann AA. Balancing genetic uniqueness and genetic variation in determining conservation and translocation strategies: a comprehensive case study of threatened dwarf galaxias, Galaxiella pusilla (Mack) (Pisces: Galaxiidae). Molecular Ecology. 2013; 22(7):1820–35. doi: 10.1111/mec.12227 PMID: 23432132 7. Caballero A, Rodrı ´guez-Ramilo ST. A new method for the partition of allelic diversity within and between subpopulations. Conservation Genetics. 2010; 11(6):2219–29. 8. Bell M, Foster S. Introduction to the evolutionary biology of the threespine stickleback. In: Bell M, Foster S, editors. The Evolutionary Biology of the Threespine Stickleback. New York: Oxford University Press; 1994. p. 1–27. 9. Foster S, Baker J, Bell M. The case for conserving threespine stickleback populations: Protecting an adaptive radiation. Fisheries. 2003; 28(5):10–8. 10. Go ´mez Caruana F, Sa ´nchez Artal S, Peiro ´Go ´mez S. Estatus del Espinoso (Gasterosteus aculeatus L., 1758) en la Comunidad Valenciana y bases para el programa de su recuperacio ´n. 1997. Available from: http://www.mediterranea.org/cae/statuses.htm 11. Doadrio I. Atlas y Libro Rojo de los Peces Continentales de España. Madrid: CSIC—Ministerio de Medio Ambiente; 2002. 375 p. 12. Garcı ´a Novo F, Martı ´n Cabrera C. Doñana. Agua y biosfera. La biodiversidad en Doñana. Doñana: Confederacio ´n Hidrogra ´fica del Guadalquivir. Ministerio de Medio Ambiente; 2005. p. 329–50. 13. Ribeiro F, Beldade R, Dix M, Bochechas J. Carta Piscı´cola Nacional. Direcc¸ão Geral dos Recursos Florestais—Fluviatilis, Lda. 2007. Available from: http://www.cartapiscicola.org/ 14. Ministerio de Agricultura, Alimentacio ´n y Medio Ambiente. Situacio ´n actual del Listado de Especies Silvestres en Re ´gimen de Proteccio ´n Especial y Cata ´logo Español de Especies Amenazadas. 2015. Available from: http://www.magrama.gob.es/es/biodiversidad/temas/conservacion-de-especies/ especies-proteccion-especial/ce-proteccion-listado-situacion.aspx 15. Generalitat Valenciana. Informe Te ´cnico 06/2015. Conservacio ´n de peces marismeños en la Comunitat Valenciana. Balance de 25 años de trabajo. Servicio de Vida Silvestre. Direccio ´n General de Medio Natural; 2015. Available from: http://www.habitatge.gva.es/documents/91061501/161573226/IT06_2015 +25+a%C3%B1os+conservacion+peces+ribere%C3%B1os/1a6dcf5a-6833-4fe6-b96d11792493e69e 16. Araguas R, Vidal O, Pla C, Sanz N. High genetic diversity of the endangered Iberian three-spined stickleback (Gasterosteus aculeatus) at the Mediterranean edge of its range. Freshwater Biology. 2012; 57 (1):143–54. 17. Sanz N, Araguas R, Vidal O, Viñas J. Glacial refuges for three-spined stickleback in the Iberian Peninsula: mitochondrial DNA phylogeography. Freshwater Biology. 2015; 60(9):1794–809. Population Genetics of Ibero-Balearic Gasterosteus aculeatus PLOS ONE | DOI:10.1371/journal.pone.0170685 January 24, 2017 25 / 29