Genetic investigation of population structure in Atlantic chub mackerel, Scomber colias Gmelin, 1789 along the West African coast

Sampling

Fish were collected during the scientific surveys on board the R/V Dr. Fridtjof Nansen carried out in 2017 and 2019 in the area stretching from Morocco (~28°N) to Namibia (~22°S) (Fig. 1). Sampling permits were obtained through the “Application of the Ecosystem Approach to Fisheries management considering climate and pollution impacts” (GCP/GLO/690/NOR). Pelagic trawl hauls were performed after echosounders target identification of fish schools and, when possible, 30 individuals were collected per trawling station (Table 1). The sampling procedure reflected the species’ distribution (Baird, 1977; Scoles, Graves & Collette, 1988), whose southern limit of highest abundance seems to be placed at low latitudes in the southern hemisphere. Thus, some 75.5% of the individuals were collected north to Ghana whereas the southernmost countries such as Angola and Namibia only accounted for 7.5% of the individuals in line with the species biomass levels observed in those zones. Fin clips from 1,848 fish were collected and preserved in ethanol prior to DNA extraction.

DNA isolation and genotyping

DNA was extracted using the Qiagen DNeasy 96 Blood & Tissue Kit in (Qiagen, Hilden, Germany) 96-well plates; each of which contained two or more negative controls. DNA concentration was quantified using a NanoDrop 8000 (Thermo Fisher Scientific, Waltham, MA, USA). The molecular markers used in the current study were initially developed for both S. colias (one) and S. japonicus (seven) (Catanese et al., 2010; Chen et al., 2017; Zeng & Cheng, 2012) as microsatellite cross-species amplification has proven a successful tool for fisheries management and conservation (e.g., Maduna et al., 2014) formerly applied to Scomber species (Tang et al., 2009; Zeng, Cheng & Chen, 2012). Eight microsatellite loci isolated from S. colias (Sco2_1) and S. japonicus (SJ78, SJNT19, SJT5, SJT53, SJT122, SJT182, SJT199), respectively (Table S1 in Supporting Information), were genotyped in three multiplexed reactions. PCR amplifications were performed in a final volume of 10 μL containing approximately 20 ng DNA template, 1 × buffer, 2–3 mM MgCl₂, 0.2 mM dNTPs, 0.2–0.5 μM of each primer and 1U GoTaq polymerase in an Applied Biosystems Gene Amp PCR Systems 2700 thermal cycler. PCR profiles consisted of 15 min denaturation at 95 °C followed by 27–32 cycles of 30 s denaturation at 94 °C, 90 s at annealing temperature (Table S1 in Supporting Information) and 90 s extension at 72 °C with a final step of 7 min at 72 °C. PCR products were analysed on an ABI 3730 Genetic Analyser and the 500LIZ^TM size standard (Applied Biosystems, Waltham, MA, USA) was used to accurately determine the size of the fragments and allelic variation. Scoring was performed using the software Fragment Profiler 1.2 (Amersham Biosciences, Buckinghamshire, England). Automatically binned alleles were manually checked by two researchers prior to exporting data for statistical analyses.

As it was not always possible to obtain a minimum of thirty good-quality genotyped individuals per site, nearby sampling sites within countries were merged to achieve significant sampling sizes (Hale, Burg & Steeves, 2012). Sampling sizes ranged from 22 individuals (GAM in Gambia) to 58 (MOR_2 in Morocco). However, in Sierra Leone (N = 11), and Guinea-Bissau (N = 5) (Table 1), sampling sizes were too low and therefore these samples were excluded from some of the analyses.

Statistical analysis

Loci were screened for null alleles, large allele dropouts and potential scoring errors due to stuttering bands using the software Micro-Checker v.2.2.3 (Van Oosterhout et al., 2004). The frequency of null allele(s) was estimated with the maximum likelihood method using the EM algorithm of Dempster, Laird & Rubin (1977) implemented in the software Genepop 7 (Rousset, 2008). To assess the effect of including loci with possible null allele(s) on population differentiation estimates, the software FreeNA (Chapuis & Estoup, 2006) was used to calculated both uncorrected and corrected F_ST values. Confidence intervals (95%) of null frequencies were based on 1,000 bootstraps.

The statistical power of this microsatellites set to detect genetic differentiation was estimated using POWSIM 4.1 (Ryman & Palm, 2006), which uses χ² and Fisher tests to assess whether the observed data set carries enough statistical power to detect a Nei’s F_ST value significantly larger than zero. The analyses were conducted maintaining effective population size constant (N_e = 15,000) following Henriques et al. (2017), and changing the number of generations of drift since divergence (t: 15, 35, 100, 160, 200, 350 and 500 generations) to model different levels of F_ST (0.0005, 0.0012, 0.0033, 0.0053, 0.0066, 0.0117, and 0.0165, respectively). The probability of rejecting the null hypotheses of no population differentiation (type I error) was also estimated (F_ST = 0, t = 0). Statistical power was determined as the proportion of tests indicating statistical significance (p < 0.05), with 1,000 replicates. In addition, to assess if this suite of twelve microsatellites would accurately discriminate between individuals in a population, the genotype accumulation curve was built using the function genotype curve in the R (R Core Team, 2020) package poppr (Kamvar, Tabima & Grünwald, 2014) by randomly sampling x loci without replacement and counting the number of observed multilocus genotypes (MLGs).

The total number of alleles, number of private alleles, and allelic richness per sample was calculated using MSA 4.05 (Dieringer & Schlötterer, 2003). The observed (H_o) and unbiased expected heterozygosity (uH_E), inbreeding coefficient (F_IS) as well as the number of deviations from Hardy-Weinberg expectations (HWE) were computed per sample using GenAlEx v6.1 (Peakall & Smouse, 2006). Linkage disequilibrium (LD) between pairwise loci per sample was computed using the program Genepop (Rousset, 2008). The false discovery rate (FDR) correction of Benjamini & Hochberg (1995) was applied to p-values to control for Type I errors.

In marine fish species, loci carrying signatures of locally divergent selection might function as powerful markers to evaluate spatially explicit genetic structure, as well as to outline fisheries stocks for sustainable management (Russello et al., 2012). Thus, two approaches were combined to identify loci eventually departing from neutrality: BayeScan 2.1 (Foll & Gaggiotti, 2008) and Arlequin v.3.5.1.2 (Excoffier, Laval & Schneider, 2005). In BayeScan, sample size was set to 10,000 and thinning interval to 50; and loci with a posterior probability over 0.99, corresponding to a Bayes Factor >2 (i.e., “decisive selection” (Foll & Gaggiotti, 2006), were retained as outliers. In Arlequin, analysis was simulated based on 1,000 demes with 50,000 simulations under a hierarchal island model.

Pairwise F_ST (Weir & Cockerham, 1984) was computed using GENEPOP 4.0.6 (Rousset, 2008), whereas population differentiation was tested via G exact tests calculated using the following Monte Carlo Markov Chain (MCMC) parameters: 10,000 steps of dememorisation, and 5,000 iterations for 100 batches, also with GENEPOP 4.0.6.

The Bayesian clustering approach implemented in STRUCTURE v.2.3.4 (Pritchard, Stephens & Donnelly, 2000), and conducted using the software ParallelStructure (Besnier & Glover, 2013), was used to identify genetic groups under a model assuming admixture and correlated allele frequencies, both using and without using geographic LOCPRIORS to assist the clustering. Ten runs with a burn-in period consisting of 100,000 replications and a run length of 1,000,000 MCMC iterations were performed for K = 1 to K = 10 clusters. To determine the number of genetic groups, STRUCTURE output was analysed using two approaches: (a) the ad hoc summary statistic ΔK of Evanno, Regnaut & Goudet (2005) and (b) the Puechmaille (2016) four statistics (MedMedK, MedMeanK, MaxMedK and MaxMeanK), specially recommended for uneven sampling sizes, both conducted using StructureSelector (Li & Liu, 2018). Finally, the ten runs for the selected Ks were averaged with CLUMPP v.1.1.1 (Jakobsson & Rosenberg, 2007) using the FullSearch algorithm and the G′ pairwise matrix similarity statistic, and graphically displayed using barplots. Furthermore, the relationship among samples was also examined using the Discriminant Analysis of Principal Components (DAPC) (Jombart, Devillard & Balloux, 2010) implemented in the R (R Core Team, 2020) package adegenet (Jombart, 2008) in which groups were defined a priori using geographically explicit samples. To avoid overfitting, both the optimal number of principal components and discriminant functions to be retained were determined through cross validation using the xvalDapc function from adegenet (Jombart & Collins, 2015; Miller, Cullingham & Peery, 2020).

The relationship between genetic (F_ST) and geographic (km) distance was examined to investigate if it followed the expectations of an “Isolation by Distance” pattern (IBD), i.e., increasing genetic differentiation with geographic distance as a result of restricted gene flow and drift (Rousset, 1997; Slatkin, 1993; Wright, 1943). A two-tailed Mantel (1967) test was conducted using PASSaGE v2 (Rosenberg & Anderson, 2011) and significance was assessed via 10,000 permutations. The matrix of pairwise shortest distance by water was created by calculating least-cost distances via seas (avoiding landmasses) between sampling sites using the lc.dist function from the R (R Core Team, 2020) package marmap v1.0 (Pante & Simon-Bouhet, 2013).

Quality control

Some 679 individuals out of the 1,848 genotyped ones were discarded due to exceeding the threshold of acceptance of missing genotypic data, which was set at 25% or due to unreliable scoring. From the 1,169 retained individuals, no missing markers were reported for 84% of the fish, whereas only 4% of the individuals displayed the maximum missing data allowed (see distribution of missing data per sample and locus in Fig. S1 in Supporting Information).

Data validation conducted with Micro-Checker did not find evidence of scoring error due to stuttering or large allelic dropout but suggested the possibility of null alleles in two of the loci, which showed a significant deficit of heterozygotes (p < 0.0001) after global Hardy-Weinberg tests. The frequency of null alleles per sample in these loci took a maximum value of 0.216 (locus STJ5 in ANG_2), whereas it was <0.05 in 71% of the cases (see Table S2 in Supporting information). As global and pairwise divergence calculated with and without using the ENA correction implemented in FreeNA displayed similar values (i.e., a difference of 0.00045 in the overall F_ST values); the two loci suggestive of displaying null alleles were retained for further analyses.

The simulation analyses revealed that the dataset carried enough statistical power to detect genetic differentiation as low as F_ST = 0.0066 in 80% of the cases, whereas only in 11% of the cases when levels of differentiation were five-fold lower (F_ST = 0.0012). The probability of rejecting the null hypothesis of panmixia if true (type I error) was estimated at 5% (Fisher test) and thus considered appropriate (Fig. S2 in Supporting Information). In addition, the number of microsatellite markers with the capacity of differentiating unique individuals can be inferred from the plateau of the genotype accumulation curve and this was achieved with some 37% of the markers used (Fig. S3 in Supporting Information).

Summary statistics

The eight loci showed high rates of amplification success (93–100%) and suitable levels of polymorphism for population genetic analysis (100% polymorphic loci per sample). Allelic richness based on a minimum sample size of 20 diploid individuals took similar values throughout all the geographic range and moved between 11.5 (IVCO) and 13.6 (GAB_1), whereas the number of alleles for the same set of samples ranged between 96 and 137, respectively (Table 1). A total of 17 private alleles were found; they were distributed across all loci but SJT122 and ranged from 1 (Sco2-1, SJ78) to 4 (SJNT19, SJT199 and SJT53) and across 12 samples ranging from 1 (MOR_6, MOR_7, MOR_9, MAU_4, GAM, SIER, IVCO,GHA_1, GAB_1) to 4 (SEN_1). H_o and uH_E took similar values across samples and in the range 0.73–0.9. uH_E was consistently higher than H_O in all samples; in consequence, multilocus test detected heterozygote deficit in every sample but SIER (p-val = 0.257). The three loci displaying deviations from HWE in a larger number of samples after correcting for FDR were: SJT5 (in 20 samples), and Sco2-1 and SJ78 (in six). Both SJT5 and Sco2-1 exhibited homozygote excess and likely presence of null alleles.

Genetic differentiation

No locus was found to be deviate from neutral expectations according to Arlequin, whereas BayeScan flagged all of them as putative candidates to balancing selection (Table S3 in Supporting Information). Due to the lack of consensus of neutrality between methods, the full set of microsatellites was retained for subsequent analyses.

Overall, no significant population structuring was detected across the study area as indicated by a global F_ST of 0.001 (p = 0.587). Pairwise F_ST ranged from 0.000 to 0.011 (Table 2) with no significant differentiation whatsoever across the investigated geographic range.

STRUCTURE conducted both with and without LOCPRIORS reported a decreasing trend of LnP(D) across subsequent values of K with the highest average likelihood at K = 1 (Figs. S4A and S4B in Supporting Information). Evanno test, which by definition cannot yield K = 1, reported K = 2 (no priors) and K = 5 (LOCPRIORS), respectively but with extremely little support as both ΔK values were low (ΔK < 7). Puechmaille (2016)’s four statistics (MedMedK, MedMeanK, MaxMedK and MaxMeanK) implemented in StructureSelector confirmed one single genetic group in the model using priors, and K = 2 when not using them (Figs. S4A and S4B in Supporting Information). The barplot for K = 2 for the model without using LOCPRIORS revealed inferred ancestry to cluster ranging around 0.5 for each individual (Fig. S5 in Supporting Information), which is consistent with one single genetic unit.

Cross validation determined that 200 was the optimal number of PCs to be retained for the DAPC analysis using the 33 geographically explicit locations. In agreement with STRUCTURE, the corresponding DAPC plot revealed a large overlap of the individuals sampled across all the studied geographic range (Fig. 2) with none of the axes reaching 7% of the variation. Although some some individuals of SEN_1 partially deviated from the cloud, SEN_1 was not significantly different from any of the remaining samples according to pairwise F_ST.

Mantel test revealed lack of correlation between genetic differentiation (pairwise F_ST) and geographic distance along the latitudinal gradient ranging from Morocco to Namibia (r_xy = 0.012, p = 0.401).

Genetic diversity

Many marine fishes exhibit high levels of genetic diversity due to historically large effective population sizes (Hauser et al., 2002; Hauser & Carvalho, 2008). All samples analysed in the current study showed high levels of genetic diversity in terms of allelic richness (11.5 to 13.6) and heterozygosity (≥0.7) in agreement with results obtained for other medium and small pelagic species, e.g., S. australasicus in the western North Pacific (Tzeng et al., 2009), European sprat Sprattus sprattus (Glover et al., 2011), European anchovy Engraulis encrasicolus (Zarraonaindia et al., 2009) or European pilchard Sardina pilchardus (Baibai et al., 2012) in the North-East Atlantic. Overall, all samples (with the exception of the one from Sierra Leone represented by just eleven individuals) deviated from HWE expectations showing a deficit of heterozygotes. Heterozygote deficiency is not new to Scomber literature as it has been shown in microsatellite studies of S. colias in the East Atlantic Ocean and the Mediterranean Sea (Medina-Alcaraz, 2016) as well as in S. japonicus along the Chinese coast (Cheng, Zhu & Chen, 2014; Zeng, Cheng & Chen, 2012). Heterozygote deficiency in non-selfing, diploid populations is relatively common (Waldman & McKinnon, 1993) and in largely outcrossing species, it has commonly been interpreted as indirect evidence of the mixing of differentiated gene pools (Wahlund effect). However, selection and the presence of null alleles are common cause of heterozygote deficiency. In agreement with the latter, in the current study, two out of the eight genotyped loci were likely to display null alleles according to Micro-Checker assessment. Finally, and given the S. colias stocks are considered to be overexploited along West Africa (FAO, 2020), the third hypothesis to explain heterozygote deficiency could be overexploitation as overfishing reduces population size drastically, leading to genetic bottlenecks and loss of alleles (Smith, 1994).

Genetic population structure

In this study, the lack of genetic differentiation in S. colias revealed by eight microsatellites in the geographic area from Morocco (27.39°N) to Namibia (22.21°S) supporting a theory of panmixia would be a foreseen outcome due to pertinent ecological features of pelagic fish, such as high fecundity, large effective population sizes that impose limitations to genetic drift as well as a notable dispersal capacity as planktonic larvae and free-swimming adults that promotes gene flow (e.g., Carlsson et al., 2004; Martínez et al., 2006; Waples, 1987, 1998). The combination of both features acts as a major homogenizing force hindering genetic differentiation and eventually leading to panmixia (Canales-Aguirre et al., 2018; De Bie et al., 2012; Roy et al., 2014; Weersing & Toonen, 2009).

S. colias is a highly migratory species with individuals capable of covering extensive distances during their annual migrations favoured by a body shape well fitted for swift motion and free swimming mode of life (Collette, 1999; Uriarte et al., 2001). In addition, its planktonic larval stages also display considerable dispersal capabilities as pelagic larvae can be transported over vast distances by ocean currents and eddies (Hernández-León, Gómez & Arístegui, 2007) effectively leading to genetic homogeneity among populations. In agreement, other studies have revealed a notable absence of genetic differentiation on S. colias in the Northeast Atlantic and adjacent waters suggesting a large panmictic unit (Rodríguez-Ezpeleta et al., 2016; Stroganov et al., 2023; Zardoya et al., 2004). Likewise, congeneric species such as S. scombrus reveal population differentiation only at a transatlantic scale (Nesbø et al., 2000), whereas S. japonicus displays weak genetic structure across the NW Pacific (Cheng et al., 2015).

The hydrological conditions along the West African coast play a crucial role in various ecological processes and are influenced by factors such as ocean currents, eddies, temperature gradients, and upwelling phenomena. The region experiences the influence of several major ocean currents (see Fig. 1), including the Canary Current, the Guinea Current, and the Benguela Current, which affect the distribution of nutrients and plankton, which in turn impacts the abundance and distribution of marine life, including fish populations and seem to shape genetic differentiation of coastal species in this region (Alpers et al., 2013; Barton, Field & Roy, 2013; Mittelstaedt, 1991; Nielsen et al., 2018; Pelegrí & Peña-Izquierdo, 2015; Shannon, 2001). In consequence, distant sites may be connected by a strong current between them, so, the absence of the genetic structure that can result from the long-distance dispersal of the early development stages (Vásquez et al., 2013; White et al., 2010).

However, regional differences in various life-history traits with latitudinal trends have been detected in both European and NW African waters, with the Strait of Gibraltar as changing point (ICES, 2021). Several studies have addressed the delineation of S. colias population structure in the East Atlantic using phenotypic characters, and some latitudinal trends seem to occur in length composition and life-history traits from the Strait of Gibraltar northwards into European waters and southwards, through African waters (Domínguez-Petit et al., 2022; ICES, 2021; Landa et al., 2022; Perrotta, Carvalho & Isidro, 2005). Differences include size at first maturity, which seems to be larger in the South of Africa than in the North (ICES, 2021). Likewise, the spawning period, a crucial evolutionary driver in pelagic species (e.g., Lamichhaney et al., 2017; Petrou et al., 2021), corresponds to January-March in the region comprising Morocco to Senegal, in contrast to June-September in the area from Ghana to South Africa (ICES, 2021). Regarding the Atlantic coast of North America, the population structure of the Atlantic chub mackerel has not been extensively studied; however, differences in growth (Daley & Leaf, 2019), habitat suitability (Chen et al., 2009), spawning seasons (Weber & McClatchie, 2011), size at maturity (Cerna & Plaza, 2014) and morphology (Erguden et al., 2009) suggest that different sub-stocks might exist. The latitudinal trends observed in the East Atlantic, with the Strait of Gibraltar as inflection point (ICES, 2021), indicate that significant differences might be expected if individuals from the whole East Atlantic distribution of S. colias (including Mediterranean samples) were jointly analysed. This study is lacking samples from North of Morocco as well as from the Iberian Peninsula, which should be included in future genetic assessments.

Management implications

The lack of geographically explicit structure found here suggests the existence of one single putative genetic stock for S. colias inhabiting the West African waters from Morocco (27.39°N) to Namibia (22.21°S). Panmixia in a transboundary stock can present both opportunities and challenges. First, if a transboundary stock truly acts as a single panmictic population, then a single stock assessment might be sufficient, and the complexity of modelling population dynamics and setting catch quotas becomes simplified. In this case, all countries sharing the stock would benefit from implementing consistent management practices and should work together to set catch quotas and implement regulations that benefit the entire population.

However, a determination of panmixia based on a limited dataset might mask the presence of subtle subpopulations with distinct spawning grounds or migration patterns. In this case, management might need to consider specific population components, which if unnoticed could lead to inadequate management strategies and potential overexploitation of vulnerable subpopulations.

Although the sampling strategy of the present study comprises most of the West African geographical distribution of the species, this did not include the entire NE Atlantic distribution and the seasonal migrations could not be considered in the study. Therefore, although genetic divergences are not noticeable, some kind of local adaptations seem to exist, more when seasonal migration processes of the species have been described in the study area (García, 1982). From this perspective, and considering the current global warming scenario in which the northwards expansion of the species across the East Atlantic Ocean from regions of higher abundance off northwest Africa to the waters of the Atlantic Iberian and the Mediterranean Sea is a fact (Jurado-Ruzafa et al., 2024), fisheries management of Atlantic chub mackerel should account not only with punctual assessment statuses of political-based stocks but consider all the variables mentioned along the present section.

Although our data suggests that a nearly panmictic model is the most plausible scenario for the genetic population pattern of S. colias, implying high level of genetic connectivity across the Western African countries ranging from 28°N to 22°S. However, phenotypic variability and life history of S. colias, which follows a latitudinal gradient, support the evidence of regional demographic units as suggested by Sbiba et al. (2024) based upon otolith shape. Therefore, relying solely on genetic data may not provide a complete snapshot of the demographic connectivity, especially in wild marine fish populations characterized by high migration rates and large population sizes. It is crucial to consider demographic connectivity in addition to genetic connectivity when managing an important commercial species to ensure effective and sustainable measures. Finally, it cannot be neglected that genetic differentiation could have eluded the scrutiny of the set of markers used here as well as the possible impact of ascertainment bias driven by markers isolated from sampled from the Pacific ocean (Zeng & Cheng, 2012). Future studies would benefit from the genomic tools currently available such as the genome assembly of this species recently published (Machado et al., 2022).