Introduction
============

Patterns of population connectivity in marine fish depend not only on their dispersal abilities and historical demography, both of which are influenced by the proximity of putative source populations, but also on a complex interplay of these characteristics with environmental factors, such as habitat fragmentation, oceanic currents and climatic fluctuation ([@bib8]). Patterns of genetic structure are, therefore, best understood through the integration of ecological, spatial, environmental and historical information.

The archipelago of the Azores has long represented an attractive target for phylogeographic studies of littoral fishes: the volcanic origin and continued dynamism of these islands, in combination with their remote location, provide an ideal system to investigate the colonisation history of insular species, and pinpoint the traits favourable to their dispersal, establishment and expansion. The majority of studies report genetic isolation in Azorean fishes, because of limited gene flow with the neighbouring island of Madeira, as is the case of the shanny, *Lipophrys pholis* ([@bib53]) and the Montagu\'s blenny, *Coryphoblennius galerita* ([@bib16]; [@bib24]). Studies on other species in the Azores, for example, the ornate wrasse *Thalassoma pavo* ([@bib13]) and the white seabream *Diplodus sargus* ([@bib17]), report low genetic diversity, most probably attributable to the peripheral position of these populations.

The common two-banded seabream, *Diplodus vulgaris* (Geoffroy Saint-Hilaire, 1817; Family Sparidae), is an interesting candidate for a phylogeographic study involving the Azores, as its occurrence has only been reported within the past 20 years. Before that, the only *Diplodus* species known to inhabit Azorean coastal waters was the white seabream, *D. sargus* ([@bib52]). *D. vulgaris* is a benthopelagic fish occurring in rocky and sandy bottoms between the Bay of Biscay, Senegal and the Canary Islands, including the Mediterranean and the Black Sea. Its bathymetric range extends from the shallow subtidal zone down to the depth of ∼90 m ([@bib58]). The distribution area of *D. vulgaris* shows a large degree of overlap with the congeneric species *D. sargus* (Linnaeus, 1758), in which they coexist in relatively segregated niches ([@bib36]).

A detailed checklist of the ichthyofauna of the Azores ([@bib51]) does not include records of *D. vulgaris*, although there is mention of its presence in the nineteenth century ([@bib30]). Informal interviews with fishermen and local divers suggest that this species was not known in these islands before the late 1990s (S Stefanni, personal communication). *D. vulgaris* does not have a common/vernacular name in the Azores and the nineteenth century records were probably because of misidentification or misnaming ([@bib51]). The authors of the present study were able to empirically follow the demographic increase of *D. vulgaris*, with individuals becoming more abundant in successive years, starting with the first recorded occurrence of early life stages in 1997 ([@bib1]) followed by the progressive establishment of all age classes from juveniles to large adults in the following years.

In the present study we test two different scenarios to explain the presence of *D. vulgaris* in the Azores: (1) a recent colonisation resulting from dispersal events and (2) the existence of an ancient relic population that remained undetected until a recent expansion over the past 20 years. With the first scenario, we hypothesise that the population from the Azores has settled as a result of one or more relatively recent waves of colonisers. The apparent absence of this species in the Azores before such events could be explained by the generally unfavourable directions of surface currents that might have prevented dispersal to the Azores. Under this scenario we would expect (1) the Azores population to display a demographic imprint of a recent genetic bottleneck and (2) modest genetic diversity of mitochondrial and nuclear markers, with haplotypes and microsatellites alleles representing a small subset of the main phylogroup of origin (for example, founder effect). Under the second scenario, the expectations are: (1) the old established Azores population to display stronger genetic divergence as a result of long-term isolation, (2) the accumulation of new mutations leading to higher levels of genetic diversity and (3) larger frequencies of private alleles and haplotypes in the refugial population. On the other hand, low genetic diversity may also be expected if the size of the relict population was very small, and this would lead to even more pronounced divergence and lineage sorting, as a result of genetic drift. To test these hypotheses, we compared samples from 11 locations in the NE Atlantic and Mediterranean using mitochondrial DNA (mtDNA; control region), one coding nuclear gene (first intron of the S7 ribosomal protein) and nine neutral microsatellites.

Materials and methods
=====================

Sampling and laboratory procedures
----------------------------------

A total of 307 individuals of *D. vulgaris* were collected from 11 locations in the Mediterranean Sea and the NE Atlantic (mainland coast and islands) ([Figure 1](#fig1){ref-type="fig"} and [Table 1](#tbl1){ref-type="table"}). Specimens from the Azores were collected by spearfishing at different sites around the island of Fayal during the years 2007 and 2008, fish from the Balearics were purchased directly from local fishermen in Palma de Mallorca in 2007, and all other tissue samples came from finclips donated by either fishermen or scientists, as part of collaborative projects between 2007 (Italy), 2008 (Canaries, Mauritania, central and southern Portugal, Cadiz and Malta) and 2009 (Madeira). All samples were stored in 96% ethanol until DNA isolation. This was conducted using commercial kits and following the manufacturers\' instructions: REDExtract-N-Amp kit (Sigma-Aldrich, St Louis, MO, USA) or Mag-Bind Tissue DNA Kit (Omega Biotek, Norcross, GA, USA). Voucher specimens (ethanol preserved tissues or finclips) have been deposited in the Institute for Marine Research/Department of Oceanography and Fisheries of the University of the Azores, Portugal (IMAR/DOP) and Superior Institute of Applied Psychology - Research Unit, Portugal (ISPA-IU) collections.

The 5\' domain of the mitochondrial control region (CR) was successfully amplified and sequenced from 227 specimens for a portion of 344--345 bp in length (Genbank accession numbers: HQ997776--HQ997840 and JQ655470--JQ655598). The first intron of the S7 ribosomal protein gene (S7) was amplified and sequenced for 137 specimens for a portion of 269--271 bp in length (Genbank accession numbers: HQ997841--HQ997893 and JQ624680--JQ624764). Nine polymorphic microsatellites were also screened in all 307 specimens. The loci used were Dvul4, Dvul33, Dvul84 ([@bib48]), Omel2, Omel27, Omel38, Omel54, Omel58 ([@bib47]) and SL26 ([@bib8]).

PCR amplifications of mtDNA CR and nuclear DNA S7 were performed independently with the following pairs of primers: LPro1 (5′-ACTCTCACCCCTAGCTCCCAAAG-3′) and HDL1 (5′-CCTGAAGTAGGAACCAGATGCCAG-3′) ([@bib38]) for the mtDNA control region and S7RPEX1F (5′-TGG CCTCTTCCTTGGCCGTC-3′) and S7RPEX2R (5′-AACTCGTCTGGCTTTTCGCC-3′) ([@bib9]) for the S7 nuclear intron. PCR amplification reactions were performed in a 20 μl total volume with 10 μl of REDExtract-N-ampl PCR reaction mix (Sigma-Aldrich), 0.8 μl of each primer (10 μ[M]{.smallcaps}), 4.4 μl of Sigma water and 4 μl of template DNA. An initial denaturation at 94 °C for 7 min was followed by 35/30 amplification cycles (denaturation at 94 º C for 30/45 s, annealing at 55 °C for 30/45 s and extension at 72 °C for 1 min) and a final extension at 72 °C for 7 min on a Bio-Rad (Hercules, CA, USA) Mycycler thermal cycler (values CR/S7, respectively). The same set of primers was used for sequencing that was carried out at the BMR Genomics (Padua, Italy) facility.

PCR amplification for the nine polymorphic microsatellites was performed in two multiplex reactions: one included Dvul33, Omel2, Omel27, Omel38, Omel54 and Dvul84, and the second reaction pooled Omel58, Dvul4 and SL26. The volume of reagents in the first multiplex reaction were 0.2 μl of each primer, 5 μl of master mix buffer, 0.8 μl of water and 2 μl of 25 ng μl^−1^ DNA. Volumes of reagents in the second reaction were: 0.2 μl of Omel58 primer, 0.25 μl of each Dvul4 and SL26 primer, 5 μl of master mix buffer, 1.2 μl of water and 2 μl of 25 ng μl^−1^ DNA. Reaction conditions for both multiplexes were as follows: 15 min at 95 ºC; 34 cycles at 94 °C for 45 s, 58 °C for 45 s, 72 °C for 45 s; and 72 °C for 45 min.

Genotyping was performed by allele sizing on an ABI 3130xl Genetic Analyser (Applied Biosystems) using forward primers labelled with NED, PET, FAM and VIC dyes, and an internal size standard labelled with LIZ 600 (Applied Biosystems, Foster City, CA, USA). The software GENEMAPPER version 4x (Applied Biosystems) was used to score alleles.

Data analysis
-------------

Sequences, including *D. sargus* (CR: EF468585; S7: EF467745) as outgroup, were aligned using Clustal X ([@bib56]) and CLC Main Workbench (CLC Bio, Aarhus, Denmark).

Genetic diversity indices for CR and S7 markers were estimated using the software ARLEQUIN 3.5 ([@bib22]). In S7, which has a biparental mode of inheritance, most individuals showed heterozygosity, and therefore we used the ELB algorithm ([@bib21]) to phase alleles. Nucleotide diversity (*π*) and haplotype diversity (*h*) were calculated independently for CR and S7 sequences.

Phylogeographic patterns among haplotypes of CR and alleles of S7 were analysed using the Median Joining Network method ([@bib3]) implemented in Network (version 4.6.1.0; [fluxus-engineering.com](http://www.fluxus-engineering.com/)) that infers the most parsimonious branch connections between sampled haplotypes or alleles.

ARLEQUIN was also used to assess deviations from Hardy--Weinberg expectations using randomised distributions of *F*~IS~ values with sequential Bonferroni correction ([@bib46]), and to explore historical demography inferred by mismatch distributions. *P*-values were adjusted to avoid false positives using an optimised false discovery rate approach implemented in *Q*-value ([@bib55]).

In order to assess the statistical power of our data sets to detect population genetic structure, a retrospective power analysis was performed using POWSIM 4.1 ([@bib45]). This program simulates a large number of populations that have been allowed to drift to user-defined levels of genetic divergence (*F*~ST~). The program then samples these populations, and tests the null hypothesis of genetic homogeneity among populations using the χ^2^ approach and Fisher\'s exact test. The percentage of significant outcomes (*P*\<0.05 and *P*\<0.01) from a large number of such samplings is interpreted as the power of the data to detect the defined level of genetic divergence. The level of genetic divergence is user defined by varying combinations of effective population size (*N*~e~) and generations of drift (*t*). For CR, statistical power to detect different levels of *F*~ST~ by varying *t* with a fixed *N*~e~ was tested using 1000 dememorisations, 100 batches and 1000 interactions per batch. Analyses were performed to estimate the risk for false significance (type I error) and to identify the lowest level of genetic divergence that could be detected among our sampling sites with a power ⩾95% and ⩾99% using our actual study data.

Hierarchical analysis of molecular variance ([@bib23]) was employed as a framework to test for alternative geographic groupings among the samples (specific designs are illustrated in the Results). In the case of the S7 intron, the analyses were also run in ARLEQUIN, after allowing the program to reconstruct the haplotypes present.

The software MicroChecker ([@bib57]) was used to check whether microsatellite allelic frequencies were likely to be affected by the presence of null alleles, and ARLEQUIN was employed to test for loci under selection, using the *F*~ST~ outlier approach, under the assumption of hierarchical population structure ([@bib20]). Subsequently, expected and observed heterozygosity (*H*~e~ and *H*~o~) and allelic richness, standardised on the smallest sample size (*R*~A~), were calculated using FSTAT 2.9.3 ([@bib27]), and Kruskal--Wallis nonparametric analysis of variance test was used to evaluate whether different populations had significantly different diversity levels. The *F*~ST~ genetic fixation ([@bib60] and *D*~est~ genetic differentiation ([@bib32]) statistics were estimated with the diveRsity package 1.9.5 ([@bib33]). The variance of each statistic was assessed through the calculation of 10 000 pairwise bootstrapped 95% confidence limits using a bias corrected method that basically recentres the confidence interval (CI) around the initial parameter estimate. We employed both genetic estimators as they present advantages and drawbacks in quantifying population structure (for a discussion see [@bib7]; [@bib32]; [@bib44]; [@bib37]; [@bib59]).

STRUCTURE 2.1 ([@bib42]) was employed for the microsatellite data set to identify the most likely number and composition of population clusters in the data set irrespective of geographic location. An admixture model was used and 10 independent runs were performed for each *K*-value (1--10) using 500 000 iterations with a burn-in period of 100 000. Support for the optimal cluster structure was sought by examining both the averaged log-likelihood across runs for each *K*-value, and the relevant Δ*K* distribution according to [@bib19].

Nonmetric multidimensional scaling (MDS) was inferred from pairwise *F*~ST~ and *D*~est~ matrices using the SMARCOF package 1.5-0 ([@bib15]) and MDS was employed to visually represent the genetic relationships among population samples, using microsatellite-, mtDNA- and S7-based matrices of differentiation. Correlation tests were used to assess the degrees of association between pairwise microsatellites and mtDNA using *F*~ST~ and *D*~est~ values. The ade4 R package 1.6-2 ([@bib18]) was used to perform correlation tests running comparisons between simulated (10 000 permutations) and observed data. BOTTLENECK 1.2 ([@bib41]) was used to test for a recent expansion after a population reduction. The program tests for an excess of heterozygosity (*H*~e~\>*H*~eq~) based on the assumption that, following a bottleneck, the number of alleles in a population (from which the heterozygosity expected at equilibrium *H*~eq~ is inferred) drops faster than classical Wright/Nei expected heterozygosity, *H*~e~ ([@bib10]). The program was run with 5000 iterations assuming a two-phase mutation model with a proportion of mutations that follow the stepwise mutation model set at 80%. In each population, the significance of *H*~e~\>*H*~eq~ was assessed with a Wilcoxon sign-rank test.

We used approximate Bayesian computation (ABC) as implemented in the software DIYABC ([@bib11]) to compare a more elaborate range of colonisation scenarios and identify the most likely historical reconstruction of the origin of *D. vulgaris* in the Azores, combining mtDNA CR and microsatellite data sets (see Results section for explanation of scenarios and relevant probabilistic support). The 'reference table\' was built with 1 000 000 simulated data sets per scenario under three different simulations, the Beaumont, Cornuet--Miller and Guillemaud models, using a combination of summary statistics for the microsatellite data, following [@bib54], whereas for the mtDNA CR we used the following summary statistics: one-sample statistics for mean number of pairwise differences and variance of pairwise differences; two-sample statistics for number of segregating sites, mean of within-sample pairwise differences, mean of between-sample pairwise differences and *F*~*ST*~ between samples. In the Beaumont model, the summary statistics for the mtDNA CR were enriched taking into account Tajima\'s *D*, number of haplotypes, number of segregating sites for one sample and two samples. For the microsatellites we used the default generalised stepwise mutation model with a mean mutation rate across loci set at 10^−3^--10^−4^, whereas for the mtDNA CR we implemented the Hasegawa--Kishino--Yano nucleotide substitution model with nonuniform evolutionary rates among sites (*G*=0.30 and *I*=1) and set a mean mutational rate at 4--12% per million years as a reliable estimate for sparid CR ([@bib5]). Population groupings were informed by prior analysis of population structure (see Results), and focussed specifically on the testing of Azorean colonisation. We used uniform priors for the three scenarios, and initially gave no constraints to effective population sizes (*N*~i~) and coalescent events (*t*~i~) (see [Supplementary Material online](#sup1){ref-type="supplementary-material"}). Then we applied stricter ranges for both *N*~i~ and *t*~i~ to verify robustness of estimates, Azores (AZ) *N*=1--500, Madeira+Portugal South+Cadiz (MAD+PTS+CAD) *N*=500--50 000, Portugal Central (PTC) *N*=500--50 000 and time of divergence for the Azorean sample (*t*1)=1--100 generations. Reliability of scenarios was visualised through principal component analysis, with summary statistics as variables, and posterior probabilities were compared by both logistic regression and direct estimates convergence using 1% of closest simulated data sets (details of additional comparative fits of principal component analysis plots can be found in the [Supplementary Material online](#sup1){ref-type="supplementary-material"}).

Results
=======

Mitochondrial DNA control region
--------------------------------

A 345-bp segment of the mtDNA control region was obtained from 227 individuals, yielding 116 haplotypes with 118 variable sites, 48 of which were parsimony informative. The percentage of private haplotypes (haplotypes only detected in a single location) was ca. 20% in the Azores sample, whereas in other sampling sites ranged between 50% (southern Portugal) and 78% (Cadiz and Mauritania). Overall haplotype diversity was high (*h*=0.983±0.004), whereas the corresponding nucleotide diversity was low (*π*=0.018±0.010). Haplotype diversity ranged from 0.688 to 1.000, whereas nucleotide diversity ranged from 0.012 to 0.017 ([Table 1](#tbl1){ref-type="table"}). Sequences from the Azores showed the lowest haplotypic diversity (0.668±0.035) and the second lowest nucleotide diversity (0.014±0.004). The 23 mtDNA CR sequences of *D. vulgaris* from the Azores contained only 5 haplotypes: the most common haplotypes (SH11, SH5 and SH14) were represented by 19 Azorean fish, 6 from Madeira, 1 from Cadiz, 3 from southern Portugal and 1 from the Balearics.

Estimates of α (type I) error calculated using power analyses indicated that the risk of false significance was below 5% (3.2% χ^2^) or just slightly above (5.7% Fisher\'s exact test). When re-running the simulations using different values for drift in terms of generations, the power analyses indicated that the amount of genetic divergence detectable among our sampling sites with ⩾95% confidence (94.9% χ^2^, 96.8% Fisher\'s exact test) was *F*~ST~ ⩾0.0129 (*N*~e~=2500, *t*=65) and with ⩾99% confidence (98.9% χ^2^, 99.7% Fisher\'s exact test) was *F*~ST~ ⩾0.0157 (*N*~e~=2500, *t*=79). The overall *F*~ST~ value observed was 0.0524. Because of the high number of CR haplotypes, *F*~ST~ values were always smaller than *D*~est~ values ([Figure 2](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}). When testing for pairwise comparisons involving the Azores, both measures of genetic differentiation were significant. For all other comparisons between sampling localities, the two indices agree only in Madeira vs Palermo, Canaries vs southern Portugal, Canaries vs Cadiz and Canaries vs Madeira. For the remaining comparisons, *D*~est~ values were significant and *F*~ST~ were not ([Figure 2](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}).

A series of hierarchical analyses of molecular variance were performed to identify the geographical groupings that best represented the spatial distribution of sequence diversity. The *F*~SC~/*F*~CT~ ratio indicates the ratio between the level of variance within and between groups for each grouping level ([Table 2](#tbl2){ref-type="table"}). The grouping of localities with the lowest ratio was assumed to be the structured grouping that explained the greatest amount of genetic variation. The two-cluster option, in which only the Azores was separated from the other sampling locations, resulted in the largest between-group variance (*F*~CT~=0.107, *P*\<0.001), with little differentiation among populations within the same geographical region (*F*~SC~=0.017, *P*\<0.05). In the three-cluster option, in which the Mediterranean was separated, little improvement was observed: whereas no difference in the within-group variance component was found (*F*~SC~=0.015, *P*\<0.05), the between-group variance was reduced by half (*F*~CT~=0.052, *P*\<0.001). When Madeira was also separated from the rest of the Atlantic sites (four-cluster option), the *F*~SC~/*F*~CT~ ratio dropped considerably mainly because of the negligible within-group variance (*F*~SC~=0.002, *P*=not significant) and no further reduction of the between-group structure (*F*~CT~=0.056, *P*\<0.001). The mismatch analysis showed a unimodal distribution for all locations except the Azores ([Supplementary Figure S1](#sup1){ref-type="supplementary-material"}). The Tau (τ) values, which provide a rough estimate of when rapid changes in the populations have occurred, were similar in all instances ([Supplementary Figure S1](#sup1){ref-type="supplementary-material"}).

The haplotype network did not show any geographic pattern with regards to the Atlantic vs Mediterranean divide, or for any of the individual sampling sites ([Figure 3](#fig3){ref-type="fig"}). The most common haplotype found in the Azores was shared with Madeira, southern Portugal, Cadiz and the Mediterranean, and was included into a large reticulated clade. An independent lineage was found that included haplotypes from the Azores (one shared with Madeira and southern Portugal, and one unique) along with others from Portugal (centre and south), Canaries and the Mediterranean ([Figure 3](#fig3){ref-type="fig"}).

S7 nuclear ribosomal DNA
------------------------

The S7 nuclear region was represented by 127 unique alleles from a total of 274 sequences (138 individuals) and consisted of 271 bp, including 36 variable sites. Estimates of type I error calculated using power analyses indicated that the risk for false significance was below 5% (3.8% χ^2^, 4.6% Fisher\'s exact test). When rerunning the simulations using different values for drift in terms of generations, the power analyses indicated that the amount of genetic divergence detectable among our sampling sites with ⩾95% confidence (94.0% χ^2^, 96.4% Fisher\'s exact test) was *F*~ST~ ⩾0.0035 (*N*~e~=10,000, *t*=70) and with⩾99% confidence (99.0% χ^2^, 99.4% Fisher\'s exact test) was *F*~ST~⩾0.0041 (*N*~e~=10,000, *t*=82). The observed overall *F*~ST~ value was −0.0146.

*F*~ST~ pairwise values were generally smaller for S7 than for CR, and only four comparisons were significant (Torre Guaceto vs Madeira, Canaries, central and southern Portugal; [Figure 2](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}). *D*~est~, on the other hand, proved to be significant in all comparisons involving the Azores and in large part of the other comparisons ([Figure 2](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}). The haplotype network did not reveal any geographical association among haplotypes and it appeared as a single large clade ([Figure 3](#fig3){ref-type="fig"}). Representatives from the Azores were widely spread and shared with all other localities ([Figure 3](#fig3){ref-type="fig"}).

Microsatellites
---------------

No null alleles or departures from neutrality were detected. Overall, private allele counts ranged between 0 (Torre Guaceto) and 26 (southern Portugal) ([Table 1](#tbl1){ref-type="table"}). In the Atlantic, the number of private alleles was smallest in the Azores, with similar low values also in central Portugal, Cadiz and Madeira ([Table 1](#tbl1){ref-type="table"}). Estimates of observed heterozygosity ranged from 0.784 in the Azores to 0.863 in Palermo, and none of the sampling sites were significantly more or less diverse from one another (Kruskal--Wallis test, *P*=0.792; [Table 1](#tbl1){ref-type="table"}). Similarly, allelic richness varied between 7.2 in Azores and 10.1 in Palermo and South Portugal and yet was not significantly different among locations (Kruskal--Wallis test, *P*=0.387). Tests for deviations from Hardy--Weinberg expectation were significant after Bonferroni correction only for Cadiz and Portugal Central. No significant signals of recent bottlenecks were detected, in particular in the focal Azorean population (one-tailed Wilcoxon sign-rank test for *H*~e~ excess, *P*=0.45). Estimates of type I error calculated using power analyses indicated that the risk of false significance was slightly above 5% (6.9% χ^2^, 7.1% Fisher\'s exact test). When rerunning the simulations using different values for drift in terms of generations, the power analyses indicated that the amount of genetic divergence detectable among our sampling sites with ⩾95% confidence (96.1% χ^2^, 93.4% Fisher\'s exact test) was *F*~ST~ ⩾0.0023 (*N*~e~=10 000, *t*=47) and with ⩾99% confidence (98.8% χ^2^, 98.2% Fisher\'s exact test) was *F*~ST~ ⩾0.0028 (*N*~e~=10 000, *t*=57). The overall *F*~ST~ value across populations was however much larger: 0.0240.

The *L*(*K*) distribution obtained with STRUCTURE showed the largest likelihoods for *K*-values between 1 and 4 ([Supplementary Figure S2A](#sup1){ref-type="supplementary-material"}) and the highest Δ*K* value is for *K*=2 ([Supplementary Figure S2B](#sup1){ref-type="supplementary-material"}). The compositional bar plot for *K*=2 provided a hint for the Azores divergence, but for *K*=3--4 the pattern was more evident ([Supplementary Figure S3](#sup1){ref-type="supplementary-material"}). All pairwise *F*~ST~ and *D*~est~ comparisons involving Azores returned the highest values when compared with all other groups, whereas no significance was detected by both indices for intra-Mediterranean pairwise comparisons ([Figures 2 and 4a](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}).

Combined marker analyses
------------------------

MDS plots based on pairwise *F*~ST~ matrices from microsatellites and mtDNA markers returned a strong spatial pattern compared with S7 ([Figures 4a--f](#fig4){ref-type="fig"}). On the other hand, *D*~est~-based MDS plots displayed a geographical pattern only with microsatellite data ([Figure 4b](#fig4){ref-type="fig"}). The finer resolution provided by microsatellite markers highlighted the Atlantic--Mediterranean spatial divide and the arrangement displayed in the *D*~est~ plot clustered locations into two groups accordingly to their oceanographic basin ([Figure 4b](#fig4){ref-type="fig"}). Contrasting patterns were shown by the MDS plots based on CR data, using both measures of differentiation ([Figures 4c and d](#fig4){ref-type="fig"}). The plot based on pairwise *F*~ST~ values resembled the pattern seen for the microsatellites but with a greater distance for the Azores and less discrimination among all Atlantic-Mediterranean locations ([Figure 4c](#fig4){ref-type="fig"}). On the contrary, the MDS plot computed from the *D*~est~ matrix did not show any indication of explicit spatial clustering, with the Azores positioned in proximity of the Balearics, then Cadiz, Madeira and, further distant, southern Portugal ([Figure 4d](#fig4){ref-type="fig"}). There was a pronounced correlation observed between the *F*~ST~ values of microsatellite and CR (*r*=0.851; *P*=0.034), but no correlation was observed when comparing the *D*~est~ values (r=0.166; *P*=0.204) for the same two markers ([Figure 5](#fig5){ref-type="fig"}). No correlation was detected in matrix comparisons between microsatellites and S7 (*F*~ST~: *r*=0.009, *P*=0.354; *D*~est~*: r*=0.243, *P*=0.132) and CR against S7 (*F*~ST~: *r*=0.166, *P*=0.160; *D*~est~: *r*=0.052, *P*=0.369).

The three most realistic ABC scenarios consisted of comparing the Azores with two other potential source populations that comprised Madeira+Portugal and South+Cadiz, on one hand, and Portugal Centre on the other ([Figure 6a](#fig6){ref-type="fig"}). The best supported scenario consisted of an initial split between Madeira+Portugal South+Cadiz and Portugal Centre, followed by dispersal and colonisation of the Azores from Madeira+Portugal South+Cadiz ([Figure 6b](#fig6){ref-type="fig"} and [Supplementary Figure S4](#sup1){ref-type="supplementary-material"}). The summary statistics for this scenario were in the range of the simulated values with the exception of the Cornuet--Miller model, in which the Garza Williamson\'s *M* observed values for the Madeira+Portugal South+Cadiz sample fell in the 0.1% distribution tail. The simulations using the Beaumont, Cornuet--Miller and Guillemaud models provided similar estimates of long-term effective population sizes (*N*~e~) and times of divergence ([Table 3](#tbl3){ref-type="table"}). For the three most realistic scenarios, the prior intervals of *N*=1--500 for Azores and *t*1=1--100 allowed for an optimised fit of the posterior distributions of the observed data ([Figure 6](#fig6){ref-type="fig"}). For the Azores, the *N*~e~ median value varied from 283 to 312 and for Madeira+Portugal South+Cadiz values ranged between 36 700 and 39 600. The median value of the posterior distribution for times of coalescence between these two populations varied between 41.2 (95% CI: 5.5--83.6) and 52.4 (95% CI: 6.32--89.0) generations ([Figure 6b](#fig6){ref-type="fig"} and [Table 3](#tbl3){ref-type="table"}) equivalent to roughly 80--150 years, assuming generation times of 2--3 years ([@bib28]).

Discussion
==========

This phylogeographic study revealed that the recent appearance of *D. vulgaris* in the Azores is unlikely to have been derived from an ancient population, as, despite evidence for strong differentiation, there is no indication of lineage divergence or a conspicuous presence of private alleles. We propose that the origin of the species in the Azores is compatible with the occurrence of rare dispersal events that led to a recent successful establishment of the fish around these oceanic islands. Before discussing these findings in detail, we address two possible caveats: the small sample size for some locations and the short fragment length of mtDNA CR and nuclear S7 markers, 345 and 271 bp long. Although the number of specimens collected for most of the sites was ^\<^50, they provided a sufficient level of diversity to accurately estimate haplotype frequencies ([@bib29]). POWSIM indicated that a type II error because of low power was unlikely (probability \>0.95) (CR: *F*~ST~=0.0129--0.0157, S7: *F*~ST~=0.0035--0.0041 and microsatellites: *F*~ST~=0.0023--0.0028). These values are far lower than the *F*~ST~ values empirically detected between the Azores and any other sampling sites (CR: *F*~ST~=0.1410--0.2370, S7: *F*~ST~=0.0069--0.0583 and microsatellites: *F*~ST~=0.0477--0.0687). In terms of sequence length, besides the POWSIM indications, it is noteworthy that the mtDNA CR and nuclear S7 have displayed adequate informative variation to tackle phylogeographic problems and those exact fragments have been successfully used in other comparable studies (see, for example, [@bib4]; [@bib16], [@bib17]). We therefore conclude that the molecular toolkit employed here is fit to address the focal questions.

Population structure: Azores vs other locations
-----------------------------------------------

Although evidence provided by the three markers is not always congruent, there is nonetheless support for the genetic distinctiveness of the Azores population: (1) analysis of molecular variance on mtDNA CR indicates a four-cluster structure (\[AZ\] \[MAD\] \[Atlantic\] \[Mediterranean\]) as the grouping that explained the greatest amount of genetic variation ([Table 2](#tbl2){ref-type="table"}); (2) mtDNA, S7 and microsatellite *F*~ST~ and *D*~est~ values are significant in all pairwise comparisons involving the Azores ([Figure 2](#fig2){ref-type="fig"} and [Supplementary Table S1](#sup1){ref-type="supplementary-material"}); and (3) microsatellite *F*~ST~ and *D*~est~ and mtDNA CR *F*~ST~ MDS plots return a consistent pattern of Azores separation ([Figures 2](#fig2){ref-type="fig"} and [4](#fig4){ref-type="fig"}). Similarly, this spatial arrangement is supported by the correlation between CR and microsatellite *F*~ST~ matrices ([Figure 5](#fig5){ref-type="fig"}). Only the mtDNA CR *D*~est~ matrix and the S7 marker fail to differentiate the Azorean fish from the others ([Supplementary Table S1](#sup1){ref-type="supplementary-material"} and [Figures 4d--f](#fig4){ref-type="fig"}). This can possibly be due to the presence of several shared haplotypes between the Azores and the other sampling sites for the mtDNA CR ([Figures 1](#fig1){ref-type="fig"} and [3](#fig3){ref-type="fig"}) and the typically slow-evolving nature of the S7 nuclear region ([@bib34]).

The population structure for *D. vulgaris* might be explained by the mechanism that has led this species to reach the Azores. Pioneer expansion, as described in [@bib31], can have a substantial role in marine connectivity. This type of expansion occurs primarily with rare, long-distance dispersal events (that is, small contingents of colonisers settled in the Azores), and populations diverge through a 'reshuffling effect\' where available polymorphisms may be stochastically sorted in the new habitat through coloniser admixture (see [@bib31] and references therein), and not only through genetic drift.

Recent colonisation or demographic expansion of an old population?
------------------------------------------------------------------

The ABC results point to a contemporary (80--150 years ago) divergence between the Azores and the most likely populations of origin. The lack of long-term signature of independent lineages at the most evolutionarily conservative markers constitutes further support of a recent origin of the Azorean population. Moreover, the trend for lower diversity in the Azores provides some indications of reduced effective size in this population (confirmed by the estimates obtained through ABC; the number of microsatellite private alleles is 4 and private mtDNA CR haplotypes is only 1). On the other hand, microsatellites indicate no significant signals of recent bottlenecks in the Azorean population, although the comparatively low number of both private haplotypes and private alleles in that population further points to a small 'reshuffled\' pioneer unit (sensu, [@bib31]) that has diverged from continental source populations, in isolation for a relatively short time.

Simulations using ABC provide estimates that are consistent with the hypothesis that the Azorean *D. vulgaris* reached the archipelago as a result of rare dispersal events from the neighbouring Atlantic islands and continental localities over the past century ([Figure 6b](#fig6){ref-type="fig"}). The assumption that the other oceanic island of Madeira might have acted as a stepping-stone for most of the warm water fish fauna of the Azores ([@bib50]) is a plausible framework to interpret the demographic dynamic of *D. vulgaris*, and the ABC approach offers support to this hypothesis. The broad confidence intervals (95% CI: 11--267 years ago) around the ABC estimates (80--150 years ago) do not exclude the possibility of an even more recent colonisation event that appears to be supported by the lack of visual/catch records of this species. On the other hand, we cannot discount the possibility that *D. vulgaris* sightings may have been favoured by the recent boom in scuba-diving tourism in the Azores and the intensive visual census monitoring programs supported by the Department of Oceanography and Fisheries of the University of the Azores since the 1990s that has increased the focus on coastal species compared with a traditional emphasis on the deeper water species targeted by commercial fisheries and scientific surveys.

Scenario reconstruction
-----------------------

The exhaustive probabilistic approach offered by ABC supports the hypothesis of rare dispersal events that led to the establishment of the Azorean population, providing estimates of *N*~e~ and times of divergence that are similar in the three models ([Figure 6b](#fig6){ref-type="fig"} and [Table 3](#tbl3){ref-type="table"}). The time elapsed between initial colonisation and the onset of rapid population growth is a process that involves ecological and evolutionary adaptations, and thus the time required to overcome these constrains varies among species ([@bib25]). Prolonged time lags might occur during any phase of the invasion process, subsequently affecting their range of expansion ([@bib49]; [@bib14]). It seems therefore possible that *D. vulgaris* reached the Azores in low numbers perhaps as early as in the nineteenth century. Populations may have been small and/or ephemeral and possibly the environmental conditions were not favourable to the expansion of larger demographic units. Nevertheless, an even more recent colonisation of *D. vulgaris* within the past two decades (as faunal records suggest; [@bib1]), or a sudden population expansion following an earlier dispersal events, could be linked to a recent rising of the sea surface temperature (SST). This species is less adapted to cool waters than some of its congenerics, for example, *D. sargus*, that is reported to extend its distribution to higher latitudes (for a review see [@bib39]). The increased SST may have created more suitable habitat conditions for this species in the Azores and favoured its establishment, as temperature-driven distributional shifts in marine organisms have not been uncommon occurrences in the past two decades ([@bib2]; [@bib40]; [@bib6]). The global mean SST values for 1980--1994 and 1995--2005 were 11.1±0.1 °C and 11.7±0.1 °C, respectively ([@bib43]). The SST over the latter period averaged ∼0.6 °C warmer than the preceding 15 years (consistent with Figure 14d in [@bib43]).

The colonisation of an oceanic archipelago by a littoral species in historical times represents a notable example of the effectiveness of sporadic larval dispersal events in suitable climatic conditions ([@bib35]). The temperature rise that affected the North-East Atlantic sea surface waters caused an impact on the circulation patterns around the Azores. In this stretch of ocean, currents are predominantly south-eastern, therefore heading in opposite directions to favour any passive dispersal of larvae from the neighbouring Madeira toward the Azores. However, historical data series indicate episodic anomalies of the SST (mainly in autumn and winter) during which surface currents moved in the opposite direction, and could have enhanced the dispersals of this species to the Azores ([@bib50]). Episodic anomalies are likely to determine temporal 'pulses\' of drops and surges of species diversities with different adaptations to temperature (see, for example, [@bib53]). Moreover, the reproductive period for *D. vulgaris* is between November and March with a pelagic larval duration of 29--58 days ([@bib26]), endowing this species with at least the potential for 'oceanic leaps\' if conditions are suitable.

Implications for fundamental and applied tasks
----------------------------------------------

Collective evidence indicates that the Azorean population of *D. vulgaris* may have originated from waves of colonisers from neighbouring Atlantic localities within the past centuries, and demographically expanded more recently. What is unclear is the number of generations incurred between the initial period of colonisation to the time at which the presence of this species reached a substantial number of individuals to become conspicuous and therefore detectable by local fishermen and divers. Ecological, historical and environmental data support the view that stochastic events, coupled with favourable climatic conditions, may lead to the establishment of new founding populations at the edge of the distribution range. However, the discrepancy between the centuries vs decadal time frames suggested by molecular data and observational records respectively indicate that there may be a time lag between founding at the range edge and a successful population establishment (see, for example, [@bib14]; [@bib25]). This will likely depend on the interplay between the life history traits of the coloniser and the ecosystem structure of the colonised habitat ([@bib49]). Although *D. vulgaris* is not considered a valuable fishery resource, the same dynamics that led to its colonisation (over centuries) and increase (over decades) may well apply to other species with similar life histories, but greater attractiveness for exploitation, such as the case of recently expanding populations of gilt-head seabream (*Sparus aurata*) in northern Europe ([@bib12]). Such demographic processes should be taken into consideration in models that aim at improving both exploitation management and designing protected areas.

Data archiving
==============

Sequence data have been submitted to GenBank: accession numbers HQ997776--HQ997893, JQ655470--JQ655598 and JQ624680--JQ624764.

Genotype data have been submitted to Dryad: [10.5061/dryad.bd6h3](10.5061/dryad.bd6h3).

Data available from the Dryad Digital Repository: <http://dx.doi.org/10.5061/dryad.bd6h3>.

We appreciate the skilful technical assistance provided by S Chenu, A Pereira and Carlotta Sacchi. This work was supported by the MarinERA project 'Marine phylogeographic structuring during climate change: the signature of leading and rear edge of range shifting populations\' and by the Eco-Ethology Research Unit (PEst-OE/MAR/UI0331/2011, *Fundação para a Ciência e Tecnologia*, partially FEDER funded). The collection of fish and the monitoring census benefit from work done through the past projects: CLIPE, MAREFISH, MARÉ, MAROV, OGAMP, MARMAC, MARMAC II. SS is a research fellow supported by the Marie Curie Grant co-funded by the EU under the FP7-People-2012-COFUND; Co-funding of Regional, National and International Programmes, GA no. 600407 and the Bandiera Project RITMARE. SMF was supported by an FCT Grant (SFRH/BPD/84923/2012). IMAR-Univ Azores was funded through the pluriannual and programmatic funding scheme as research unit 531 and associate laboratory LARSyS 9, with financial support from OE and COMPETE, through FCT and DRCTC/GovAzores. We are also extremely grateful to Norberto Serpa, Azores, Marta Linde, Balearic Islands, Molly Buchholz-Sørensen, Malta, Paolo Guidetti and Marco Arculeo, Italy, for collecting samples. This work has benefited greatly from the comments of Per Erik Jorde from the University of Oslo and three anonymous reviewers.

We dedicate this paper to the memory of our dear professor, colleague and friend Vítor Almada who prematurely died during the course of this work.

[Supplementary Information](#sup1){ref-type="supplementary-material"} accompanies this paper on Heredity website (http://www.nature.com/hdy)

The authors declare no conflict of interest.

Supplementary Material {#sup1}
======================

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

![Sampling locations for *Diplodus vulgaris*. Colour-coded pie charts indicate relative contribution of mtDNA CR shared haplotypes organised by sampling areas. Private haplotypes are in orange, and the dimension of each pie is proportional to sample size.](hdy201555f1){#fig1}

![Bar plots of population pairwise *F*~*ST*~ (in black) and *D*~est~ (in grey) for each of the marker type used in this study. ATL, Atlantic; AZO, Azores Is.; CAD, Cadiz; CAN, Canaries Is.; MA, Malta; MAD, Madeira Is.; MAU, Mauritania; PA, Palermo; PM, Palma de Mallorca (Balearic Is.); PTC, Portugal Centre (Sesimbra, Tejo and Sado estuary); PTS, Portugal South (Algarve, Arade and Lagos), TG, Torre Guaceto. \*Statistical significance of differentiation was assessed through the calculation of 95% confidence limits out of 10 000 bootstrap replicates.](hdy201555f2){#fig2}

![Haplotype relationships of *D. vulgaris* represented in a median joining network for mtDNA (CR) and nuclear DNA (S7). Haplotypes are represented by circles whose sizes are proportional to absolute frequencies in the total sample and the colours used represent the locations as indicated by the legend on the side.](hdy201555f3){#fig3}

![Multidimensional scaling plot based on the matrix of *F*~ST~ pairwise comparisons from microsatellites (**a** and **b**), mtDNA CR (**c** and **d**) and S7 (**e** and **f**) data sets, respectively. PTC, Portugal Centre (Sesimbra, Tejo and Sado estuary); PTS, Portugal South (Algarve, Arade and Lagos). Dark grey dot refers to Atlantic and light grey triangles refer to Mediterranean samples.](hdy201555f4){#fig4}

![Microsatellite against mtDNA CR for the *F*~ST~ and *D*~est~ matrices using all sampling localities. A significant correlation using all data is reported only for the *F*~ST~ matrices (*r*=0.851; *P*=0.034), whereas it is not significant when comparing the *D*~est~ matrices (*r*=0.166; *P*=0.204). When data are treated separately as 'only Azores comparisons\' (dark grey dot, upper trendline) and 'excluding Azores comparisons\' (light grey triangles, lower trendline), the pairwise comparisons within each one of these two 'clouds\' also show lack of correlation.](hdy201555f5){#fig5}

![(**a**) The three most realistic ABC scenarios comparing the Azorean sample (AZ) with two other potential source populations. (**b**, left) Reliability of all scenarios under the three models (I=Cornuet--Miller; II=Beaumont; III=Guillemaud) visualised through posterior probabilities as the relative proportions (direct estimates) and logistic regressions (red, green and blue represent the parameters associated to the scenarios 1, 2 and 3, respectively). (**b**, right) Principal component analysis (PCA) showing the good fit of posterior distributions for scenario 1 under the three models (the observed data---large yellow spots---fall in the clouds of the dots representing the posterior (larger red dots) and the prior (smaller red dots) distributions). CAD, Cadiz; MAD, Madeira Is.; PTC, Portugal Centre (Sesimbra, Tejo and Sado estuary); PTS, Portugal South (Algarve, Arade and Lagos).](hdy201555f6){#fig6}

###### Sample locations, sample sizes and summary statistics for the mtDNA CR, nuclear S7 and nuclear microsatellites of *Diplodus vulgaris*

  *Basin*                        *Locality*                 *Mitochondrial CR*   *S7*  *Microsatellites*                                                                                                     
  --------------- ---------------------------------------- -------------------- ------ ------------------- ------------- ---- ---- ------- ------- ------- -------- ---- ---- ------- ------- -------------- --------
                                 Azores Is.                         23            23   0.668±0.035         0.014±0.004    15   30   0.013   0.235   0.174   0.275    28   4    0.784   0.766      0.041       7.167
                                Madeira Is.                         30            30   0.940±0.006         0.017±0.002    19   38   0.017   0.276   0.266   0.035    34   7    0.816   0.784   [0.054]{.ul}   8.911
  Atlantic                      Canaries Is.                        21            19   0.983±0.006         0.021±0.003    11   22   0.016   0.264   0.209   0.218    22   12   0.843   0.809   [0.064]{.ul}   9.939
                                 Mauritania                         10            10   0.978±0.018         0.016±0.003    4    8    0.015   0.464   0.306   0.377    27   13   0.798   0.762      0.070       8.872
                   Portugal C.[b](#t1-fn4){ref-type="fn"}           34            29   1.000±0.002         0.022±0.002    24   48   0.016   0.229   0.215   0.061    46   5    0.842   0.786    **0.078**     9.521
                   Portugal S.[c](#t1-fn5){ref-type="fn"}           25            22   0.961±0.006         0.017±0.002    26   52   0.017   0.208   0.185   0.111    60   26   0.839   0.826      0.025       10.069
                                   Cadiz                            20            20   0.984±0.006         0.013±0.002    4    8    0.013   0.514   0.600   −0.333   26   5    0.829   0.774    **0.087**     9.478
                                 Malta Is.                          11            3    1.000±0.157         0.023±0.011    5    10   0.017   0.391   0.283   0.299    11   2    0.792   0.798      0.040       9.778
  Mediterranean                 Balearic Is.                        14            11   1.000±0.012         0.017±0.003    9    18   0.016   0.366   0.333   0.094    14   7    0.829   0.857      0.003       9.642
                                  Palermo                           28            19   1.000±0.004         0.018±0.002    10   20   0.015   0.295   0.314   −0.070   28   11   0.863   0.885      −0.007      10.104
                               Torre Guaceto                        11            8    0.964±0.029         0.020±0.005    10   20   0.013   0.289   0.208   0.290    11   0    0.804   0.859      −0.020      9.556

Abbreviations: CR, control region; *F*~IS~, Wright\'s inbreeding coefficient within populations; *h*, haplotype diversity (±s.e.); *H*~e~, unbiased expected heterozygosity; *H*~o~, observed heterozygosity; mtDNA, mitochondrial DNA; *N*, number of individuals screened; *n*, number of haplotypes; π, nucleotide diversity (±s.e.); PA, number of private alleles; *R*~A~, allelic richness

Values that significantly deviate from the expected Hardy--Weinberg equilibrium are in bold, whereas underlined are reported the values that are no longer significant after performing a Bonferroni test.

Number of alleles.

Sesimbra, Tejo and Sado estuary.

Algarve, Arade and Lagos.

###### Fixation indices for the mtDNA CR corresponding to the groups (*k*) of populations inferred by AMOVA for the 11 areas of sampling

  *k*                 *Group composition*                F~*SC*~     F~*ST*~      F~*CT*~     F~*SC*~*/*F~*CT*~
  ----- ----------------------------------------------- ---------- ------------ ------------ -------------------
  2          \[AZ\] \[MAD+Atlantic+Mediterranean\]       0.0166\*   0.1219\*\*   0.1070\*\*        0.1551
  3        \[AZ\] \[MAD+Atlantic\] \[Mediterranean\]     0.0146\*   0.0657\*\*   0.0519\*\*        0.2813
  4      \[AZ\] \[MAD\] \[Atlantic\] \[Mediterranean\]    0.0022     0.0578\*    0.0558\*\*        0.0394

Abbreviations: AMOVA, analysis of molecular variance; Atlantic, Canaries Is., Mauritania, Cadiz, Portugal Centre and Portugal South; AZ, Azores Is.; CR, control region; MAD, Madeira Is.; Mediterranean, Malta, Palermo, Balearic Is. and Torre Guaceto; mtDNA, mitochondrial DNA.

*F*~ST~ (variance among populations relative to the total variance), *F*~SC~ (variance among populations within groups) and *F*~CT~ (variance among groups relative to the total variance).

\**P*\<0.05 and \*\**P*\<0.001.

###### Summary results from approximate Bayesian computation (ABC) approach obtained for the most plausible scenario

  *Populations*                             *Cornuet--Miller*            *Beaumont*               *Guillemaud*
  ------------------------------------- ------------------------- ------------------------- -------------------------
  Azores *N*~e~                              312 (65.5--474)           300 (54.3--467)           283 (54.2--467)
  Madeira+Portugal South+Cadiz *N*~e~    38 400 (20 300--48 700)   36 700 (18 000--48 400)   39 600 (22 100--48 900)
  Portugal Centre *N*~e~                  17 400 (5260--35 000)     16 100 (5060--33 700)     18 400 (5690--35 400)
  Time of divergence *t*1                   48.5 (7.9--88.9)           52.4 (6.32--89)          41.2 (5.5--83.6)
  Time of divergence *t*2                     154 (26--596)            521 (127--920)             296 (62--775)

Median effective population size (*N*~e~) estimates (and 95% confidence interval) and time of coalescence in terms of generations (*t*1 and *t*2) from each of the three models implemented in the analysis.

[^1]: These authors contributed equally to this work.
