Genetic and Morphological Variation of the Forkbeard, Phycis phycis (Pisces, Phycidae): Evidence of Panmixia and Recent Population Expansion along Its Distribution Area

The knowledge of population structure of a species is essential to effectively assess and manage fisheries. In the present study, genetics, by mitochondrial DNA cytochrome b sequence analysis, and body geometric morphometrics were used to evaluate the existence of distinct populations of the forkbeard (Phycis phycis), an important commercial species in several European countries, especially Portugal and Spain. For geometric morphometric analysis, specimens were collected in the Northeast Atlantic Ocean—Azores, Madeira and mainland Portugal, and for genetic analysis, these samples were complemented with samples collected in the Mediterranean Sea—Spain, Italy and Croatia, in order to cover the entire distribution area of the species. Body shape of the forkbeard from the Northeast Atlantic was found to be highly variable. This variation was probably associated with different environmental factors between the study areas. Despite morphological variation, a low genetic differentiation between samples from different areas was found, most likely due to gene flow that occurred in the past or with the demographic history of the species. Moreover, the presence of unique haplotypes in the Northeast Atlantic and in the Mediterranean suggests that recent gene flow between populations from these areas should be limited. Altogether, a high haplotype diversity, a low nucleotide diversity, a “star-like” network and the results of the mismatch distribution, indicate a possible signature of recent population expansions, which probably started during the end of the Last Glacial Maximum and led to the colonization of the Northeast Atlantic and the Mediterranean.


Introduction
of fragments of mitochondrial DNA (mtDNA) cytochrome b gene (cyt b) and of the first intron of the nuclear S7 ribosomal protein gene (S7). These genes have been widely used in both genetic structure of fish populations (e.g. [23][24][25]) and fish phylogeographical and demographic studies (e.g. [25][26][27][28]), proving to be useful in revealing historical and present-day barriers to gene flow in widespread fish species [29]. Variation in body shape of the forkbeard was analysed using geometric morphometrics (by landmark analysis), another tool commonly used in fish population differentiation (e.g [30][31][32][33][34][35][36]). Through both genetic and morphological data, patterns of variation between populations could be interpreted.

Ethics statement
Sampling for the present study focused on the forkbeard from the Northeast Atlantic Ocean and the Mediterranean Sea. Ethical or government approval, specific permissions or licenses for sample collection were not required for this study, as all specimens were collected as part of routine fishing procedures by fishermen of commercial fleet. As the forkbeard is a deep-water species, fish are killed during the hauling of fishing gear due to differences in atmospheric pressure. No animals were killed specifically for this study. Dead specimens were acquired in the commercial ports of Horta (Azores, Portugal), Funchal (Madeira, Portugal), Peniche (mainland Portugal), Fuengirola (Spain), Fiumicino (Italy) and Split (Croatia). Sampling was performed in the laboratory. The forkbeard is not an endangered or protected species.

Sampling
For geometric morphometric analysis, forkbeard specimens were collected from commercial landings of fishing vessels operating off Azores and Madeira archipelagos, and mainland Portugal-Northeast Atlantic (Fig 1), between November 2011 and March 2013. All specimens were deep-frozen (-18˚C), stored in a horizontal position to avoid any deformation of the body until the time of the analysis, and thawed before being analysed. Only fish with mouth closed and with no body deformation were analysed. For each fish, total length (TL, to the nearest 0.1 cm) and sex were recorded. A total of 268 adult and sexually mature fish (TL > 40 cm) were analysed: 79  For genetic analysis, a piece of fin of sampled specimens from each population from the Northeast Atlantic-Azores (n = 29), Madeira (n = 26) and mainland Portugal (n = 27) was collected and preserved in 96% ethanol. In order to cover the entire distribution area of the forkbeard, fins of specimens from 3 populations from the Mediterranean Sea-Spain (n = 17), Italy (n = 2) and Croatia (n = 26) were also collected (Fig 1), preserved in 96% ethanol and shipped to the laboratory.
To remove non-shape variation, a generalized Procrustes superimposition was applied to the landmark configurations [43]. This corresponds to fitting a very simple model, taking into account only translation, rigid rotation, scale and possibly uniform shape change for the differences in landmark configurations [43], in order to minimize the sum-of-squared distances (Procrustes distances) between homologous landmarks of all specimens. Given the different length ranges of the samples from the three study areas, a multivariate regression of the Procrustes coordinates on centroid size was performed [44] to remove the size effect and possible allometric relationships between variables. The residuals of this regression were used as "sizefree" variables in following statistical analyses.  The existence of sexual dimorphism within the study areas was determined by the quantification of differences between the mean shape of males and females using the Procrustes distance. A permutation test using 10,000 runs was used to assess the null hypothesis of no differences between shape of males and females from the three study areas.
Canonical variate analysis (CVA) was performed to detect morphometric differences in the body shape of forkbeard from the three study areas and to investigate if body shape could be used to classify samples in terms of area of origin. Pairwise comparisons of mean shapes of forkbeard from different areas were based on Procrustes distances, and a permutation test with 10,000 runs was used to assess the null hypothesis of no difference between samples [35].
The mean shapes of the forkbeard from the three study areas were estimated by warped outline drawings, which were performed using the wireframe function [45].
All geometric morphometric analyses were done in the software MorphoJ version 1.05c [45]. A significance level of 0.05 was set for all the statistical analysis.
Statistical analysis. Canonical discriminant analysis (CDA) was also performed to detect morphometric differences in body shape of the forkbeard from the three study areas. The Jackknife cross-validation was used to calculate an unbiased estimation of classification success, being each specimen allocated to the group with the nearest centroid and the proportion of correctly allocated specimens calculated [46]. The discriminatory effectiveness of the analysis was also calculated using the Wilks' lambda (λ) [47] and the Cohen's Kappa (κ) [48] (the latter provides more reliable results than the former [49]). These statistical analyses were performed in IBM SPSS Statistics 20. The significance level was set at 0.05 for all the statistical tests used.
PCR amplifications were performed in a 20 μl total reaction volume with 0.6 μM of each primer, 0.175/0.2 mM dNTPs, 0.925/1.125 mM MgCl 2 , 1.4/1.2 μl BSA (10 ng/μl), 4.0 μl 5x Colorless GoTaq Flexi Buffer, 0.05 U GoTaq DNA Polymerase (Promega) and 2.0 μl of template DNA (values for cyt b/S7, respectively). PCR conditions used were composed by an initial denaturation step at 95˚C for 5 min, followed by 35 cycles of denaturation at 95˚C for 45 sec, annealing at 55˚C or 68˚C, for cyt b or S7, respectively, for 35 sec and extension at 72˚C for 1 min, with a final extension period at 72˚C for 10 min. Amplification results were confirmed by 1% TBE agarose gel electrophoresis. All PCR products were purified with the SureClean kit (Bioline) following the manufacturer's protocol, and sequenced in the forward and the reverse directions (using the same primers) by an ABI PRISM 3700 DNA analyser at Macrogen (http://www.macrogen.com). The existence of haplotypes present in one single specimen was confirmed by performing independent PCR amplifications and sequencing (as previously described). The obtained cyt b and S7 sequences were compared against a database in Gen-Bank using the basic local alignment search tool (BLAST) (http://blast.ncbi.nlm.nih.gov/Blast. cgi) to confirm that they belong to both genes.
Fragment sequences of both genes were verified and edited using the software SEQUENCHER version 4.0.5 (Gene Codes Corporation) and BIOEDIT version 7.0.9 [51]. For each specimen, the existence of unique nucleobases for each sequence position was confirmed by analysing the chromatograms, excluding any possible occurrence of errors in the base calling and any evidence of pseudogenes. All sequences of both genes were aligned using CLUSTALX version 2.1 [52] and converted into the appropriate format with CONCATENATOR version 1.0.1 [53].
Genetic data analysis. Data analysis of nuclear S7 was suspended because this gene did not show genetic variation.
Two different geographical groups were defined a priori: Northeast Atlantic Ocean, composed by samples of Azores, Madeira and mainland Portugal, and Mediterranean Sea, constituted by samples of Spain, Italy and Croatia.
The number of polymorphic sites and haplotypes, haplotype diversity (h), and nucleotide diversity (π) were calculated for each sampled population, as well as for both geographic groups and for the total samples, using the software DNASP version 5.10.01 [54]. Maximum parsimony median-joining haplotype network [55] was created using NETWORK version 5.0.0.0 (Fluxus Technology Ltd.). Analysis of molecular variance (AMOVA), with 10,000 permutations, was performed using ARLEQUIN version 3.5.1.2 [56] to assess population genetic structure. This analysis produces estimates of variance components and F-statistic analogues, reflecting the correlation of haplotypes at different levels of hierarchical subdivision [57].
Mismatch distribution (frequency distribution of the number of pairwise differences among haplotypes) was used to detect and estimate the timing of population expansion [58]. Estimated expansion values were calculated in ARLEQUIN software and graphics of frequency distribution were performed using DNASP. Populations that have been historically stable are predicted to have multimodal mismatch distributions, whereas those that have undergone a recent expansion are predicted to be unimodal [59]. To test the observed mismatch distribution goodnessof-fit to the expansion model and to obtain confidence intervals around the estimated mode of mismatch distribution, 1,000 permutation replicates were used [60]. Statistically significant differences between observed and expected distributions were evaluated with the sum of the square deviations (SSD) and Harpending's raggedness index (raggedness) [61,62]. Neutrality tests of Tajima's D [63] and Fu's F-statistics [64] were performed in ARLEQUIN to detect changes in population size and/or estimating deviations from neutrality, assuming a constant population size at mutation-drift equilibrium [2]. Thus, evidence of expanding populations was assumed when significant negative values of Tajima's D and Fu's F-statistics were obtained [65,66].
Time since population expansion in years (t), as well as the 95% confidence interval, was calculated from τ = 2ut, where τ is the expansion time parameter generated by ARLEQUIN and u is the mutation rate per nucleotide per year multiplied by sequence length (i.e. number of nucleotides) [58,62]. A generation time of one year was assumed, as the forkbeard spawns every year [20], and a mutation rate of 3.86% per million years was considered, based on other gadiform species [67], since no information on mutation rate of the forkbeard is available.

Geometric morphometric analysis
In this study, a similar length range of specimens from the Northeast Atlantic was used, based only in adult mature specimens. However, and despite the fact that juveniles and very large specimens were excluded from analyses to minimize the effect of size on shape, a very small amount, but statistically significant (p < 0.031) of shape variation due to allometric growth was detected (1.16% of total shape variation). Therefore, data were allometrically corrected to remove the size effect between variables and the residuals were used as "size-free" variables in the statistical analyses.
No statistically significant differences were found between males and females shape within the three study areas from the Northeast Atlantic (p Azores = 0.094; p Madeira = 0.075; p Mainland = 0.268), so sexes were analysed together for each geographical area. Statistically significant differences were found in the mean body shape of forkbeard between the three geographical areas (Table 1). Fig 3 shows the differences between the mean shapes of forkbeard from the three areas and the overall mean shape. Specimens from Azores and mainland Portugal showed a relative displacement of anatomic landmarks 2, 3 and 7. A relative expansion of body height was observed in Madeira specimens. Both Madeira and Azores fish also showed a relative displacement of landmark 10 (posterior limit of the operculum).  The two first canonical variate functions were significant for the discrimination of specimens from the three study areas (p < 0.0001). Their score plot (Fig 4) shows a separation between the three areas, although some overlapping can be observed. The first two canonical variate functions explained 56.7% and 43.3% of between-group variance, respectively.
Assignment of forkbeard specimens to geographical areas correctly classified 77.6% of the total analysed specimens (λ = 0.234; κ = 0.663) ( Table 2), meaning that a low misclassification occurred between the three study areas. The highest classification rate was found for mainland Portugal specimens with 81.2% classification success, while Azores showed the lowest classification success with 74.7%.

Genetic analysis
Within the 127 samples of forkbeard from the Northeast Atlantic and the Mediterranean, a total of 16 haplotypes of the mitochondrial cyt b gene fragment (GenBank accession numbers: KM252661-72 and KX215149-52) and only one haplotype of the nuclear S7 gene fragment (GenBank accession number KP005452) were found. The results of the latter gene were not shown due to its lack of genetic variation. Pseudogenes were not found in any of the samples sequenced, since two real peaks in the same position during the base calling in the chromatograms were never obtained. The 846 bp amplified fragments of cyt b showed 17 polymorphic nucleotide positions and 8 of the 16 haplotypes were unique ( Table 3). The Atlantic group showed 4 unique haplotypes (H3, H14 and H15 from Mainland Portugal, and H8 from Azores), all derived from the most commons haplotypes, and the Mediterranean group also exhibited 4 other unique haplotypes (H9, H10 and H16 from Croatia, and H12 from Spain), all being derived haplotypes. The most common haplotype (H1) was present in 67 specimens (~53% of total samples) and was found in all sampled geographic areas. The overall values of haplotype diversity (h) and nucleotide diversity (π) for the entire sample were 0.675 ± 0.038 and 0.00144 ± 0.00013, respectively. The Northeast Atlantic and the Mediterranean groups showed similar haplotype and nucleotide diversities (0.626 ± 0.054 and 0.00133 ± 0.00016; 0.749 ± 0.046 and 0.00165 ± 0.00022, respectively). In the Atlantic group, the highest haplotype richness was found in Mainland Portugal, where, for 27 samples, nine different haplotypes were present, and Azores showed the lowest one (5 different haplotypes were detected in 29 samples). In the Mediterranean Sea, the highest haplotype richness was found in Croatia, with 8 different haplotypes in 26 samples. A summary of the distribution of haplotypes per population and group, as well as a summary of haplotype and nucleotide diversities indices, are presented in Table 3. The haplotype network did not show a phylogeographical structure among the sampled geographical region, presenting a "star-like" shape where the most common haplotype (H1) is shared between all sampled geographic areas, and the less frequent and unique haplotypes are connected to the central most common haplotype by only one or two mutation steps (Fig 5).
To assess the distribution of mtDNA variation, sampled geographic areas were pooled according to the two defined regional groups (Northeast Atlantic Ocean and Mediterranean Sea) and AMOVA was performed. AMOVA also revealed a lack of genetic structure, indicating that the overall source of variation was within populations (100.95%, p = 0.721) instead of among groups (0.43%, p = 0.301) ( Table 4). The pairwise F-statistics (Table 5) showed no significant differentiation between the two regional groups.
The distribution of pairwise nucleotide differences (mismatch distribution) showed a unimodal shape (Fig 6), suggesting that population range expansion may have occurred. The observed raggedness index was low and both P SSD and P RAG showed that the observed   distributions did not differ significantly from those expected under a sudden expansion model ( Table 6). Tajima' D and Fu's F tests were negative, although not statistically significant (Table 6), which do not support the population expansion hypothesis. The timing of demographic expansion for the forkbeard was estimated to have occurred at approximately 17,245 years ago (5,298-32,731 years ago).

Population structure and morphological differentiation
A substantial variation in forkbeard body shape was found, suggesting a clear separation between the three populations from the Northeast Atlantic. However, this morphological differentiation is not in accordance with the genetic population structure evidenced by the molecular marker used. A large component of morphological variation in fishes is induced by exposure to different environmental factors during their development [68], namely food availability, temperature, salinity and prolonged swimming [33,69]. Thus, the phenotypic plasticity of organisms may determine the phenotypic differentiation in populations of a species experiencing specific environmental conditions [33,69]. As fish morphology is particularly dependent on   Table 6.
doi:10.1371/journal.pone.0167045.g006 environmental factors during early life stages, phenotypic differentiation may indicate that most fish from each population spent their lives in separate regions [33,70] and provide an indirect assessment of population structure. However, it does not give direct evidence of genetic isolation between populations [35,70]. Body shape and other morphological characters have long been used to delineate populations and continue to be used successfully [30][31][32][33][34][35][36], based on the fact that morphological features can relate to fitness and strong selective pressures that may underlie rapid genetic divergence between groups of fish [10].
Results of body geometric morphometric analysis obtained in this study suggests its usefulness as a population differentiation tool and indicate the existence of at least three phenotypically distinct local populations in the Northeast Atlantic, as reported by Vieira et al. [19] using otolith shape analysis. This population differentiation suggests a link between the levels of phenotypic divergence and geographical distance, and also a limited adult fish migration between populations from the Azores, Madeira and mainland Portugal. It is possible that different conditions (e.g. sea water temperature, depth, sediment type, food resources, predation risk, etc.) between the three study areas affect fish growth and behaviour [71], and thus, fish morphology. The discriminant analysis indicated that morphometric differentiation between the Azores and mainland Portugal was mainly due to differences in head characters of fishes, which may reflect differential habitat use [33] and may also be related to prey size [72,73]. Madeira specimens showed a change in body depth that are mostly related to hydrodynamic and swimming performance [74][75][76], and locomotor adaptations to food capture and escape from predators [75,76]. It is known that a higher body depth in benthic fish provides them a steeper power curve, furthering the prey capture and escape from predators (e.g. [74,77,78]).
Genetic data from cyt b suggest that the forkbeard forms a panmictic population along the Northeast Atlantic and the Mediterranean. F-statistics and AMOVA indicated no genetic structure over space, probably due to the high frequency of the most common haplotype (H1) (which minimizes the effect of rarer haplotypes), and also due to similarities between sequences caused by the low nucleotide diversity. These facts suggest that high levels of gene flow between geographically separated areas must have occurred in the past, or that the current distribution is an expanded area of an ancestral refugial and panmitic population, with low levels of diversity but widespread common/ancestral haplotypes. Moreover, and taking into account the existence of unique haplotypes in the Northeast Atlantic and the Mediterranean, recent gene flow between populations from these groups is unlikely to occur or, if present, occurs at low levels. Although larval behaviour of the forkbeard is unknown, other species of the genus Phycis present pelagic larvae for several months [79,80] and, thus, long distance dispersal may occur during forkbeard larval phase. The putative absence or reduction of current larval dispersion of the forkbeard between the Northeast Atlantic and the Mediterranean  suggests the presence of significant barriers to gene flow in this species. In fact, oceanographic and topographic features of the Atlantic, such as the patterns of circulation flow and water masses around seamounts, cause parcels of water to become trapped around these structures (e.g. Taylor columns) and retain larvae [81]. In addition, the large geographical distance between the European continental slope, the Azores (Mid-Atlantic Ridge) and Madeira may also create a barrier to dispersion of larvae [24] of this species.

Demographic history
The genetic diversity pattern of the forkbeard might be consequence of population expansion after a period of low effective population size due to founder events or bottlenecks [82]. The phylogenetic relationships between mtDNA cyt b haplotypes was defined by a "star-like" network topology, in which rare haplotypes are derived from the most common haplotypes (H1, H2 and H6), presumably ancestors, frequently by one or two mutation steps. This could indicate that the forkbeard populations have recently expanded in size from one or from a small number of founders following a population bottleneck [59], suggesting a scenario of sudden population expansion in the Northeast Atlantic and the Mediterranean. Furthermore, mismatch distribution analysis also supports the hypothesis of population range expansion. The time of expansion estimated for the forkbeard was suggestive of a growing population since the end of the Last Glacial Maximum (LGM). Such demographic changes might have been caused by variations in the sea temperature and level during this time. Previous studies have suggested a strong historical influence on the genetic population structure of other gadiform fish species in the North Atlantic has (e.g. [83,84]). One of the major focuses in studies of phylogeography of the Northeast Atlantic specimens has been the location of refugia, where these organisms survived the successive glacial peaks during the Pleistocene [27]. These refugia most likely acted as propagation sources for the recolonization of more northerly areas during the interglacial periods [85]. Most studies consider the Pleistocene and in particular the LGM to have been critical times when populations were not present in the frozen north and went through severe bottlenecks in their refuge habitats [86]. The forkbeard has a subtropical geographical distribution, exhibiting a reduced temperature tolerance (12-15˚C,~300 m depth). During the last glaciation and due to the progression of the polar front, this species most likely has been pushed to southern limits of its distribution to seek refugia with similar temperature to its optimum. Based on this, and on the haplotype diversity presented by the forkbeard populations and the generated network, two possible scenarios of population expansion could be delineated. In the first scenario, during the LGM, forkbeard populations took refuge in the Mediterranean Sea and, after the end of this period, the species expanded to the northern and southern Northeast Atlantic. However, and according to this scenario, the Mediterranean group should show much higher haplotype diversity than Atlantic one. On the other hand, the Mediterranean LGM sea temperatures decreased drastically, being low enough to support populations of northern gadids and Atlantic salmon [85], but never to support forkbeard populations. Thus, this scenario is unlikely to have occurred. In the second and most likely expansion scenario, during the LGM the southern extreme of forkbeard distribution (Cape Verde region) served as refugium for haplotypes H1, H2 and H6, and, after the glaciation, the species expanded northward colonizing the rest of the Northeast Atlantic and the Mediterranean. The presence of the most common haplotype (H1) in all sampled geographic areas corroborates this hypothesis. Then, and although the gene flow between the Atlantic and the Mediterranean groups may have occurred, each group began its own process of divergence as suggested by the presence of four unique haplotypes in both the Atlantic and the Mediterranean groups. In fact, the oceanographic and topographic features of the Atlantic [81], and the large geographical distance between the Azores, Madeira and the European continental slope [24], seems to have limited the gene flow between the Northeast Atlantic and the Mediterranean populations of the forkbeard, leading to the emergence of unique haplotypes in each region. This hypothesis of the southern extreme refugium can also be reinforced by the similarity of sea temperature in the Western Sahara coast, Cape Verde and Mauritanian coast during the last glaciation [87,88], and the current forkbeard's thermal belt. With the rising of sea level in interglacial periods, new areas could be colonized, enabling the forkbeard to extend its ranges and resulting in spatial and demographic and expansions. In agreement, the Atlantic tropical cost of Africa has previously been proposed as possible refugium during the glacial periods for fish species of the Northeast Atlantic (e.g. bleniids [26,89], Chromis limbata [90]).

Final considerations and implications for fishery management
At the present, little information exists about the population structure of the forkbeard and from a management point of view the knowledge of this issue is crucial.
Genetic data suggest that the forkbeard forms a panmictic population in the Northeast Atlantic and the Mediterranean, and body shape variations indicate, as otolith shape analysis had already shown [19], a clear separation between the three populations from the Northeast Atlantic. The absence of a clear genetic structure and the differences in forkbeard morphology are not contradictory, since neutral genetic markers were used in this study. In fact, similar results have already been obtained by other authors studying different animal groups like molluscs (e.g. [2]) and fish (e.g. [34,91,92]). However, the existence of unique haplotypes in the Northeast Atlantic and the Mediterranean suggests that recent gene flow between populations from these groups should be limited. These facts suggest that gene flow between geographically separated areas must have occurred in the past, probably during the larval phase of the species, being more restricted in the present. Furthermore, a prolonged separation of post-larval fish in different environmental regimes leads to phenotypic plasticity in response to distinct selective pressures or environmental constraints.
The information presented in this study shows that, along its distribution area, the forkbeard presents, at least, three phenotypic stocks-mainland Portugal, Azores and Madeira, in the southern NE Atlantic, and two genotypic stocks between the NE Atlantic and the Mediterranean. Any depletion in one of these populations might not be compensated by migration from the others, at least at a sufficiently rapid rate to ensure resource sustainability. The most precautionary management approach should be to consider the different populations of the forkbeard from the Northeast Atlantic and the Mediterranean as separate stocks, to ensure sustainability of resources and genetic biodiversity maintenance [13]. A failure to employ the knowledge on the population structure in the management of complex fish stocks, such as in the case of the forkbeard, may originate incorrect management actions that may lead to overexploitation of some stocks.
Future studies using a holistic approach, which employ a broad spectrum of complementary techniques (e.g. otolith chemical composition and parasites as biological tags) [13], as well as the use of more sensitive markers such as microsatellites, or adaptive markers like candidate genes (e.g. quantitative trait loci associated with body shape, such as genes coding the a-3 subunit of Collagen VI [93] and a glycerol kinase [94]), will be necessary to complementary address this issue.
respectively. Authors are very grateful to the people who helped in the collection of fins from fishes from the Mediterranean Sea, particularly to Dr. José Luis Rueda (IEO, Spain), Dr. Luciana Sola (UniRoma, Italy) and Dr. Sanja Matić-Skoko (IOF, Croatia).