Sample collection & DNA extraction

Between 2002 and 2009 a total of 250 E. marginatus tissue samples were collected within the central Mediterranean region by line, spear and net, and for those recorded, between 1 and 64 m depth (Fig 1; S1 File). Samples from Croatia, Libya, N. Sicily, IT and Linosa Island, IT were collected by boat or directly from the reef. Samples donated from Tunisia were collected from fishermen at Mercabarna Fish Market in Barcelona, ES and are a sub-set used in a previous study by Schunter et al. [35]. In accordance with conservation of an endangered species, it was determined prior to commencement of this project that no Maltese fish would be killed for use in this local study. Maltese samples were collected from the Valletta Fish Market, artisanal fishermen or from local restaurant owners. Between 2007 and 2009, a total of 89 dusky grouper tissue samples were collected within Malta’s FMZ for which 31 samples have known catch location sites (S2 File). For all samples, GPS coordinates were provided by the collector or determined afterward based on catch location site name provided. Portions of skeletal muscle, organ, gill and/or caudal fin from both caught and live fish were collected and preserved in 20% DMSO, 250 mM Na₂-EDTA saturated in NaOH or 96% EtOH and/or the -80°C freezer [42, 43].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Epinephelus marginatus sampling sites and suitable habitat map within the Mediterranean Sea.

Collection sites of 250 E. marginatus individuals in the central Mediterranean region. Origin of samples and spawning season (summer) surface circulation [44] are overlain on an Aquamaps suitable habitat map for E. marginatus where yellow to red colours represent least to most suitable habitat [45]. Hurd Bank, MT (HB), Malta, MT (MAL), Linosa, IT (LIN), Croatia (CRO), Libya (LIB¹, LIB², LIB³), N. Sicily, IT (SCL) and Tunisia (TUN).

https://doi.org/10.1371/journal.pone.0159864.g001

Total genomic DNA was isolated from individuals either by phenol:chloroform extraction [46] or using a Qiagen ‘DNeasy Blood and Tissue’ kit (Qiagen Cat. No. 69504). Both methods consumed between 25–200 mg of tissue. Extracted genomic DNA by phenol:chloroform was re-suspended in 50–100 μl TE buffer (10 mM Tris pH 8.0, 1 mM EDTA) whilst samples processed by Qiagen ‘DNeasy Blood and Tissue’ kit were processed according to the manufacturer’s protocol for tissue and eluted in 100 μl AE buffer (10 mM Tris pH 9.0, 0.5 mM EDTA).

DNA analyses

Fourteen microsatellite primer pairs were used for intra- and interpopulation genetic analyses. The cross-species primers used were as follows: seven developed for the red hind grouper, Epinephelus guttatus (RHCA001, RHCA002, RHCA003, RHCA004, RHCA007, RHCA008 and D076 [47]), five developed for the gag grouper, Mycteroperca microlepis (GAG007, GAG010, GAG038, GAG045 and GAG049 [48]), one developed for Epinephelus merra (EM10 [49]) and one developed for Serranus cabrilla (SC06 [50]).

Microsatellite repeat motifs for each locus were determined by sequencing a random Maltese individual using a Roche 454 Genome Sequencer FLX (Roche Life Sciences, USA) at ecogenics GmbH, Zürich, CH. Annealing temperature optimization (T_a) for each locus was empirically determined using the temperature gradient option (57°C ± 5°C) on an Eppendorf Master Cycler thermal cycle (Table 1). Each 50 μl PCR reaction (10 mM Tris-HCl (pH 8.3), 50 mM KCl, 1.5 mM MgCl₂, 0.2 mM dNTPs, 0.2 μM F/R Primers, 1 U Taq (NEB Cat. No. E5000S), 0.15 mg/ml BSA, 50–200 ng genomic DNA) was amplified with the following protocol: an initial denaturation step of 94°C for 2 min, followed by 34 cycles of (94°C for 30 s, T_a 30 s, 72°C for 60 s) and a final extension step of 72°C for 15 min.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Reference data for 14 Epinephelus marginatus microsatellite loci.

https://doi.org/10.1371/journal.pone.0159864.t001

DNA amplification for microsatellite genotyping were carried out in three multiplex reactions using a Qiagen Type-it Microsatellite PCR kit (Qiagen Cat. No. 206243). Each 11 μl reaction consisted of 5.5 μl master mix (3mM; MgCl₂ [6 mM] pH 8.7, dNTPs, HotStarTaq Plus Taq), 1.1 μl primer mix (0.2 μM of each primer), 2.9 μl nuclease free H₂O and 1.5 μl genomic DNA (50–200 ng). In order to distinguish between microsatellite markers downstream during genotyping the loci pooled in each reaction differed in allelic range and 5’ flourochrome labelling in the forward direction with either 6FAM, NED, VIC or PET (Standard Dye Set DS-33; Applied Biosystems Cat. No. 450056). Using an Eppendorf Master Cycler thermal cycler the optimal PCR cycling parameters for all multiplex reactions were determined to be as follows: an initial denaturation step of 95°C for 5 min, followed by 30 cycles of (95°C for 30 s, 52°C for 90 s, 72°C for 30 s) and a final extension step of 60°C for 30 min.

Genescan^™ 500 LIZ size standard (Applied Biosytems, Inc) was added to the multiplex reactions and each sample was genotyped using a Roche 454 Genome Sequencer FLX (Roche Life Sciences, USA) at ecogenics GmbH, Zürich, CH. In addition, a small preliminary batch (30 samples with six loci) were also genotyped using an Applied Biosystems^™ 3130x1 Capillary Electrophoresis (CE) system at the Royal Zoological Society of Scotland, Wildgenes Laboratory. Several of the samples from the preliminary batch were re-run on the Roche 454 Genome Sequencer FLX for comparison and the raw allele lengths were found to be nearly (± 0.1 bp) identical and so it was decided that no adjustments between datasets needed to be made.

Scoring of alleles were done according to the quality-assurance procedures recommended by DeWoody et al. [51] and Matschiner & Salzbuger [52] where automated allele calling is combined with visual inspection of the electropherograms (FSA files) for each sample followed by automated allele binning. Alleles were sized relative to the internal size standard GeneScan^™ 500 LIZ with the aid of GeneMarker^® software ver. 1.97 (Softgenetics, LLC). The raw allele length was then converted to an integer using the parametric string operative software Tandem [52]. Micro-Checker software ver. 2.2.3 [53] was used to detect genotyping errors associated with microsatellite analysis such as stutter bands, large allele dropout, deviation from repeat motif and typographic errors. INEst (Inbreeding/Null Allele Estimation) IIM program with 10,000 iterations of the Gibbs Sampler was implemented to estimate null allele frequencies to avoid upward bias estimates seen in small and/or isolated populations such as the Maltese and Pelagie archipelagos [54].

Malta FMZ: Demographic stock assessment

Minimum and maximum weight and length ranges reported were based on the observed weight-length trend line equation obtained in this study. Because no fish were killed or purchased for conservation reasons, histological analysis was not possible and therefore age (i.e. counting sagittal otolith rings) and sex allocation (i.e. gonadal histology) information for the Maltese dusky grouper was inferred based on age- and sex-length relationships from individuals in the neighbouring Sicilio-Tunisian Channel. Sex allocation of the Maltese dusky grouper was surmised from total length (cm) data, by comparing the gender-length relationship from a histological study conducted in the Sicilio-Tunisian Channel by Chauvet [18], where based on observed data the relationship parameters were described as: juvenile (≤ 37 TLcm), female/transitional-sex (38–67 TLcm) and male (≥ 68 TLcm).

Similarly, age distribution of the Maltese dusky grouper were extrapolated from age-length relationship equations of Sicilio-Tunisian Channel individuals where x is age described in years and y is total length (cm). Chauvet [18] determined the age-length relationship of dusky groupers from the Sicilio-Tunisian Channel for individuals ≤ 21 years old with the von Bertalanffy growth curve (y = -0.1817x² + 7.9393x + 11.544, R² = 0.998), and then by the mean extrapolated linear fitting of the von Bertalanffy growth curve for years between 12 and 22 (dL/L) to produce estimates for the ages between 22 and 35 years (y = -0.0466x² + 3.5335x + 44.476, R² = 1). Based on the observed length-weight relationship of the Maltese dusky grouper in this present study, the length of Maltese individuals were then inserted into the growth equations provided by Chauvet [18] and solved with the quadratic formula to produce age estimates.

Microsatellite data analyses

Population diversity indices [55,56] and genetic variability parameters were calculated for microsatellite genotypic data using Arlequin software ver. 3.5.2 [57]. Genotypic Linkage Disequilibrium (LD) between pairs of loci were determined by a likelihood-ratio test [58–60] under the assumption of Linkage Equilibrium (LE) after applying the False Discovery Rate (FDR) method to account for multiple comparisons [61,62] with 16,000 permutations. Deviation from Hardy-Weinberg Equilibrium (HWE) was assessed using a modified version of the Markov Chain (MC) algorithm, analogous to the Fisher’s exact test, with a burn-in of 100,000 and 1,000,000 MC iterations [63,64]. The Arlequin programme was also used to carry out population specific inbreeding indices, F_IS (range: -1 to 1), where positive results were interpreted as inbreeding and negative results as outbreeding populations when compared with HWE estimates [65–67]. The effective historical population size (N_e = θ/4μ) was inferred with θ_H [68] for microsatellite data using the mutation rate (μ = 5.56X10^-4) for the common carp, Cyprinus carpio [69]. This was followed by an estimation of contemporary effective population size per generation (N_e) and effective number of breeders in one reproductive cycle (N_eb) within the FMZ using the Burrows method based on linkage disequilibrium of genotypic data [70] and the molecular coancestry method [71], respectively with the software N_eEstimator ver.2 [72].

To elucidate population subdivision of dusky groupers Geneland ver. 4.0.5 [73,74] package within R ver. 2.8.1 statistical software [75] was implemented to visualise spatial analysis of molecular clustering inferred by stochastic simulation and Markov Chain Monte-Carlo (MCMC) inference of structure using multilocus data of georeferenced individuals. Population models were tested between 1 and 10 based on a matrix of genotypes and spatial coordinates for each sample. Final parameters for spatial Bayesian analysis were set to 100,000 MCMC iterations, thinning 100 under the assumption of the Falush (correlated allele frequencies) and null allele models and repeated 10 times to verify consistency of the K estimate. Maximum rate of Poisson process (equal to sample size) and maximum number of nuclei in the Poisson-Voronoi tessellation (three times the sample size) were adjusted from the minimum default values when applicable based on recommendations from Guillot et al. [73]. Results were ranked by Geneland based on the mean logarithmic posterior probability and post-process analysis was conducted on the top run with a burn-in length period of 200. To test robustness of the results, several input parameters were varied (correlated/uncorrelated allele frequency models, null allele model, coordinate uncertainty matching monthly displacements from 20–500 m) in multiple parallel runs. Alternate models did not significantly alter results, however not enabling the ancillary null allele algorithm, slightly reduced consistency across runs inferring K under the Falush model but not under the Dirichlet model (uncorrelated allele frequencies). Management Units (MU) were assigned to the FMZ based on criteria set forth by Mortiz [76] based on significant divergence of nuclear alleles and geographic separation between groups.

In order to obtain information about the relationship between Maltese dusky groupers within the Fisheries Management Zone and neighbouring metapopulations, samples from surrounding Mediterranean locations were collected and assessed with 14 microsatellite markers. Based on recommendations by Kalinowski [77] and Pruett & Winker [78], unless otherwise indicated, a subset of up to approximately 25 individuals from each of the metapopulations were utilised in the evaluation of regional connectivity. Population statistical analyses of endangered species can be complicated due to limited population and sample size. Therefore as indicated within the methodology, various techniques of resampling were implemented to mitigate bias as a result of small or unequal sample and population sizes. The Arlequin programme was used to carry out analysis of population genetic structure inferred by AMOVA with 16,000 permutations [65–67] and an exact test of population differentiation, analogous to a Fisher’s exact test, with 100,000 MCMC and 10,000 dememorization steps with a significance level of 0.05 [79,80]. Comparison of genetic distances and coancestry within and between populations due to structure were measured by variance in allelic frequencies with Weir and Cockerham’s θ_ST, an analog of Wright’s F_ST that accounts for unequal sample sizes [67,81], Slatkin’s linearized F_ST [82] and Reynolds’ F_ST [83] with 10,000 permutations. Fixation index values (range: 0 to 1) were generally interpreted according to guidelines discussed by Wright [84] as: 0–0.05 (little genetic differentiation), 0.05–0.15 (moderate genetic differentiation), 0.15–0.25 (great differentiation) and greater than 0.25 (very great differentiation). Arlequin was also used to estimate relative ancestral population size [85], population divergence times (τ) and to produce a pairwise matrix of differentiation based on Nei’s standard genetic distance [86]. Pairwise population genetic differentiation with Jost’s D along with 95% confidence intervals [87] were calculated within R ver. 3.5.4 [75] diveRsity package [88] with 1,000 bootstrap matrices using the full data set.

Microsatellite loci under selection were assessed using 50,0000 coalescent simulations to generate joint null distributions of F_ST versus heterozygosity using the finite island model [89,90]. Loci that fell outside of the 99% quantile were identified as being under positive selection. Allelic richness (A_R) was used to estimate long-term potential for adaptability of populations with the software HP-Rare [91,92] using the hierarchal rarefication method to account for variation in sample size. Significance of results were tested using a one-tailed Sign test (α = 0.05) in order to identify populations with relatively reduced A_R and diminished capacity to withstand deterministic and stochastic events.

Demographic history was evaluated using the Garza-Williamson modified index (M value) for microsatellite data with Arlequin, for detection of population size reduction expressed as a ratio of number of alleles to allelic range, which decreases proportionally along with severity and length of population size reduction. This method was chosen because it is capable of detecting ancestral population decline several hundred generations before present time with high statistical power. Critical value (M_C) thresholds for population decline within each deme were independently calculated using Critical_M [55] with 10,000 replicates based on ancestral θ_H [68]. As a general guideline, M values greater than 0.8 are interpreted to mean there is no population size reduction, between 0.43 and 0.7 to be indicative of a recent population size reduction and anything less than 0.43 as specific to a remnant population and of a significant size reduction in the past [55]. The Bottleneck ver. 1.2.02 [93] programme was then utilised to detect evidence of recent bottleneck events within each population based on a significant difference between expected heterozygosity from the observed number of alleles and expected heterozygosity at mutation-drift equilibrium, under the assumption that a reduction in allelic diversity precedes a reduction in heterozygosity within populations that have recently experienced a decline in N_e [94]. Bottleneck is capable of detecting a population size reduction of 50 N_e within the last 25–250 generations after initial population decline. Based on recommendations for analysis of microsatellite data [95], parameters were set to 95% Stepwise Mutation Model (SMM) within the Two Phase Model (TPM, which allows for multiple-step mutations) with a TPM variance of 12% and 10,000 iterations [96]. Significance of heterozygous excess (α = 0.05) was determined by the one-tailed Wilcoxon sign-rank test [97] due to reduced probability of Type I errors when compared with alternative methods [98]. Additionally, allelic frequency mode-shift distortion from mutation-drift equilibrium (i.e. deviation from L-shape distribution) was also examined to corroborate detection of recent bottleneck events [99].

The evolutionary relationship between neighbouring Mediterranean dusky grouper metapopulations was visualised by constructing a phenogram with the programme Poptree2 [100] to generate an unrooted Neighbour-Joining [101] phylogenetic tree using δμ² genetic distances for microsatellites [102] with 1,000 bootstrap replicates. Genetic barriers to gene flow between metapopulations were investigated using Monmonier’s maximum-difference algorithm using the programme Barrier ver. 2.2 [103]. Barrier compares genetic and geographic distances under the assumption of gene flow break to produce a Delaunay triangulation where plausible genetic barriers identified by Monmonier’s algorithm are overlaid. Significance of genetic partitions were measured using bootstrap values based on 1,000 permutated θ_ST distance matrices [67], to account for variation in sample size when using the full data set with diveRsity. Mantel tests [104] were performed with 10,000 random permutations to test for statistically significant associations between pairwise genetic (F_ST/1 –F_ST) and geographic (log km) distance matrices with the software genepop ver. 4.2 [105,106]. Geographic distance between populations were measured by following the principal surface currents in June during the Mediterranean spawning season [44].

Temporal migration trends were assessed by comparing estimates of modern and ancestral gene flow. Contemporary detection of F₀ immigrant individuals within each population were directly measured with the Bayesian assignment method of Rannala and Mountain [107] using the Monte-Carlo resampling algorithm of Paetkau et al. [108] and a threshold value of 0.01 with the programme GeneClass2 [109]. This was followed by analysis with the divMigrate function within the diversity package for detection and visualisation of contemporary migratory patterns and relative magnitude with G_ST [110] based on the method of Sundqvist et al. [111]. Ancestral estimates of migration were inferred using Maximum-Likelihood (ML) simulations of coalescent genealogies with the programme migrate-n [112]. migrate-n allows for examination of past asymmetrical migration with unequal population sizes approximately 4N_e generations before present (estimated from this study to be several thousand generations). Ancestral directional migration was calculated under the Brownian microsatellite model with relative variable mutation rates estimated from the data for each loci. As recommended by Beerli [113], θ_i (mutation-scaled effective population size) and directional M_i→j (mutation-scaled effective immigration rate) were initially approximated based on 10 short and three long runs with default production parameters. Convergence analysis was then repeated four additional times using start values from the prior run to produce the final estimates. The number of effective migrants per generation (Nm) were calculated from the product of θ_i and M_i→j for diploid data (θ = 4N_eμ) where a generation is defined as 6.3 years for the dusky grouper [16].

Malta FMZ: Demographic stock assessment

Demographic results based on ancillary data collected during this study found the observed weight-length relationship of the Maltese dusky grouper was y = 0.0791x^2.63 with a correlation coefficient of R² = 0.95 (Fig 2). The length range of commercially landed dusky groupers within Malta’s FMZ between 2007 and 2009 were found to be between 23 and 103 TLcm with a mean length of 52 ± 17 TLcm (n = 84). Weights ranged between 0.3 and 16 kg with an average weight of 3 ± 3 kg (n = 41). Based on total length it was estimated that approximately 17% of the population was juvenile, 68% female/sex-changing and 15% were male with a male to female sex ratio of 1:5. Extrapolating from age-length data this would correspond to an age range between 0⁺ and 24 years, with an average age of 6 ± 4 years.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Maltese Epinephelus marginatus weight-length relationship.

Weight-length relationship of the Maltese E. marginatus based on individuals collected during this study between 2007 and 2009.

https://doi.org/10.1371/journal.pone.0159864.g002

Malta fisheries dusky grouper catch landings data (kg) is available through the Ministry for Resources and Rural Affairs (MRRA), Department of Fisheries for the years between 1947 and 2009 [114]. Whilst data for the years between 1947 and 1952 are listed under the English name of dusky perch (a synonym for the dusky grouper) and in Maltese as Ċerna, most recently, since 2000, the MRRA lists Ċerna commercially as an Epinephelus spp. or grouper spp. therefore raising concern that the landings data collected during the last decade may include additional grouper species.

The data revealed an overall decreasing trend in annual catch landings (R² = 0.69) of Ċerna in Malta between 1947 and 2009. When excluding outlier data in the years where the catch report (kg) was zero (1985, 1994–1997, 2000), estimates showed a decrease of 99% in catch landings between 1947 and 2009 (Fig 3). The monthly mean in Ċerna catch landings, averaged between 1957 and 2009, showed a peak in catch landings during the local summer spawning months between May and August (Fig 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Malta annual Ċerna catch landings data (kg) between 1947 and 2009.

Results show a decreasing trend in catch landing of Maltese Ċerna (kg) over this time period. An estimate based on catch landings from the 1940s (mean 11,583 kg year^-1) and the 2000s (mean 982 kg year^-1) show a decrease in catch landings of 92% over the last 62 years. An assessment of catch landings based on the trend line (polynomial, order 2) estimates a decline of 99% between 1947 and 2009. Outlier data (0 kg year^-1 reported) are excluded from the chart and estimates.

https://doi.org/10.1371/journal.pone.0159864.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Malta mean monthly Ċerna catch landings data (kg) between 1957 and 2009.

Mean monthly catch landings, averaged between 1957 and 2009, show Ċerna landing peaks (polynomial, order 2) during the summer spawning months between May and August.

https://doi.org/10.1371/journal.pone.0159864.g004

Malta FMZ: Genetic population structure

Within Malta’s FMZ, results show all 14 microsatellite loci highly polymorphic with the exception of RHCA003 and GAG007 which displayed two and three alleles, respectively. A total of 225 alleles were found in the 89 local individuals sampled. The number of alleles per locus ranged between 2 (RHCA003) and 27 (GAG038, GAG045), averaging 16.07 across all loci. No evidence of large allele dropout or null alleles were detected. Overall observed heterozygosity (H_O = 0.683) was lower than the expected heterozygosity (H_E = 0.757). The total effective number of adults was found to be around 130 individuals with an estimated nine (95% CI: 6.5–11.8) effective breeders per reproductive cycle (Table 2). Assessment of historical (N_e = 982) versus modern effective population size infers an eightfold decrease over the last ~4,000 generations, corresponding to roughly 25,000 years before present. Overall, genetic analysis describes the population of E. marginatus’ within the Fisheries Management Zone as decreasing in size (Ɵ_H = 2.2), which has gone through a significant size reduction in the past (M = 0.41, M_C = 0.76), with a positive and significant inbreeding coefficient (F_IS = 0.10, p < 0.001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Contemporary effective population size estimates of Epinephelus marginatus within Malta’s Fisheries Management Zone.

https://doi.org/10.1371/journal.pone.0159864.t002

Spatially explicit estimates with Geneland detected two (K = 2) unique genetic clusters within the FMZ for the basis of further genetic analysis. Post-process analysis revealed “Malta” to be a single congruent population (81% of individuals assigned) and identified a significant geographic barrier to gene flow (F_ST = 0.066, p < 0.001) between the Maltese archipelago and “Hurd Bank” (19% of individuals assigned), located 12 nm offshore on the Malta Plateau, a shallow submarine ridge connecting Malta to southern Sicily, IT (Fig 5). Global AMOVA as a weighted average over loci found 6.24% variation amongst populations, 8.51% variation amongst individuals within populations and 85.25% variation within individuals. Further analysis describes the Hurd Bank subpopulation–although results were nonsignificant–as an outbreeding population (F_IS = -0.161, p > 0.05) and Malta as an inbreeding subpopulation (F_IS = 0.091, p < 0.001), indicating a high level of self-recruitment within the island complex (Table 3). No loci were found to be under positive selection (p > 0.01) suggesting the difference between subpopulations is likely due to limited gene flow as opposed to localised adaption within an ecological niche. Based on significant genetic and geographic separation between Hurd Bank and Malta, for the purpose of conservation management, two MUs were assigned within the FMZ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Spatially explicit estimates of Epinephelus marginatus genetic clustering within Malta’s Fisheries Management Zone using Geneland.

Map of posterior probability of individuals belonging to a genetic cluster. Geneland maps assign pixels to population clusters (K = 2) where black dots represent individual E. marginatus samples. Contour lines depict the spatial change in population assignment probability where the highest probability is shown in white and the lowest in red. Map is of posterior probability of belonging to the “Hurd Bank” cluster with an overlay of the Maltese Islands. Not shown is the inverse map of posterior probability of belonging to “Malta.”

https://doi.org/10.1371/journal.pone.0159864.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Epinephelus marginatus genetic population diversity indices within Malta’s Fisheries Management Zone.

https://doi.org/10.1371/journal.pone.0159864.t003

Regional connectivity

A total of 228 alleles were found in the 117 individuals sampled from seven Mediterranean populations with a per locus range between four (GAG007) and 29 (GAG045), averaging 16.3 alleles across all loci (S3 File). All 14 microsatellite loci were highly polymorphic with the exception of GAG007 and RHCA03 which displayed only four and five alleles. The mean observed heterozygosity (H_O) ranged from 0.634 in Tunisia and 0.795 for the Hurd Bank subpopulation (S4 File). No null alleles were detected within any of the populations tested.

Global AMOVA as a weighted average over loci detected 2.1% variation amongst populations, 89.1% variation within individuals and 8.8% variation amongst individuals within populations. Results of global F_ST were nonsignificant indicating weak genetic structure of dusky grouper demes within the central Mediterranean region (F_ST = 0.021, p > 0.05). Furthermore, a global test of non-differentiation found populations to be non-differentiated (p > 0.05) from one another and analysis with the finite island model did not detect any loci under positive selection indicating reduced levels of genetic differentiation within the greater region (Fig 6). Six populations showed some deviation from HWE (between one and five loci) due to heterozygote deficit but not consistently at the same locus. Linkage disequilibrium after FDR correction was observed between one loci pair in the Linosa population and at six locus and between eight loci pair within Hurd Bank. No LD was detected within any other Mediterranean population. Linkage disequilibrium can be an effectual predictor of demographic and selective history within a population. In general, populations that are small or have experienced a recent bottleneck demonstrate a high level of LD due to genetic drift [115]. Likewise, negative selection associated with inbreeding depression contributes to lower levels of haplotype diversity, departures from HWE and smaller effective population sizes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Detection of microsatellite loci under selection in Mediterranean Epinephelus marginatus based on F_ST genome scans with Arlequin.

None of the 14 loci were shown outside the 99% quantile, the set limit for detection of directional selection.

https://doi.org/10.1371/journal.pone.0159864.g006

Weak and non-significant genetic structure was detected between Malta and Linosa (F_ST = 0.0017, p > 0.05), and Malta and N. Sicily (F_ST = 0.0084, p > 0.05) whilst the most significantly divergent population was found to be from Libya (F_ST = 0.0345, p < 0.05). All Mediterranean populations were found to be genetically differentiated from Hurd Bank (Table 4). Relative ancestral population size estimates were not significantly different between populations. Divergence times between populations generally corresponded to the level of coancestry where metapopulations within the Sicily Channel exhibited comparatively high pairwise coancestry coefficients along with relatively recent divergence times (Fig 7). Reynolds’ F_ST (coancestry coefficient) was consistent with results of the unrooted NJ phenogram based on δμ² depicting the genetic relationship amongst regional progenitor dusky grouper populations which also displayed a network of connectivity within the Sicily Channel (Fig 8).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Pairwise genetic population division (F_ST) and differentiation (Jost’s D) between seven Mediterranean Epinephelus marginatus populations.

https://doi.org/10.1371/journal.pone.0159864.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Demographic history and genetic distance within and between Mediterranean Epinephelus marginatus populations evaluated with Arlequin.

https://doi.org/10.1371/journal.pone.0159864.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Phylogenetic tree of Epinephelus marginatus within the central Mediterranean region constructed with the software Poptree2.

Shown is a neighbour-joining tree based on δμ² microsatellite genetic distance calculated between each sampling location. Results reveal the biogeographically isolated populations of Linosa, Malta and Hurd Bank (highlighted in red) on the same branch suggesting a network of connectivity located within the Sicily Channel. Bootstrap values shown at nodes.

https://doi.org/10.1371/journal.pone.0159864.g008

Findings of the Garza-Williamson modified index for the detection of historical population size reduction found that all populations have experienced a past decline in size and that more than half (Hurd Bank, Croatia, Libya, N. Sicily) were symptomatic of remnant populations (Fig 9). Hurd Bank displayed a significantly low level of allelic richness when compared to neighbouring populations and analysis with Bottleneck detected clear evidence of a recent bottleneck when measured by significant heterozygous excess (p < 0.05) and mode shift of allelic distribution at mutation-drift equilibrium (Table 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Detection of Epinephelus marginatus population size reduction with the programme Arlequin.

https://doi.org/10.1371/journal.pone.0159864.g009

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Evidence of bottleneck and long-term potential of adaptability in Mediterranean Epinephelus marginatus populations.

https://doi.org/10.1371/journal.pone.0159864.t005

Investigation into genetic barriers associated with geographic location using Barrier identified a single significant barrier to gene flow between Hurd Bank and surrounding populations (Fig 10). Interestingly marine biogeographic barriers of land bridges (e.g. Italian Peninsula, Sicily Island) and hydrological processes (e.g. surface current confluence and recirculation patterns) were not identified as barriers to pelagic larval dispersal using this method. The absence of any obvious biogeographic barriers between Hurd Bank and neighbouring populations indicates a biological or ecological barrier to gene flow may be present. Results of the Mantel test showed that genetic distance is not significantly correlated with geographic distance (Pr = 0.423, p > 0.05), suggesting that isolation by distance is not the principal factor influencing population substructuring and connectivity in Mediterranean dusky groupers (S5 File).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Geographic barriers of Epinephelus marginatus based on θ_ST by the programme Barrier.

Hurd Bank, MT (HB), Malta, MT (MAL), Linosa, IT (LIN), Croatia (CRO), Libya (LIB¹, LIB², LIB³), N. Sicily, IT (SCL) and Tunisia (TUN). Barrier significance (red) based on 1,000 θ_ST bootstraps matrices.

https://doi.org/10.1371/journal.pone.0159864.g010

Spatial Bayesian analysis with Geneland used to assess connectivity to Malta’s FMZ produced a pattern of genetic clustering similar to previous analyses where individuals from the islands of Linosa and the Maltese archipelago (Malta, Gozo, Comino), were found to be a homogeneous genetic cluster. Within the FMZ, individuals from Hurd Bank, an area 12 nm east of the main island, were again identified as a unique cluster increasing the evolutionary interest that the shallow water corridor between Malta and Sicily may play as it relates to population substructuring. The remaining neighbouring populations were generally assigned to either the Western Mediterranean Sea or North African/Eastern Mediterranean Sea with the notable exception of individuals clustered near an area of hydrological recirculation near Tripoli (Fig 11).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Spatially explicit estimates of Epinephelus marginatus to genetic clusters within the Mediterranean Sea using Geneland.

Posterior density distribution estimates with Geneland identified four (K = 4) unique genetic clusters within the central Mediterranean region. Each map represents the posterior probability of an individual belonging to one of these unique spatial clusters where black dots depict E. marginatus individuals. Probability of assignment to cluster is visualised by contour lines where the highest probability is shown in white and the lowest in red. (a) N. Sicily, IT (SCL), Libya (LIB¹); (b) Croatia (CRO), Tunisia (TUN), Libya (LIB², LIB³); (c) Hurd Bank, MT (HB); (d) Malta, MT (MAL), Linosa, IT (LIN).

https://doi.org/10.1371/journal.pone.0159864.g011

Direct estimates of dusky grouper migrants with Geneclass2 detected exchange of F₀ individuals between Malta and Linosa, located within the Sicily Channel. One potential immigrant from Tunisia was identified in the Hurd Bank subpopulation and no exchange of F₀ Hurd Bank individuals were detected between any other demes (S6 File). Evolutionary and contemporary directional gene flow per generation analysed with Migrate-n and diveRsity highlight the Sicilio-Tunisian Channel as a biogeographical crossroad between the Western and Eastern Mediterranean Sea, acting as a regional hub, in facilitating exchange of dusky grouper migrants (Table 6, Fig 12). Hurd Bank appears to have a limited pattern of migrant exchange following the recent bottleneck event.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Ancestral effective number of Epinephelus marginatus migrants per generation (Nm) in the central Mediterranean region estimated with the programme Migrate-n.

https://doi.org/10.1371/journal.pone.0159864.t006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Visualisation of contemporary gene flow between seven Mediterranean Epinephelus marginatus demes with diveRsity.

Based on projected direction of gene flow and magnitude (G_ST), the Sicilio-Tunisian Channel appears to act as a central hub of exchange between the Western and Eastern Mediterranean Sea. Hurd Bank, MT (HB), Malta, MT (MAL), Linosa, IT (LIN), Croatia (CRO), Libya (LIB), N. Sicily, IT (SCL) and Tunisia (TUN).

https://doi.org/10.1371/journal.pone.0159864.g012

Malta FMZ: Demography

Malta’s dusky grouper population decline appears to be part of a greater trend within the Mediterranean Sea where catch landings reported from seven Mediterranean countries between 1990 and 2001 showed a collective decline of 88% over this time due to overexploitation [3]. Evidence documented during this study also revealed a moderate mean catch length (52 ± 17 TLcm) where over one-third (38%) of specimens measured below the set minimum catch size limit of 45 TLcm for grouper spp. [116]. The mean catch length range of commercially caught dusky groupers in Malta raises concern of population extinction probabilities, where Individual Based Model (IBM) simulation tests found protection of dusky grouper females at a flexible life stage (52–77 TLcm, 7–17 years [16]) were the most important demographic parameter to protect in this protogynous hermaphroditic species to achieve population stability [24]. Whilst the minimum catch size limit of 45 TLcm was not effectively enforced during this study, had it been, it hypothetically may have caused a demographic shift toward targeting of the female/transitional-sex individuals where the expected outcome would have further destabilised the local population and delayed recovery. Based on this information, a minimum catch size limit as a strategy for conservation of the Maltese dusky grouper is not recommended.

The mean weight of Maltese dusky groupers sampled during this study was 3 ± 3 kg, which is considerably smaller than the weight (5 kg) at first maturity previously reported for the Maltese [22,23] and neighbouring Tunisian dusky groupers [18]. However, Ben Miled et al. [24] recognised a direct correlation between large populations of dusky groupers at carrying capacity (37.3 g m^-2 [117]) and an increased growth rate and age at sexual maturation and inversion. Therefore, it can be expected that a decreasing trend in the overall population size would result in individuals within that population to grow smaller and mature at an earlier age as a coping mechanism.

Sex allocation estimates of the Maltese dusky grouper population between 2007 and 2009 (17% juveniles, 68% female/transitional-sex and 15% male) were found to be comparable with a proximal population of dusky groupers from the Pelagie Islands, IT between 1994 and 1997 before it was established as an MPA (12% juveniles, 69% female/transitional-sex and 19% male [21]). The male-to-female sex ratio estimate in Malta of 1:5 is also consistent with six previous studies throughout their native global range which varied between 1:1.1 and 1:7.4 [16,19,21,29,118,119].

A comparison of the weight-length relationships of Maltese dusky groupers to discrete global populations (listed in descending order of size (g/TLcm): Senegal, W. Algeria, Brazil¹, Egypt, Azores, French Lavezzi Islands, S. Italy, Malta [this study], Brazil², E. Algeria) found that the Maltese dusky grouper mean body size was modest; however it was in accordance with the body size of individuals from the neighbouring populations of S. Italy and E. Algeria in the central Mediterranean region suggesting that this may be a regional trait [13,17,20,120–125]. As discussed by Reñones et al. [16], the number of eggs laid is a function of body size where potential fecundity increases exponentially with size, indicating a decreased ability of the dusky grouper to repopulate in this region.

Malta FMZ: Population genetic structure

All 14 microsatellite loci were polymorphic and thus useful for describing the population substructure within Malta’s FMZ. Observed and expected heterozygosity for microsatellites in the Maltese dusky grouper subpopulations (H_O = 0.691–0.795; H_E = 0.695–0.759) were comparable to those in a similar Mediterranean study (H_O = 0.658–0.768; H_E = 0.703–0.762) using the same technique [35]. A moderate and significant inbreeding coefficient was observed within the FMZ (F_IS = 0.10, P < 0.001) when compared to these nine neighbouring Mediterranean populations (F_IS = -0.058–0.065 [35]) where five of these populations were outbreeding and none showed LD within any population. The high level of genotypic disequilibrium observed within Hurd Bank is of conservation concern and symptomatic of a destabilised population and consequently reduced N_e.

Within the management zone, measurement of adult effective population size was low, estimated to be around 130 individuals with fewer than a dozen effective breeders per reproductive cycle. In order for a species to retain its evolutionary potential, to avoid localised extinction and survive deterministic and stochastic events, it is recommended that immediate conservation action be put in place in Malta to restore natural dynamics. As a general guideline, in order to ensure short-term survival and minimise inbreeding depression within populations, Franklin [126] proposed that the Minimum Viable Population (MVP) should not be allowed to fall below 50 (N_e) individuals. However, as Jamieson and Allendorf [127] point out, this guideline is based on an averaged N_e/N_census ratio of 0.10 [128], where high-fecundity groups such as marine fishes tend to have lower ratios. Whilst the N_e/N_c ratio for dusky groupers is currently unresolved, for the gag grouper (Mycteroperca microlepis) which has similar biology and generation time it is estimated to be on the order of 0.01 [129], where a shift in female bias sex ratios is expected to have detrimental long-term effects that would not immediately become detectable [130]. Therefore, it is important that monitoring of sex ratios be included in the design of any FMZ conservation management plan.

Approximately 90% of dusky groupers sampled within the FMZ were identified as resident individuals with evidence of localised genetic clustering found between the island complex and Hurd Bank, located 12 nm offshore. Genetic homogenisation and consistent gene flow within the archipelago, based on observed surface current direction during spawning season (S2 File), could be maintained by larval migration originating from the northern tip of Gozo, acting as a source population within the archipelago. Rare breeding size males (≥ 85 TLcm) were present throughout the islands with a higher incidence of large size individuals (≥ 60 TLcm) observed in northern Gozo near Dwejra and between Marsalforn-Ramla Bay. Northern Gozo is an area that consists of a complex topography composed of coralligenous bioherm, small caves, tunnels, large and medium sized boulders that offer numerous cryptic shelters as well as nearby tidal pools with dense algal coverage and sea grass beds. It is an ideal environment for the Maltese dusky grouper at all stages of the life cycle. Therefore, it is suggested that conservation management include monitoring of the prospective source populations in Dwejra and Marsalforn-Ramla Bay, Gozo to better understand local dynamics.

Regional connectivity

Regional patterns of genetic structure observed were similar to results seen by Schunter et al. [35] and de Innocentiis et al. [34] that found both areas of genetic partitioning and intra-Mediterranean panmixia, but not panmixia throughout the Mediterranean Sea. Neither this study, nor Schunter et al. [35] or de Innocentiis et al. [34], detected isolation by distance as the principal factor influencing genetic structure. Spatially explicit genetic clustering analysis and temporal patterns of migration adds further detail showing an interconnected network within the Sicilio-Tunisian Channel between the Western and Eastern Mediterranean Sea basins.

Generally, in order to minimise loss of genetic diversity, between one and 10 migrants per generation is considered effective at maintaining within population heterozygosity and preventing long-term deleterious loss of polymorphisms [131]. For newly founded or severely restricted populations, such as Hurd Bank, this level of migration may be inadequate due to the loss of allelic variation where recovery of allelic richness, as opposed to frequency, would be a more meaningful demographic measurement [132]. For comparison, it was estimated 20–30 migrants per generation would be adequate to maintain homogeneity between populations of the gag grouper (M. microlepis) with an estimated regional N_e of 16,500 [130]. Although regional genetic differentiation of the Mediterranean dusky grouper was found to be weak, inadequate levels of migration leading to inbreeding depression would not likely become detectable for several hundred generations.

Influences of geographic and genetic isolation of Hurd Bank leading to significant genetic partitioning and differentiation within the region may be attributed to the founder effect. The founder effect is a type of bottleneck caused by a selective sweep or colonisation event by a few individuals. Whilst several founder models have been proposed [133–136], Hurd Bank exhibits signature signs of a population in expansion following a severe reduction in N_e [130]. Small effective populations (N_e) initially display a relative reduction in genetic and allelic variation and an increase in LD caused by genetic drift. This is often followed by population expansion and an increase in allelic frequencies, but not allelic richness (due to loss of rare alleles) which leads to divergence from progenitor populations in absence of selection and immigration. Based on life history traits and environmental needs of dusky grouper settlers within shallow seagrass beds, the existence of an open water population at 40–50 m depth is somewhat surprising. Due to genetic divergence from nearby demes, juvenile stage migration related to intraspecific interaction, ontogenetic changes or seasonal micro-migration does not appear to be a plausible explanation [14]. Investigation into translocated individuals from a 2005 pilot E. marginatus restocking programme [137,138] that released 95 F1 juveniles (3–4 years old) from SW Sicily, located 120 nm NW of Hurd Bank, roughly 1.5 years prior to sampling as an explanation for genetic differentiation of the Hurd Bank population was also ruled out due to dissimilar genetic profiles (S7 File). Hurd Bank, formed of coral and sand, is an area that is heavily trawled and currently serves as an anchorage for large ships, leading to concerns of coral destruction and ecosystem imbalance [139]. Therefore, further genetic sampling coupled with an ecological survey of Hurd Bank is recommended.

Conservation strategy

Effects of geographic productivity and environment on dusky grouper size and population density have been documented from several protected areas throughout the Western Mediterranean Sea, including Ustica Island Marine Reserve, IT established in 1986 [117,140–151]. These studies indicate an established positive correlation between the degree of habitat protection and increase in biomass linked to growth rate as well as a causative homing range shift and population restructuring [150,152]. A ban on spearfishing appears to have a similar effect on biomass increase affecting density and growth rates in localised populations [144]. Unprotected spawning aggregates of the migrating (220 km [153]) Nassau grouper (Epinephelus striatus) that have been heavily targeted by fishermen are not known to re-establish once they have been fished out [154–156]. Although dusky groupers only display minor migrations to spawning sites, targeting of trophy males within spawning aggregations around Malta is of conservation concern [157], and therefore recommendations for future studies include identifying local spawning event sites so that they may be considered in the design of MPA no-take/no-entry zones. It is estimated that in well managed MPAs, a resurgence of the dusky grouper population can be seen in as little as three years following protection, and carrying capacity (95% K) can be achieved following 21 to 29 years of protection [117,141]. Therefore, the combination of local habitat degradation along with the absence of an effective conservation plan aimed at increasing biomass, may explain the modest size of Maltese dusky groupers.