Introduction
All around the world increased fishing pressure has resulted in overexploitation of populations with declines in overall abundance of populations and average fish size along with loss of genetic diversity (Marchal et al., Reference Marchal, Ulrich, Korsbrekke, Pastoors and Rackham2003; Froese, Reference Froese2004; Pérez-Ruzafa et al., Reference Pérez-Ruzafa, González-Wangüemert, Lenfant, Marcos and García-Charton2006; Favoretto et al., Reference Favoretto, Mascareñas-Osorio, León-Deniz, González-Salas, Pérez-España, Rivera-Higueras, Ruiz-Zárate, Vega-Zepeda, Villegas-Hernández and Aburto-Oropeza2020). A major management strategy for fisheries and marine conservation has been the implementation of natural protected areas (NPAs), whose success is based on the connectivity between populations (Palumbi, Reference Palumbi2003; Galindo et al., Reference Galindo, Olson and Palumbi2006; Green et al., Reference Green, Maypa, Almany, Rhodes, Weeks, Abesamis, Gleason, Mumby and White2015; Salles et al., Reference Salles, Maynard, Joannides, Barbu, Saenz-Agudelo, Almany, Berumen, Thorrold, Jones and Planes2015), thus understanding the extent of connectivity between fragmented populations is integral to understanding the stability of the subpopulations and their patterns of dispersal (Cowen et al., Reference Cowen, Lwiza, Sponaugle, Paris and Olson2000, Reference Cowen, Paris and Srinivasan2006; Christie et al., Reference Christie, Johnson, Stallings and Hixon2010). Fisheries management tools can protect spawning biomass and intraspecific genetic diversity, thus maintaining population abundances and keeping fisheries healthy (Hellberg et al., Reference Hellberg, Burton, Neigel and Palumbi2002; Jones, Reference Jones2002; Ablan, Reference Ablan2006).
A population has been defined as a group of individuals of the same species that live together in an area of sufficient size to permit normal dispersal and/or migration behaviour and in which numerical changes are largely determined by birth and death processes (Berryman, Reference Berryman2002). In this regard, population structure has been acknowledged as essential to infer potential evolution and extinction risk derived from over-exploitation and inappropriate management of marine resources (Ihssen et al., Reference Ihssen, Booke, Casselman, McGlade, Payne and Utter1981; Cadrin et al., Reference Cadrin, Kerr and Mariani2014; Cadrin, Reference Cadrin2020). Different techniques have been used for population discrimination (or stock identification) in fisheries, including the use of parasites (MacKenzie et al., Reference MacKenzie, Campbell, Mattiucci, Ramos, Pinto and Abaunza2008; Lester & MacKenzie, Reference Lester and MacKenzie2009; Baldwin et al., Reference Baldwin, Banks and Jacobson2012; Vasconcelos et al., Reference Vasconcelos, Hermida, Saraiva, González and Gordo2017), comparison of patterns of morphometric variation in fish (Murta, Reference Murta2000; Turan et al., Reference Turan, Oral, Ozturk and Duzgunes2006); differences in the morphology (Tracey et al., Reference Tracey, Lyle and Duhamel2006; Stransky et al., Reference Stransky, Baumann, Fevolden, Harbitz, Høie, Nedreaas, Salberg and Skarstein2008a, Reference Stransky, Murta, Schlickeisen and Zimmermann2008b) />and chemical composition of otoliths (Campana, Reference Campana1999; Bergenius et al., Reference Bergenius, Mapstone, Begg and Murchie2005; Jónsdóttir et al., Reference Jónsdóttir, Campana and Marteinsdottir2006) and molecular analysis such as mitochondrial DNA or microsatellite DNA sequencing (Beacham et al., Reference Beacham, Wood, Withler, Le and Miller1999, Reference Beacham, Pollard and Le2000; Shaklee & Currens, Reference Shaklee, Currens and Hallerman2003; Corander et al., Reference Corander, Marttinen and Mäntyniemi2006). Despite progress in the different methodologies applicable to population identification, the fullest possible picture can be obtained in response to the ecological, evolutionary and operational requirements that the population structure may involve, thus a recommended protocol is the complementary application of genotypic, ecological and/or phenotypic-based approaches in a holistic manner to maximize the probability of defining population structure properly (Begg & Waldman, Reference Begg and Waldman1999; Abaunza et al., Reference Abaunza, Murta, Campbell, Cimmaruta, Comesaña, Dahle, Gallo, García Santamaría, Gordo, Iversen, MacKenzie, Magoulas, Mattiucci, Molloy, Nascetti, Pinto, Quinta, Ramos, Ruggi, Sanjuan, Santos, Stransky and Zimmermann2008; Pita et al., Reference Pita, Casey, Hawkins, Villarreal, Gutiérrez, Cabral, Carocci, Abaunza, Pascual and Presa2016).
Currently, the black grouper, Mycteroperca bonaci, is one of the species within the regional multispecies fishery of the groupers, but unfortunately remains as an unregulated resource in coastal waters of the northern coast of the Yucatan Peninsula (DOF, 2006). Although NPAs have been considered the most efficient means of maintaining any population (Jennings et al., Reference Jennings, Reynolds and Polunin1999; Palumbi, Reference Palumbi2003) the mechanisms underlying the differentiation (ergo reduced connectivity) are still not well understood for the black grouper within the region, especially taking into account the poor vagility of the black grouper (Farmer, Reference Farmer2009). Thus, we evaluated population discrimination with special emphasis on the spatial scales (whether islands or coastal locations) because Alacranes is a reef island while Dzilam and Celestun are coastal locations, as well as the effects of protection from fishing on the population structure because despite the fact that all locations are NPAs each one has different fishing protection levels. Dzilam has a state NPA category where fishing is allowed, Celestun is a Biosphere Reserve where fishing is allowed under certain parameters, meanwhile Alacranes is a National Park where fishing of this species is not carried out on a regular basis.
Our research group of the marine biology laboratory (LABIOMA) from the University Autonomous of Yucatan (UADY) had previously used parasites as biological tags and DNA microsatellites separately through collecting specimens of M. bonaci among three separated (155–250 km) localities (Celestun, Dzilam and Alacranes Reef). First, the parasite communities of M. bonaci revealed the existence of two distinct metazoan parasite communities, Celestun and Alacranes (Espínola-Novelo et al., Reference Espínola-Novelo, González-Salas, Guillén-Hernández and MacKenzie2015), however, the Dzilam sample had not been analysed until the present study. Then, by means of the hosts’ genetic similarities those results were investigated at genotypic level in a recently published study to shed light on the population structure of this species (González-Salas et al., Reference González-Salas, Villegas-Hernández, Poot-López, Pech-Puch, Guillén-Hernández and Barrera-Guzmán2020). However, the aim of the present study was a joint analysis of the parasite community and the molecular markers of their hosts, to facilitate population discrimination in the black grouper (M. bonaci) of the northern coast of the Yucatan peninsula focusing on the comparative discussion, highlighting the contradictions or coincidences between techniques to provide a final holistic overview. To this end, we combined the samples used for both approaches by including the parasites of the Dzilam sample (unpublished data) and by using the genotypic data of most of the hosts (150 out of 161 specimens) with similar size ranges.
Materials and methods
Study area and sampling
From July 2008 to February 2011, specimens of M. bonaci were collected at three geographically distinct locations in the northern coast of the Yucatan Peninsula (Figure 1). Celestun is located in the east region and Dzilam in the west region of the state of Yucatan; both areas are coastal and include two of the main ports of the state where fishing takes place by mostly artisanal fishermen. In contrast, Alacranes is a reef island located 130 km from the coast, where this species is usually caught by sportfishermen. A total of 161 individuals of M. bonaci were provided by the artisanal fishermen who caught them with hook and line or spear fishing on selected fishing grounds with depths from 1.5–20 m during the same period among the locations in order to avoid temporal bias. Immediately after capture, fish were transported on ice to the laboratory, where they were examined, measured (cm) and weighed (g).
Study area and location of three sampling sites of M. bonaci within each NPA in the northern coast of the Yucatan Peninsula. Alacranes is a reef island (National Park 130 km offshore), Dzilam (state NPA) and Celestun (Biosphere Reserve) are coastal locations. Dominant direction of ocean surface currents field in the Gulf of Mexico (determined with geostrophic currents by NOAA/AOML and NOAA/CoastWatch) are shown with dark arrows. Bathymetry isobaths (dashed lines) are shown (modified from the Mexican Navy Nautical Chart) at 10, 20, 30, 50 and 200 m of depth. Major coral reef areas (shaded light grey).
Parasite analysis
Entire fish were examined fresh for endoparasites. Recovered parasites were counted, preserved in 70% alcohol and processed for subsequent identification by means of conventional helminthological methods (Lamothe-Argumedo, Reference Lamothe-Argumedo1997). The overall black grouper sample (161 hosts) were distributed as follows: 63 at Alacranes (TL = 40.14 ± 8.25), 54 at Celestun (TL = 35.75 ± 7.52) and 44 at Dzilam (TL = 38.89 ± 6.27). The overall size of specimens used for the parasite approach ranged from 28.0–46 cm TL. According to Crabtree & Bullock (Reference Crabtree and Bullock1998) two age classes might occur, i.e. <42 cm (<2-years) and >42 cm (>2-years) in M. bonaci, however, only five fish were greater than 42 cm TL, thus most of the sample consisted of specimens of about the same age class (<2-years). Moreover, in a preliminary analysis similar values of abundance and richness of parasite taxa were observed between hosts of these two age classes, therefore the effect of fish length and age on parasite communities were disregarded. The comparative analyses between the infracommunities registered in each locality were carried out considering only the intestinal helminths (actually from whole digestive tract including intestine and stomach), in which case ecological parameters such as prevalence, mean intensity and mean abundance values were obtained per locality, and parasite communities were described at the infracommunity level following Bush et al. (Reference Bush, Lafferty, Lotz and Shostak1997).
In order to show quantitative similarity between localities, resemblance matrices among samples (hosts per locality) were estimated using Bray–Curtis similarity with parasite abundances, previously transformed by fourth root to reduce the importance of the dominant taxa and allowing low abundance taxa to become more relevant and reduce dispersion data (Anderson, Reference Anderson2001, Reference Anderson, Gorley and Clarke2005). Then, using the PRIMER package (Clarke & Gorley, Reference Clarke and Gorley2015) an analysis of similarity (ANOSIM) was used to determine if there was a significant difference in parasite assemblage structure between localities (overall and pairwise), in which case ANOSIM statistic (termed R) compares the mean of ranked dissimilarities between groups to the mean of ranked dissimilarities within groups, so that an R value close to ‘1’ suggests dissimilarity between groups while an R value close to ‘0’ suggests that dissimilarities are greater within groups than between groups (Clarke & Gorley, Reference Clarke and Gorley2015). A non-metric multidimensional scaling (MDS) was used to 2-D visualize differences in the parasite assemblages between localities, and a canonical discriminant analysis of principal coordinates (sub-routine PERMANOVA CAP for Bray–Curtis similarity) was used to identify the percentage of correctly classified individuals per sampling locality. Finally, a similarity percentage (SIMPER) analysis was also performed to identify the taxa mainly responsible for similarities and differences in parasite assemblage structure between localities.
Genetic analysis
In order to avoid possible genetic structure bias, only specimens whose size ranged between 30.0 and 45.0 cm TL were sampled to ensure that those individuals were about the same age and equal size sample, that is to say 150 (50 per sampling locality) out of the 161 hosts whose size ranges were about the same (28.0–48 cm TL). The molecular procedures (DNA extraction, PCR reactions) are detailed in González-Salas et al. (Reference González-Salas, Villegas-Hernández, Poot-López, Pech-Puch, Guillén-Hernández and Barrera-Guzmán2020), in which five microsatellite loci from the black grouper genome were amplified by PCR based on sequences in Zatcoff et al. (Reference Zatcoff, Ball and Chapman2002). Alleles were scored, and a matrix with the microsatellite genotype data was obtained with the number and size of alleles per locus by sampling sites.
As stated in González-Salas et al. (Reference González-Salas, Villegas-Hernández, Poot-López, Pech-Puch, Guillén-Hernández and Barrera-Guzmán2020), different measures of genetic variability, such as the number of alleles (Na), number of genotypes (Ng), observed (Ho) and expected heterozygosity (He) were calculated either per locus and multilocus per locality. Wright's inbreeding coefficient (F IS) was used to test whether the black grouper localities were in Hardy–Weinberg equilibrium, whereas Slatkin's linearized fixation index F ST (Slatkin, Reference Slatkin1995) was used as a measure of interpopulation variability. Analysis of molecular variance (AMOVA) was also used to assess the relative partitioning of genetic variation among localities both on a multilocus basis and locus by locus. Pairwise analyses of Nei's genetic distances (Ds) and Wright's gene flow (Nm) among localities were estimated. Additionally, the genetic structure was visualized in a multivariate space performing a Principal Coordinate Analysis (PCoA) on a matrix of Nei's individual pairwise genetic distances. Finally, a Bayesian clustering (STRUCTURE analysis) was estimated to identify any potential cryptic genetic variation within and among samples.
Results
Parasites of M. bonaci
A total of 21,632 individual intestinal parasites belonging to 15 taxa, 6 of which were identified to species level, were found in the 161 black groupers examined (Table 1). Thirteen taxa were found as adults and only 2 as larvae (Tetraphyllidea gen. sp. and Pseudoterranova sp.). Digeneans and nematodes were the groups with the highest number of taxa (6 each), followed with 2 acanthocephalans and only one taxon of cestodes. The greatest number of taxa were found at Celestun (12) followed by Alacranes (9) and Dzilam (9). Alacranes and Celestun, as well as Dzilam and Celestun, shared 47% (7 out of 15 taxa) of the parasites, while Alacranes and Dzilam shared 5 (33%) of the 15 taxa found in the intestine of this host species (Table 1).
Infection levels of the metazoan intestinal parasites, such as prevalence (P%), mean abundance (MA), and mean intensity (MI) of M. bonaci in three different localities of Yucatán, Mexico
N, number of hosts examined per locality.
Average richness and abundance of parasites per host were observed to be higher in Dzilam: 2.91 (±1.25) taxa and 349.32 (±620.86) total parasites in comparison with Celestun and Alacranes with lower values (Table 2). In Dzilam abundance was ranked first by Prosorhynchus atlanticus (digenean), secondly by Prosorhynchus sp. and thirdly by Gorgorhynchus sp. (acanthocephalan). In Alacranes, the dominant parasites were P. atlanticus, Hamacreadium mutabile and Prosorhynchus sp., all of them digeneans. Prosorhynchus atlanticus was also the most abundant taxa in Celestun, followed by Lepidapedoides epinepheli (digenean) and Prosorhynchus sp.
Average values of the parasite infracommunity descriptors from specimens of M. bonaci collected from Celestun, Alacranes and Dzilam, with standard deviations in parentheses. The rank of parasites with highest abundance per locality is also shown
Significant differences in the infracommunity compositions were found, indicating that parasite assemblages were structured locally (ANOSIM, global R = 0.289, 0.1% significance or P = 0.001, Figure 2A). Pair-wise comparisons revealed significant differences between all localities (Figure 2A), but Alacranes–Dzilam were apparently more dissimilar to each other (R = 0.447 and 0.1% significance), followed by Celestun–Dzilam (R = 0.398 and 0.1% significance), whereas Celestun–Dzilam were the least dissimilar ones (R = 0.094 and 0.1% significance). Non-metric MDS evidenced this parasite assemblage structure, where Dzilam is clearly differentiated in the multivariate space (Figure 2B). In terms of composition, the SIMPER analysis revealed that among the parasites found in M. bonaci, the taxa mainly responsible for the differences among localities were the digeneans Prosorhynchus atlanticus, Prosorhynchus sp., Lepidapedoides epinepheli and Hamacreadium mutabile as well as the acanthocephalan Gorgorhynchus sp. The comparison of the assemblage structure between Alacranes and Dzilam found the highest average dissimilarity (84.61), mainly due to three taxa: Prosorhynchus sp., Prosorhynchus atlanticus and Gorgorhynchus sp. (Figure 2C). Celestun and Dzilam showed an estimated 79.8 average dissimilarity, in this case the same three taxa contributed most to this dissimilarity. In contrast, the assemblage structure between Celestun and Alacranes was observed with a lower average dissimilarity (74.52), but along with Prosorhynchus atlanticus, another two taxa (Lepidapedoides epinepheli and Hamacreadium mutabile) contributed to the abundance differences between these localities (Figure 2C).
Graphical abstract of population structure of M. bonaci using intestinal parasites as biological tags for differentiation among localities based on the global and pairwise ANOSIM (A) with R (significance in percentage), MDS plot (B), and SIMPER (C). Symbols indicate position of individual samples in the multivariate space in the MDS.
Finally, discriminant analysis for Bray–Curtis similarity using the fourth root-transformed abundances showed significant differences between samples (tr = 0.7535, P <0.0001), with 76.40% of individuals correctly classified to their locality of origin (Table 3), in which 87.30% of individuals were correctly classified to Celestun, 62.96% in Alacranes and 77.27% to Dzilam.
Cross validation of the discriminant analysis of principal coordinates (sub-routine PERMANOVA CAP for Bray–Curtis similarity) with fourth root transformed abundances showing the number and percentage of hosts correctly classified to their locality of origin (Celestun, Alacranes and Dzilam) based on the parasite infracommunity
Genetics of M. bonaci
A total of 35 alleles (varying from 80 to 262 base pairs long) and 69 genotypes were found across all five loci and localities (Table 4). The greatest allelic richness was observed in Dzilam (Na = 30), followed by Celestun and Alacranes Reef with 29 alleles in both (Figure 3A), whereas the genotypic richness was observed to be very similar among localities: Dzilam (Ng = 49), Alacranes Reef (Ng = 52) and Celestun (Ng = 53). Both values of observed (Ho) and expected heterozygosity (He) reflected the high rate of heterozygosity that occurs in five microsatellite loci (Table 4), all of them giving similar estimates for each locality. In a multilocus manner all three expected heterozygosity values (0.82–0.85) were lower than observed values (Ho ≈ 1.0) (Table 4); however, it is worth mentioning that Dzilam showed the highest expected heterozygosity (0.85) in comparison with both Celestun and Alacranes (He = 0.82) (Figure 3A).
Graphical abstract of population structure of M. bonaci using the genetic connectivity for differentiation among localities based on the genotypic estimates (A) of number of alleles (Na), genotypic richness (Ng), observed heterozygosity (Ho), and expected heterozygosity (He). The pairwise Nei's distances (Ds), gene flow (Nm), and the fixation index (F ST) over loci at localities are also given (B), whose genotypic variations are also shown in the STRUCTURE (C) and the PCoA plots (D). Symbols indicate position of individual samples in the multivariate space in the PCoA.
Resume of genetic characteristics of the five microsatellite loci used for M. bonaci such as range of alleles (base pairs), number of alleles (Na) and genotypic richness (Ng) per locus
Observed heterozygosity (Ho), expected heterozygosity (He), inbreeding coefficient (F IS) and fixation index (F ST and its P-value) are also shown per locality and overall. All probability values indicated in bold were significant after false discovery rate. Multilocus estimates are also given for each locality.
Multilocus structure analysis showed a significant excess of heterozygotes (FIS = −0.141, P < 0.001) and a substructured population (F ST = 0.056, P < 0.0001), which suggests that it may not be panmictic. The pairwise F IS multilocus analysis showed that all three black grouper localities were significantly in Hardy–Weinberg disequilibrium (F IS = −0.1410, P < 0.001), in which case on a pairwise basis the higher F IS values were observed between Alacranes and Dzilam (F IS = −0.2007, P < 0.001) and lower values observed between Alacranes and Celestun (F IS = −0.1791, P < 0.001). Additionally, pairwise locality differentiation using the F ST indices showed that in fact there is a significant (P < 0.001) structuring among all three localities (Figure 3B). AMOVA confirmed this later by identifying a significant genetic structure among localities (F ST = 0.056, P < 0.0001) with 5.6% of the total genetic variation explained by the variation among localities. Nei's (Reference Nei1978) genetic distances (Ds) showed that Celestun and Dzilam were least similar localities to each other (Ds = 0.364) while the most similar ones were Celestun and Alacranes (Ds = 0.174) (Figure 3B). The lowest gene flow (number of effective migrants per generation) among pairs of localities was observed between Celestun and Dzilam (Nm = 4.29) and highest between Celestun and Alacranes (Nm = 8.42).
STRUCTURE analysis supported Dzilam as the most genetically distinct locality with most of its samples assigned around 80–90% to the same differentiated genetic cluster (Figure 3C), in contrast, all Celestun and Alacranes samples show an admixed genetic composition. PCoA confirmed the genetic structure in M. bonaci populations with the first two principal coordinates (Figure 3D) explaining most of the variation (14.78%). Dzilam samples stand out as the most differentiated, mainly along coordinate 1, while most Celestun and Alacranes samples overlap in the multivariate space.
Discussion
The present work represents the first study using both ecological and genotypic approaches used on the same specimens to explore the population structure of the black grouper Mycteroperca bonaci at a regional level in Gulf of Mexico, where three sampling sites (Celestun, Alacranes and Dzilam) were significantly different from each other.
The use of intestinal parasites as biological tags led to the first conclusion that a structured population of M. bonaci occurs in the northern coast of the Yucatan peninsula, since parasite infracommunities revealed the parasite assemblage structure between these distinct localities (Celestun, Alacranes and Dzilam), where the Dzilam sample was the most differentiated one, as observed altogether with the ANOSIM, MDS and SIMPER analyses. Intestinal helminth assemblages represent the available pool of parasite taxa a fish can acquire over time in a given location, thus differences between localities may be related to specific characteristics of the host ontogeny, or to particular biotic and abiotic environmental conditions (Poulin & Mouillot, Reference Poulin and Mouillot2004; Quiroz-Martínez & Salgado-Maldonado, Reference Quiroz-Martínez and Salgado-Maldonado2013). Therefore, hosts with high vagility might be expected to show high similarity values among infracommunities from different spatially separated locations, since high vagility may be an important determinant in the spread or exchange of parasites between host populations and thus determine similarity among parasite communities (Poulin, Reference Poulin2003). However, a reduced availability and abundance of suitable intermediate hosts to complete parasite life cycles might be key in the observed differences in the parasite infracommunity assemblages.
The parasite taxa mainly responsible for the differences among localities were the digeneans Prosorhynchus atlanticus, Prosorhynchus sp., Lepidapedoides epinepheli and Hamacreadium mutabile as well as the acanthocephalan Gorgorhynchus sp. These parasites may infect M. bonaci via the food chain, so differences between parasite assemblages may be explained by differences in the geographic distribution and/or abundance of organisms used as intermediate hosts by the parasites. We believe that the digeneans P. atlanticus and L. epinepheli may be short-lived intestinal parasites (probably between 10 and 18 months) that could be eliminated at the end of their lifespans. There is no clear relationship between the lifespan of the parasite and the characteristics of the host, but it seems that intestinal parasites from homeothermic hosts live longer than those from ectotherms (Kennedy, Reference Kennedy2006). Regarding this issue, some parasites (acanthocephalans and cestodes) are known to have attachment structures such as a proboscis which may leave evidence of their presence even after their death and expulsion from the host (de Buron & Nickol, Reference de Buron and Nickol1994). However, other helminths such as digeneans or nematodes do not have these kinds of structures and their detection within the host is determined only by their lifespan within it, which varies according to the characteristics of their life cycle, i.e. some intestinal parasites can live up to 2.5 years inside their bird host, but some live less as 2 months (Ellis & Williams, Reference Ellis and Williams1973). The latter period of permanence is shorter than those seen in the larval stages, which can remain in the host from host infection to its death (Williams et al., Reference Williams, MacKenzie and McCarthy1992; Marcogliese, Reference Marcogliese2004). In addition, since fish length has been shown to influence the richness and abundance of parasite infracommunities (Díaz & George-Nascimento, Reference Díaz and George-Nascimento2002; Flores & George-Nascimento, Reference Flores and George-Nascimento2009; Muñoz-Muga & Muñoz, Reference Muñoz-Muga and Muñoz2010), in our study comparisons between localities were made using hosts of the same length, to reduce the impact that may have in the discrimination between populations using these parasite species as biological markers (MacKenzie & Abaunza, Reference MacKenzie and Abaunza1998; MacKenzie, Reference MacKenzie2002).
Parasite infracommunities have been used previously as tools for host population discrimination, however, these studies used predominantly long-lived larval parasites, such as digenean metacercariae, cestode plerocercoids and juvenile acanthocephalans (George-Nascimento, Reference George-Nascimento2000; Sardella & Timi, Reference Sardella and Timi2004; Timi & Lanfranchi, Reference Timi and Lanfranchi2009), since the residence time in the host has been considered the most important factor in the use of parasites as biological tags (MacKenzie, Reference MacKenzie2002). In this regard, although ectoparasites (such as copepods and isopods) were also observed in M. bonaci they were excluded because the intestine was the site of infection where the largest number of different parasite taxa was found, besides ectoparasites (especially copepods) did not show differences in prevalence among locations. In addition, the acknowledged criteria for the use of parasites as biological tags, such as the geographic variation in prevalence or abundance along with longevity of infection and absence of reproduction directly in the host (Williams et al., Reference Williams, MacKenzie and McCarthy1992), encouraged us to use only intestinal helminth parasites excluding ectoparasites. Although the use of short-lived endoparasites (such as the tetraphyllideans) in population discrimination should be used with care (Lester & MacKenzie, Reference Lester and MacKenzie2009), when integrated with a complementary method, such as microsatellites or the chemical composition or the shape of otoliths, it can be a tool that reinforces the results (MacKenzie & Abaunza, Reference MacKenzie and Abaunza1998, Reference MacKenzie, Abaunza, Cadrin, Kerr and Mariani2014). For example, spatial discrimination of the deep-water redfish (Sebastes mentella) has been found within the Gulf of St. Lawrence and Newfoundland using fish molecular markers (Roques et al., Reference Roques, Sévigny and Bernatchez2002), and then by examining parasite communities host populations were further subdivided into four smaller groups (Marcogliese et al., Reference Marcogliese, Albert, Gagnon and Sévigny2003). In the latter study the authors only used three species, including an ectoparasite copepod (Sphyrion lumpi) along with the nematode larvae of Anisakis simplex and Hysterothylacium aduncum, since these fulfilled the criteria as biological tags. In this sense, in our study ectoparasites did not fulfil this criteria, nonetheless intestinal parasites did, and in fact detected similar spatial population discrimination, being that Celestun and Alacranes were the least dissimilar to each other, as also shown with the genetics of the hosts. As seen in our results, the lowest percentage of correct allocation of fish to their localities, based on their parasites, was recorded in Alacranes (62.96%), where 33.33% of the remaining host would be allocated instead in Celestun (Table 3). The genetic results prove that Alacranes, together with Celestun, were the localities with the highest allelic and genotypic richness as well as the pairwise highest gene flow and lowest genetic distance between them. It is likely that the connectivity between these two localities is driven by the marine currents of the Gulf of Mexico (Figure 1) enhancing therefore higher black grouper migration-emigration processes than those between the other two localities (Alacranes and Dzilam).
Population discrimination using parasites as biological tags represents a process in ecological time in all the necessary host and environmental conditions for parasite differences between localities, whereas differences at the host's genetic level would represent long-term changes in evolutionary time. Therefore, it seems that the environmental conditions and the absence/presence of the hosts necessary to complete the parasite life cycles, regulate their presence and therefore the composition of parasite communities present in fish from all locations (Thieltges et al., Reference Thieltges, Dolch, Krakau and Poulin2010), as well as the poor vagility of the black grouper which could be influencing the composition of parasite communities (Farmer, Reference Farmer2009). Either way, the use of parasites confirmed the sensitivity of this technique in the detection of different hosts’ population subunits, which would help to confirm the genotypic-based population discrimination of the black grouper by integrating our findings and reducing the type I error in identifying its subpopulations.
Similar to the parasitological results, the population discrimination based on the genotypic approach indicated that M. bonaci showed an excess of heterozygotes within each locality and the estimates of the F IS (−0.141) and F ST (0.056) suggested that all three localities deviated significantly from Hardy–Weinberg equilibrium and that there was a structuring within the study area, suggesting that the black grouper population cannot be considered a panmictic. In addition, heterozygote excess, usually revealed by negative inbreeding coefficient (F IS) such as ours, has several potential causes. For instance it may result from small reproductive population size (Pudovkin et al., Reference Pudovkin, Zaykin and Hedgecock1996; Coscia et al., Reference Coscia, Chopelet, Waples, Mann and Mariani2016), small sexual or self-incompatible populations (Balloux et al., Reference Balloux, Amos and Coulson2004), negative assortative mating (also known as disassortative mating or heterogamy) between individuals carrying different alleles (Stoeckel et al., Reference Stoeckel, Grange, Fernández-Manjarres, Bilger, Frascaria-Lacoste and Mariette2006). That is to say individuals with dissimilar genotypes or phenotypes mate with one another more frequently than would be expected under random mating, which in turn would reduce the genetic similarities within the population.
Significant genetic differences among localities showed by the multilocus AMOVA, then supported by the PCoA and STRUCTURE analyses, highlighted a gradient-like arrangement of all the localities and a singular genetic distinctness of M. bonaci samples from Dzilam. Despite Dzilam and Alacranes being geographically the closest (~155 km) localities, genetically the least distant ones were Celestun and Alacranes, whose small genetic differentiation might be explained by the non-linear distribution (but triangle shaped) of the sampling sites as well as the north-western flow (known as the Yucatan Current) on the western side of the Yucatan peninsula (Merino, Reference Merino1997). The Yucatan Current flows into the Gulf of Mexico through the Yucatan Channel and it eventually separates from the Campeche Bank and becomes the Loop Current forced by the easterly winds (Figure 1), which primarily prevail throughout most of the year and the hydrodynamics (surface currents) of the shelf have been mostly associated with this easterly wind stress pattern commonly observed in the entire region with little seasonal and spatial (depth) variation (Molinari & Morrison, Reference Molinari and Morrison1988; Enriquez et al., Reference Enriquez, Mariño-Tapia and Herrera-Silveira2010; Ruiz-Castillo et al., Reference Ruiz-Castillo, Gomez-Valdes, Sheinbaum and Rioja-Nieto2016). Ocean currents are known to have a major influence on larval dispersion with genetic consequences in fish populations (Cowen et al., Reference Cowen, Lwiza, Sponaugle, Paris and Olson2000, Reference Cowen, Paris and Srinivasan2006; Purcell et al., Reference Purcell, Cowen, Hughes and Williams2006). In this regard, a recent study (Johnston & Bernard, Reference Johnston and Bernard2017) showed that the larval dispersal of several species (lionfish, red grouper and marine species in general) decreases markedly from 89°W (close to Dzilam), seriously delaying the recruitment in the area between 91.1–91.5°W (close to Celestun), which would help explain our results in terms of genetic connectivity.
Deviations from Hardy–Weinberg equilibrium at each locality have been associated with inbreeding due to low larval dispersal in most marine species (Planes et al., Reference Planes, Parroni and Chauvet1998; Taylor & Hellberg, Reference Taylor and Hellberg2003; Purcell et al., Reference Purcell, Cowen, Hughes and Williams2006). Doherty et al. (Reference Doherty, Planes and Mather1995), as well as Riginos & Victor (Reference Riginos and Victor2001), found that fish species with long planktonic larval duration (> 40 days) usually have low levels of genetic differentiation among localities because their larvae have a greater chance of migration and gene flow, and therefore may contribute to the population gene pool. In this regard, the genus Mycteroperca normally presents an egg phase of 2 days, and a larval phase that lasts on average 43 days (Keener et al., Reference Keener, Johnson, Stender, Brothers and Beatty1988), however, given the oceanographic characteristics found by Johnston & Bernard (Reference Johnston and Bernard2017) for the Yucatan coast, its dispersion can indeed be limited by the restricted dispersion reported in the area, limiting the connectivity between populations, which suggests that self-recruitment is responsible for maintaining this locality with a limited genetic flow. Thus, the population structure observed in this study could indeed be due to local larval retention and a limited gene flow between localities, especially for Dzilam.
Regional fisheries might threaten not only the targeted population, but also entire communities of organisms, including their parasites, by reducing host abundance and food web complexity, which affect their transmission efficiency, population density and assemblage taxa richness (Wood & Lafferty, Reference Wood and Lafferty2015; Braicovich et al., Reference Braicovich, Irigoitia, Bovcon and Timi2021). Furthermore, it is important to preserve black groupers from Dzilam since they may act as reservoirs for rare genotypes (7.2% were unique), which gives a certain degree of distinctness to this particular subunit. Fishing is known to be one of the causes of reduced heterozygosity (Bergh & Getz, Reference Bergh and Getz1989). For example, fishing authorities should determine catch quotas or even regulate the fishing effort. In addition, in order to define a definitive spatial structure for the further population assessment for the regional fishery of this species, as shown in Figure 4, more information related to its distribution, dispersal, geographic variations and boundaries, needs to be gathered in line with recommendations by Cadrin (Reference Cadrin2020), along with temporal stability of spatial differences over multiple years.
Overall highlights of the exploratory population discrimination for Mycteroperca bonaci from the northern coast of the Yucatan peninsula in the present study using parasite assemblages and host genetics, and further challenges for the population assessment, in line with best practice suggested by Cadrin (Reference Cadrin2020), in a previous proposal for management strategies for this fishery.
NPAs, where samples of the black grouper were collected, have distinctive features. For instance, Dzilam was acknowledged as a State NPA in January 1989 and controlled by a dependency of the State government called Secretary of Sustainable Development. In the management plan of this NPA, M. bonaci is listed as a commercial fish. In contrast, the NPA of Celestun is known as a Biosphere Reserve controlled since November 2000 by a dependency of the Federal government called CONANP (the acronym in Spanish of Comisión Nacional de Áreas Naturales Protegidas). Its management plan only mentions the presence of 140 species of fish in the NPA, in which the red grouper (E. morio) is mentioned as a commercially important species and heavily exploited, but M. bonaci is not mentioned at all. Finally, the NPA of Alacranes, known as a National Park also controlled by the CONANP since June 1994, is located 140 km from the northern coast of the Yucatan Peninsula, in which M. bonaci is one of the species that occurs in the area, but the management plant of this NPA does not provide fishing regulations specifically for the black grouper, instead the fishing regulations declared in the management plan are only oriented to limit the number of fishermen and boats to certain areas (core zone or no-take zone) where fishing is not allowed (in comparison with Celestun and Dzilam where there are no-take zones). However, due to its insular nature and distance from the coast, the fishing pressure even in the areas where the extraction of some organisms is allowed (buffer zone and zone of sustainable use of natural resources) has not been particularly high (except in the summer months when this activity is intensified mainly by sport fishing). In these latter areas, the fishing regulation is aimed only at the use of low-impact fishing gears (e.g. hooks, traps, or collection by free diving and spearfishing), in addition to respecting the provisions for the temporary ban on groupers, which have a closed season for fishing (from 15 February to 15 March) in the State of Yucatan but this is related to the main target species (red grouper Epinephelus morio) (DOF, 2014), and the spawning season for M. bonaci usually extends from December to March (Brulé et al., Reference Brulé, Renán, Colés-Marrufo, Hauyon, Tuz-Sulub and Déniel2003). The present recognition of population subunits would bring beneficial effects to preserve more important gene pool variability of the coastal Dzilam subunit where fishing regulations are the least restrictive compared with the other studied NPAs.
Regional currents are expected to allow dispersion between distant localities to form a single panmictic population, however, both approaches (parasites and genetics) used in our study were consistent in the overall population discrimination. Moreover, the present study shows that although island NPAs (such as Alacranes) are expected to have some degree of distinctness, either ecological (Moore et al., Reference Moore, Buckworth, Moss and Lester2003; Lester & Moore, Reference Lester and Moore2015) or genetic (Ruzzante et al., Reference Ruzzante, Taggart, Cook and Goddard1997; Hemmer-Hansen et al., Reference Hemmer-Hansen, Nielsen, Grønkjær and Loeschcke2007), our results suggest that even though the black grouper population is spatially structured into at least three units (Celestun, Dzilam and Alacranes), in this case we identified the coastal locality of Dzilam as the more differentiated (genetically and lowest parasite diversity) probably caused by oceanographic-driven low levels of larval dispersal and adults' migration.
Data
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Acknowledgements
Special thanks to all the anglers from Dzilam, Alacranes Reef and Celestun. This study was financed by PROMEP 2007 and Universidad Autónoma de Yucatán (UADY). Dawrin Pech-Puch received his postdoctoral fellowship from the National Council of Science and Technology (CONACYT) of Mexico.
