Genetic population structure of the Blackspot seabream (Pagellus bogaraveo): contribution of mtDNA control region to fisheries management

Abstract Marine fisheries management models have traditionally considered biological parameters and geopolitical boundaries. The result is the existence of fisheries management units that do not match genetic populations. However, this panorama is changing with the contribution of genetic and genomic data. Pagellus bogaraveo is a commercially important sparid in the northeast Atlantic, with three stock components being considered by ICES: the Celtic Sea and Bay of Biscay, Atlantic Iberian waters and the Azores. The northern stock collapsed (1975–1985) and is essential to characterize the genetic makeup of the species, particularly in the Iberian Peninsula, where it is managed as a single stock. The mitochondrial control region was used to screen the intraspecific diversity and population structure of individuals from six locations across the species range. The genetic diversity found is similar among sites, and there is differentiation between the Azores and the remaining locations.


Introduction
The assessment of the exploitation status of marine resources requires information on various aspects of their populations, such as structure, spatial and temporal dynamics, as well as aspects of their biology, namely growth and reproduction. Marine fisheries management models have traditionally taken into account biological parameters such as weight and length, larvae abundance, fishing mortality, recruitment, spawning stock biomass and landings, through a scientific process that also involved political considerations, such as boundaries. The result is the existence of fisheries management units that do not match the distribution of biological (in the genetic sense) natural populations (Kerr et al. 2016). However, the panorama is changing with the ever-increasing genetic and genomic data contributing to the identification of natural populations (Benestan et al. 2016;Kerr et al. 2016;Benestan 2019;Hohenlohe et al. 2020) and the detection of connectivity patterns in species showing low levels of population structuring Zang et al. 2020).
Despite the increasing number of studies of genomics in stock determination (Pecoraro et al. 2018;Schulze et al. 2020) mitochondrial DNA keeps being used as a genetic tool mainly due to its molecular characteristics. MtDNA is a small (around 16,000 bp in fishes) maternally inherited circular DNA molecule contained in large quantities inside the mitochondria, and it is relatively stable in terms of gene content. Mitochondrial DNA accumulates differences faster due to a high mutation rate and smaller effective population size (Ne) comparatively to nuclear DNA, making mitochondrial markers excellent for differentiating divergent species or populations. Nevertheless, within the mtDNA molecule, genes evolve at different rates, with slowly evolving genes (like cytochrome oxidase I, COI) being used for comparisons at the species (e.g. barcoding, allowing for species identification, either in whole species or fish products) and intraspecific level, proving its value as a stock assessment tool (Ovenden 1990;Antoniou and Magoulas 2014).
The blackspot seabream (Pagellus bogaraveo Br€ unnich 1768) is a common benthopelagic sparid in the northeast Atlantic, ranging between 150 and 700m deep, from Norway to Mauritania, Madeira, Canaries, and the western Mediterranean (Froese and Pauly 2019). It is a high-value commercially important species across its distribution area, with three stock components being considered by the International Council for the Exploitation of the Sea (ICES): a) Subareas 6, 7, and 8 (Celtic Seas and the Bay of Biscay); b) Subarea 9 (Atlantic Iberian waters); and c) Subarea 10 (Azores region) (ICES 1996). In the northern component (ICES Subareas 6-8), the stock collapsed during the period 1975-1985 in the Bay of Biscay (Lorance 2011). With the establishment of Total Allowable Catches (TAC) in 2003, the fishery was restrained, resulting in a steep decline in catches (currently 1% of the historical TAC level) but the recovery of the stock will be slow due to the slow growth and late maturity of the species (ICES 2020). This situation led ICES to repeatedly recommend, since 2016, a precautionary approach with zero catch advice along with the establishment of a recovery plan (ICES 2016(ICES , 2018(ICES , 2020. Presently, no directed fisheries are permitted in that area. Hence all catches are bycatches according to EU Regulations (ICES 2020).
In the southern component of the fishery, landings originate from Portuguese vessels that only target this species seasonally in the mainland coasts, and from the Spanish target fishery that operates in the Strait of Gibraltar (the 'voracera'), corresponding to 35% and 65% landings in ICES Subarea 9, respectively (ICES 2020). The 'voracera' fishing grounds encompass areas under management or advice from different Regional Organizations/Commissions, namely, ICES, the General Fisheries Commission for the Mediterranean (GFCM), and the Fishery Committee for the Eastern Central Atlantic (CECAF). Still, the EU TAC only applies to ICES Subarea 9. However, there is a clear declining trend in landings and in Catch Per Unit Effort (CPUE) levels from the target fishery in the Strait of Gibraltar, which might be consistent with an overexploited population, a different observed trend in CPUE for a reference fleet operating in mainland Portugal (ICES 2020).
In the Azores, the species is targeted by a demersal multispecies and multi-gear fishery (ICES 2020). Here, landings were at their highest between the 1990s and 2010 (fluctuating around 1000 t) and decreased significantly onwards (600 t) due to the introduction of the TAC. Currently, the fishery is still highly constrained by management measures.
Before establishing a recovery plan for P. bogaraveo, it is essential to adequately characterize it particularly addressing the putative population structure where the species is globally considered and managed as a single stock (ICES Subarea 9: south of the Iberian Peninsula, adjacent to the Mediterranean and North Africa). The differences reported between the Portuguese CPUE and that observed for the Strait of Gibraltar further supports the need to reassess the population structure of the species within the area (ICES 2020).
Since the beginning of this century, several cruises have marked individuals in the waters of the Strait of Gibraltar, which have revealed an absence of significant migrations for this species (Gil et al. 2018;ICES 2018). In parallel investigations, Portuguese researchers have observed evidence of heterogeneity in the distribution of this species, which appears to prefer aggregating in specific locations (Farias et al. 2018). These authors reported an uneven distribution in-depth and latitude, with individuals sampled in southern Portugal preferring habitats between 200 and 400m and those in the north up to 100m (Farias et al. 2018). This difference in depth, associated with gregarious characteristics and low migration, may promote the existence of population differentiation in areas without apparent barriers to species dispersion.
The characterization of genetic diversity of blackspot seabream and the corresponding population structure has been the subject of a small number of studies using few sampling locations to compare larger areas, i.e. the Atlantic and the Mediterranean. Published studies involved mitochondrial DNA and microsatellites, and in general did not reveal population structure among Atlantic and Mediterranean locations, except for the Azores (Bargelloni et al. 2003;Stockley et al. 2005;Piñera et al. 2007). Nevertheless, the division of the species distribution area into three components (Azores, Iberian Peninsula, and the Mediterranean) has been supported by morphometric measurements (Palma and Andrade 2004) and parasitology (Hermida et al. 2013).
Considering that the studies developed to this date did not include sampling locations considered representative of the distribution of the species, this is the most geographically comprehensive genetic study on the blackspot sea bream done so far and it aims to unravel the population structure of the species, contributing to its management.

Sampling
Collection of P. bogaraveo specimens (N ¼ 127) covered six locations across the species distributional range in the Northeast Atlantic and the Mediterranean (Figure 1(A) and Table 1): Azores (Azo -Faial; Portugal), France (FRA -Ile-de-Sein), Cantabria (CAN -San Vicente de la Barquera, Luarca; Spain), Peniche (PEN -Portugal), Cadiz (CAD; Spain) and Mediterranean (MED -Malaga, Almeria; Spain). Fin clips were preserved in 96% ethanol until DNA extraction. The authors contacted colleagues and fishermen to find specimens in the extreme north and south of the species distribution area, but the species is absent or very rare. No additional sequences were retrieved from Genbank because only haplotypes with no associated frequencies were available (and sometimes only the most common is deposited) and/or the sequences were not geo-referenced.

DNA extraction, amplification, and sequencing
Total genomic DNA was extracted with the REDExtract-N-Amp Kit (Sigma-Aldrich) following the manufacturer's instructions. The mitochondrial control region (CR) was amplified, in a Bio-Rad Mycycler thermal cycler, using L-pro1 (5 0 -ACTCTCACCCCTAGCTCCCAAAG-3 0 ) and H-DL1 (5 0 -CCTGAAGTAGGAACCAGATGCCAG-3 0 ) (Ostellari et al. 1996). The PCR protocol was performed in a 20 ll total reaction volume with 10 ll of REDExtract-N-ampl PCR mix (Sigma-Aldrich), 0.8 ll of each primer (10 lM), 4.4 ll of Sigma water and 4 ll of template DNA, using the following PCR conditions: initial denaturation at 94 C for 7 0 , followed by 35/30 cycles (denaturation at 94 C for 30 00 , annealing at 55 C for 30 00 , and extension at 72 C for 1 0 ) and a final extension at 72 C for 7 0 . The same primers were used for the sequencing reactions. PCR products were purified and sequenced in STABVIDA (http://www.stabvida.net/). Chromatograms were edited with Codon Code Aligner (Codon Code Corporation, http://www.codoncode.com/index. htm), and sequences were aligned with Clustal X 2.1 (Larkin et al. 2007). All sequences were deposited in GenBank (Accession numbers MW251181 -MW251307).

Genetic analyses
PopART software (Leigh and Bryant 2015) was used to build a TCS haplotype network (Clement et al. 2000) based on the statistical parsimony algorithm of Templeton et al. (1992). The ARLEQUIN software package v3.5.1 (Excoffier and Lischer 2010) was used to estimate standard descriptive measures of genetic diversity, including the number of haplotypes and private haplotypes, haplotype diversity h (Nei 1987) and nucleotide diversity p (Nei 1987). To assess population structure, the same software was used to perform analyses of molecular variance (AMOVA) (Excoffier et al. 1992) and compute pairwise F ST estimates. The correlation between geographic distance and F ST was calculated with the Mantel test (Mantel 1967;Smouse et al. 1986), using 10,000 permutations and geographic distances estimated with Geographic Distance Matrix Generator v.1.2.3 (Ersts 2020).

Results
The final alignment contained 428 bp in length. The 127 individuals analysed were grouped into 38 distinct haplotypes, 12 (32%) of which were shared between sampling sites. The resulting network reveals a central high-frequency haplotype, which is shared by individuals from all locations (Figure 1(B)). The haplotype network has a star shape and does not show any association with the geographic origin of the samples.
The genetic diversity found is similar for the six sampling sites (Table 1). The percentage of haplotypes shared with other locations was higher for Cantabria and Peniche, with the latter being distinguished by another parameter at the level of mutations, the number of insertions/deletions (indels). The Mediterranean and Azores were the locations with the lowest percentage of shared haplotypes, with the Azores still presenting a smaller number of mutations (only six transitions).
Concerning the analysis between populations, a significant population structure was found for the entire sample area (AMOVA: F ST ¼ 0.052, p <0.001). The genetic structure found is attributed to the individuals from the Azores, which showed significant differences with the remaining five sites ( Table 2). The Mantel test did not support a significant correlation between the geographical distance and F ST (r ¼ 0.702, p ¼ 0.112), thus ruling out the hypothesis of isolation-by-distance for the blackspot seabream.

Discussion
The lack of high-resolution genetic data for P. bogaraveo and other commercial species is not the only obstacle to the estimation of the population structure of the species. The fact that different studies have used different markers also makes it impossible to compare published sequences directly. Another challenge that arises is that some studies do not make available online all the sequences (or do not reveal the geographic origin of the sample) so that the meta-analysis of published data becomes impossible (Stockley et al. 2005). Therefore, although providing only mitochondrial data, this is the most geographically comprehensive genetic study on the blackspot sea bream.  Table 1. (B) Haplotype network for the mitochondrial control region of Pagellus bogaraveo. Colours refer to sampling locations. The area of the circles is proportional to each haplotype frequency. In the case where haplotypes are shared among sampling locations, shading is proportional to the frequency of the haplotype in each sampling location. The genetic diversity of Pagellus bogaraveo displayed in the six locations is comparable, and there is genetic differentiation between the Azores archipelago and the remaining Iberian locations. Before discussing the results, we address the study's main caveats: the use of a single mitochondrial marker with a fast mutation rate and/or a relatively small number of specimens sampled. Concerning the first, the mitochondrial control region displays a high mutation rate and a somewhat irrelevant genetic recombination, which promotes regional variation within species, making it a highly used marker for a preliminary approach to stock delimitation (Gao et al. 2019;Song et al. 2020). The fact that CR corroborates the results previously obtained even using nuclear markers (i.e. microsatellites) adds strength to the validity of its use (Piñera et al. 2007). Concerning the second caveat, the number of individuals sampled follows similar studies for the same geographical area, allowing for comparisons among relevant studies, using the same marker (Francisco et al. 2014;Almada et al. 2017;Robalo et al. 2020).
The blackspot seabream displays a marked genetic structure between the Azores archipelago and the remaining Atlantic continental locations, corroborating the results obtained by Stockley et al. (2005). The differentiation of the Azorean population is common to other species of the Northeast Atlantic, namely the ballan wrasse Labrus bergylta (Almada et al. 2017), the black scabbardfish Aphanopus carbo (Stefanni and Knutsen 2007), the blenny Coryphoblennius galerita (Francisco et al. 2014) and the two-banded seabream Diplodus vulgaris (Stefanni et al. 2015). The haplotype network does not reveal any geographical pattern or association, with a very abundant central haplotype shared by individuals from all sampling locations. A similar haplotype network shape was previously found with samples from the Azores and three Mediterranean locations (Bargelloni et al. 2003). Using microsatellites, Piñera et al. (2007) found genetic homogeneity between blackspot seabream from the northern and eastern coasts of Spain, i.e. between the Spanish Atlantic and Mediterranean coasts. Here, we have extended the scope of the previous work by including samples from the French coast and the western and southern coasts of the Iberian Peninsula (Atlantic), closer to the entrance to the Mediterranean. Results across two different marker types, mitochondrial and nuclear, were consistent: no significant differences between the four Atlantic locations and the Mediterranean locations were found.
Despite the conclusions of this study, our results also highlight the need for a genomic survey. In this context, the massive sequencing of thousands of genetic markers simultaneously would be an essential tool for clarifying the population structure of this species, namely the use of single nucleotide polymorphisms (SNP) obtained by DNA sequencing associated with the restriction site (RADseq) (Andrews et al. 2016, Benestan 2019. Recent advances in the field of genomics have allowed for a more refined and powerful resolution of the spatio-temporal structure of populations of diverse fishing resources (Helyar et al. 2011, Larson et al. 2014, Pecoraro et al. 2018). Just as an example, recent studies on yellowfin tuna (Thunnus albacares; Pecoraro et al. 2018) using SNPs obtained by RADseq managed, for the first time, to identify discrete populations of this species in the Atlantic Ocean, stressing the need to integrate this type of information in stock management.