Gars of the family Lepisosteidae are the sole survivors of an ancient group of ray-finned fishes that dates back to the late Jurassic, over 155 million years ago (Brito et al. 2017). There are seven extant species classified into two genera, Lepisosteus (4 species) and Atractosteus (3 species), which primarily occur in freshwater and coastal habitats of eastern North America, Central America, and Cuba (Berra 2001;Wright et al. 2012;Echelle and Grande 2014). All gars are carnivorous and have cylindrical, elongated bodies covered in ganoid scales. As is typical for sit-and-wait predators, their dorsal and anal fins are set far back on the body near the abbreviate heterocercal caudal fin, and they have long jaws that are densely packed with sharp teeth (Echelle and Grande 2014). In nature, gars are voracious predators that capture invertebrates, fish, and occasionally even small birds and mammals with sideway strikes of their head (Lauder 1980;Porter and Motta 2004;Lemberg et al. 2019).
All gars are large-bodied fish, with longnose gar (L. osseus), Florida gar (L. platyrhincus), and alligator gar (A. spatula) growing more than 2 m in length (García de León et al. 2001;Allan et al. 2005;Mendoza Alfaro et al. 2008;Buckmeier et al. 2016). Accordingly, these species play important ecological roles as the top predators in their habitats (Fry et al. 1999;Akin and Winemiller 2006;Williams and Trexler 2006;Fletcher et al. 2015). However, despite their considerable size, gars are not typically harvested commercially or for subsistence in the United States where most of the species occur. Rather, fishermen and fisheries managers have long targeted gars as nuisance species that were thought to negatively affect more valuable sport fish populations (Scarnecchia 1992). Small commercial fisheries, primarily for A. spatula, have only persisted regionally in some southern states (Aguilera et al. 2002), although there is a growing interest in gar as trophy fish that are angled for sport (Smith et al. 2020). Population declines and even local extirpation of gar caused by overfishing (i.e., nuisance animal removal) and habitat loss (Scarnecchia 1992;David et al. 2018), in conjunction with increased demands for trophy quality fish, have led to more widespread efforts to assess and manage gar populations (Glass et al. 2015;Bohn et al. 2017;Smith et al. 2018Smith et al. , 2020)).
An important exception with considerable economic importance is the tropical gar, Atractosteus tropicus (regionally known as catán, gaspar, machorra, and pejelagarto). The smallest member of the family (maximum length of ~ 1.25 m), A. tropicus occurs from southern México to Costa Rica (Fig. 1A), where it inhabits backwaters of large rivers, oxbow lakes, lagoons, flood ponds, and swamps, including some with brackish water (Villa 1982;Bussing 1998;Miller et al. 2005;Sauz-Sánchez et al. 2021). Especially in southern México and Guatemala, tropical gar has historically been a sought-after species of cultural importance (Emery 2017), and it is traditionally consumed because of its nutritional qualities (Manuel et al. 2013;Márquez-Couturier et al. 2015;Barrientos et al. 2018). In 1996 alone, over 530 tons of A. tropicus were harvested by fisherman setting gill nets in the Mexican state of Tabasco, making it one of the most important fisheries in the region (Aguilera et al. 2002). However, natural population sizes and catch rates of this species have been declining (Aguilera et al. 2002;Siebe et al. 2005), and some regional populations have even been extirpated (Villa 1982;Mora Jamett et al. 1997; Barrientos-Villalobos and Espinosa de los Monteros 2008). Accordingly, A. tropicus is listed as endangered in El Salvador (MARN 2015), but not in other countries. Population loss and decline have been attributed to overexploitation, habitat degradation, and habitat loss (Siebe et al. 2005;Sauz-Sánchez et al. 2021). Although systematic studies assessing the health of natural A. tropicus populations are missing and the fishery of this species remains unmanaged, there is also a growing interest in its aquaculture production (Mendoza Alfaro et al. 2008;Márquez-Couturier et al. 2015). Aquaculture of A. tropicus is facilitated by its high tolerance to nitrogenous wastes (Aranda-Morales et al. 2021), its ability to breathe atmospheric air and tolerate hypoxia (Burggren et al. 2016;Martínez-Bautista et al. 2022), a thorough understanding of its reproductive biology (Márquez-Couturier et al. 2015), its acceptance of artificial feeds (including Artemia nauplii for larvae and juveniles and pellets for adults; Márquez-Couturier et al. 2015;Martínez-Cárdenas et al. 2018), and its excellent food conversion rates and profitability (Palma-Cancino et al. 2019).
The increased pressure on natural populations of A. tropicus in conjunction with its increasing use in aquaculture production poses potential conflicts. While increased aquaculture production can relieve fishing pressure on declining natural populations, sharing of brood stocks across aquaculture cooperatives in different geographic regions has the potential to interfere with the natural population genetic structure and patterns of local adaptation (Ward 2006;Östergren et al. 2021). For example, a previous study on mitochondrial DNA (mtDNA) variation in A. tropicus has revealed significant differentiation among populations (Barrientos-Villalobos and Espinosa de los Monteros 2008). Specifically, analyses of mtDNA haplotypes indicated a differentiation between populations in river drainages flowing eastward toward the Gulf of México and the Caribbean Sea (from México to Costa Rica) and populations in river drainages flowing westward into the Pacific Ocean (in México and El Salvador; see cladogram in Fig. 1B). In addition, the identification of divergent and reciprocally monophyletic mtDNA haplotypes from Guatemala has fueled some speculation about the presence of a potentially undescribed species (Barrientos-Villalobos and Espinosa de los Monteros 2008).
In this study, we used nuclear markers (microsatellites) to revisit patterns of population genetic variation in A. tropicus. We used samples from nine natural populations in three geographic regions that span the natural distribution of the species-including samples from the Grijalva and Usumacinta River basins that drain into the Gulf of México, the San Juan River basin that drains into the Caribbean Sea, and multiple drainages along the Pacific versant of Central America-to test whether population genetic structure is in concordance with geographic barriers and reflects previous analyses of mtDNA. We predicted significant population genetic differentiation between populations from the Gulf of México/ Caribbean slope of Central America and those that occur on the Pacific versant (Barrientos-Villalobos and Espinosa de los Monteros 2008). Since tropical gar are increasingly farmed in different regions of Mexico, we also investigated the provenance of cultured A. tropicus from an aquaculture cooperative on the Pacific coast of México to assess the potential for human-mediated movement of gar across biographic boundaries.
Samples of tropical gar were collected between 2009 and 2017 from 10 locations that span much of the range of the species, including 9 natural populations and 1 population associated with an aquaculture facility (see Fig. 1 and Table 1 for an overview). On the Atlantic versant, sampling sites included the Laguna Chilapa in the Río Grijalva drainage, México (N = 32 individuals); the Río Chacamax in the Río Usumacinta drainage, México (N = 20); three proximate sites in the Laguna Canitzán and the main stem of the Río Usumacinta, also in the Río Usumacinta drainage, México (N = 22 + 9 + 16); and one site in the Laguna Caño Negro and nearby rivers, Río San Juan drainage, Costa Rica (N = 25). On the Pacific versant, sampling sites were in the Río Zanatenco, México (3 individuals); the Canal de Chiquimulilla associated with the Río Los Esclaves, Guatemala (N = 34); and a site in the Barra de Santiago wetland associated with Ríos Cara Sucia and San Pedro of El Salvador (N = 39). In addition, we obtained 12 samples from an aquaculture cooperative located in La Palma, Chiapas, on the Pacific versant of México. Hence, the total sample size included in this study was 212 individuals.
All individuals for this study were adults (> 60 cm total length) caught with cast nets. A small piece of the caudal or a pelvic fin was removed from each of the individuals with dissection scissors upon capture, and tissues were preserved in 95% ethanol for later processing in the lab. All individuals were immediately released at the original collection site once the tissue sample was secured.
The total genomic DNA was extracted from ethanolpreserved fin clips with the DNeasy Blood & Tissue Kit (Qiagen, Inc., Valencia, CA) following the manufacturer's protocol. Extracted DNA was used for polymerase chain reaction (PCR) amplification of 11 tri-nucleotide microsatellites using primers Asp007 (annealing temperature, T a = 56 °C), Asp035 (T a = 56 °C), Asp040 (T a = 56 °C), Asp053 (T a = 53 °C), Asp054 (T a = 56 °C), Asp057 (T a = 56 °C), Asp066 (T a = 56 °C), Asp072 (T a = 56 °C), Asp122 (T a = 60 °C), Asp159 (T a = 58 °C), and Asp168 (T a = 56 °C) (Moyer et al. 2009;Bohn et al. 2013). Microsatellites of all loci were amplified using PCR mix A from Moyer et al. (2009): 50 mM KCl, 10 mM Tris-HCl at pH 8.3, 0.01% gelatin, 2.0 mM MgCl 2 , 200 μM dNTPs, 0.1875 U Taq polymerase, 0.3 μM of both the forward and reverse primer, 20-100 ng of template DNA, and RNase-free water to a final volume of 12.5 μl. PCR conditions consisted of an initial denaturing at 94 °C for 2 min, and then 35 cycles of 94 °C for 30 s, 56-60 °C (depending on the primer specific annealing temperature) for 1 min, and 72 °C for 1 min. At the end, a 10-min final elongation step concluded the PCR. The PCR products were separated using capillary electrophoresis on an ABI 377 automated DNAsequencer (Applied Biosystems, Inc.). Genotypes were determined using Genotyper 2.1 software (Applied Biosystems), and Micro-checker v. 2.2.3 (Van Oosterhout et al. 2004) was used to assess the microsatellite data for null alleles and scoring errors.
Data analysis was conducted in R v. 4.1.1 (R Core Team 2021), unless otherwise noted. We summarized missing data using the info_table() function from the poppr package (Kamvar et al. 2014). We then tested for departures from Hardy-Weinberg equilibria using the hw.test() function from pegas (Paradis 2010); P-values were derived analytically from a χ 2 statistic and from a Monte Carlo test with 1000 permutations. We also tested for departures from linkage equilibrium using the pair.ia() function from poppr (Kamvar et al. 2014). This function calculates a standardized index of association, r̄d, which is less biased to variation in sample sizes compared to the classical index of association, I A (Agapow and Burt 2001). P-values for deviations from linkage equilibrium were based on 1000 permutations. In addition, we calculated a series of descriptive statistics, including allelic richness using the allelic.richness() function from the hierfstat package (Goudet 2005), the number of alleles per locus and the number of polymorphic loci within each population using the summary() function in poppr (Kamvar et al. 2014), the observed and expected heterozygosity using the basic.stats() function in hierfstat (Goudet 2005), and the inbreeding coefficient, F IS , using the basic.stats() function in hierfstat (Goudet 2005). We also quantified the number of private alleles in each population with the private_alleles() function from poppr (Kamvar et al. 2014). We used two complementary methods to describe and visualize clusters of genetically related populations and individuals. First, we calculated chord distances (Cavalli-Sforza and Edwards 1967) using the genet.dist() function from the hierfstat package (Goudet 2005) to generate a neighborjoining three using the nj() function from the ape package (Paradis and Schliep 2019). Second, we conducted a principal component analysis (PCA) using the dudi.pca() function from the ade4 package (Dray and Dufour 2007).
To identify population structure empirically without a priori designation of defined populations, we used the program STRU CTU RE (v. 2.3.4;Pritchard et al. 2000), and STRU CTU RE outputs were processes with the pophelper package in R (Francis 2017). Specifically, we identified the optimal number of genetically distinct clusters (K) using the Evanno method (Evanno et al. 2005), as implemented in the evan-noMethodStructure() function, and visualized individual cluster assignment probabilities using the plotQ() function.
Finally, differentiation between all population pairs was evaluated using F ST (Weir & Cockerham 1984), as implemented in the pairwise.WCfst() function from hierfstat (Goudet 2005). We tested for by distance among all natural populations (i.e., excluding the aquaculture population, P10) using a Mantel test with F ST as the dependent variable and geographic distance (in kilometers) as the independent variable. The Mantel test was conducted with 1000 permutations of the mantel() function from the ecodist package (Goslee and Urban 2007).
Note that we conducted all analyses with and without the aquaculture population (P10, because it is a non-natural population) and the population from Río Zanatenco, Chiapas (P7, because of low sample size). However, there were no qualitative differences in the outcome of analyses with or without these populations. Since both populations are interesting from a conservation management perspective, we only present the full analyses using all 10 populations.
A total of 212 individual fish were successfully genotyped for 11 microsatellite loci. The missing data rate across all loci and populations was less than 2.3%, and no single population or locus had a missing data rate over 7.7% (Table S1). All loci were polymorphic, and the number of alleles ranged from 3 (Asp072) to 18 (Asp066; Table S2). Global analyses including all populations revealed significant deviation from both Hardy-Weinberg (Table S3) and linkage equilibria (Figure S1), which is expected from spatially structured populations with non-random mating. However, populationspecific analyses found no consistent evidence for deviations from Hardy-Weinberg expectations for specific loci or populations (Table S3), or for deviations from linkage equilibrium (Table S4).
Compared across natural populations, allelic richness varied from 1.493 in the population from El Salvador to 1.971 in the population from Laguna Canitzán in the Río Usumacinta basin (Table 2). Across the same populations, observed heterozygosity varied from 0.209-0.411 and expected heterozygosity from 0.244-0.444 (Table 2). The number of polymorphic loci ranged from 5 to 11 (Table 2). The population from Costa Rica exhibited the highest number of private alleles (10), exceeding the range for all other populations (0-4). Finally, the population from the aquaculture facility exhibited the highest values for allelic richness (2.099) and observed and expected heterozygosity (0.459 and 0.490), exceeding the metrics of all natural populations. Despite these high levels of genetic diversity, the aquaculture population exhibited no private alleles. We used a neighbor-joining tree based on populationlevel allele frequencies and a principal component analysis (PCA) to describe and visualize clusters of genetically related populations and individuals. The neighbor-joining tree (Fig. 2A) indicated a clear division between populations from the Atlantic and Pacific versants of Central America. Within the Atlantic cluster, the population from Costa Rica was divergent from the populations in México. Most notably, the population from the aquaculture facility clustered with the Atlantic populations from México, even though the facility is located on México's Pacific versant. The PCA confirmed the patterns of the neighbor-joining tree and recovered three distinct clusters (Fig. 2B), with a clear separation of the Atlantic and Pacific populations along axis 1 and separation of the Costa Rican and Mexican populations along axis 2. Again, the aquaculture population clustered with populations from the Atlantic side of México.
Finally, the identification of three distinct population clusters was supported empirically by the STRU CTU RE analysis (best supported number of clusters K = 3; Table 3; Fig. 3). The three clusters included (1) populations from the Atlantic drainages of the Grijalva and Usumacinta in México, ( 2) the population from the Atlantic versant of Costa Rica, and (3) the populations from the Pacific drainages in México, Guatemala, and El Salvador. The aquaculture population once again clustered with populations from México's Atlantic drainages (Fig. 3).
Pairwise F ST values between population pairs ranged across two orders of magnitude, from 0.007 between two populations on México's Atlantic versant (P3 and P5) and 0.722 between a Mexican population on the Atlantic versant (P4) and the population from El Salvador (P9; Table 4). With a single exception (P3/P5), all pairwise F ST values were statistically significant. There was a significant correlation between geographic and genetic distance (Mantel test: r = 0.320, P = 0.003), indicating isolation by distance (Fig. 4).
The tropical gar (Atractosteus tropicus) is the smallest species within the family Lepisosteidae, yet this emergent aquaculture species has perhaps the highest socioeconomic impact both in terms of subsistence and commercial fisheries. Our study of the population genetic structure in this species revealed deep divergences among populations in different regions of México and Central America. We identified three genetic clusters that likely represent ancient genetic lineages with distinct geographic distributions. Most prominently, there was a divide between populations inhabiting rivers on the Atlantic versus the Pacific versant, a pattern commonly observed in other freshwater fishes of this region (McCafferty et al. 2012;McMahan et al. 2013; Picq et al. Note that all comparisons were significantly different from zero (P < 0.001), except for P3/P5 (marked with an asterisk). The populations IDs (P1-P10) correspond to 2014; Palacios et al. 2016;Ward et al. 2022). In addition, populations from the Atlantic versant were further split into a northern cluster from the Grijalva-Usumacinta River basins in México and a southern cluster from the Río San Juan basin in Costa Rica.
Our findings based on the analysis of nuclear microsatellite markers exhibit commonalities and differences with a prior population genetic study that used analyses of mitochondrial DNA (Barrientos-Villalobos and Espinosa de los Monteros 2008). Both studies recovered a clear split between populations from the Atlantic and the Pacific versants, reflecting extended periods of geographic isolation and a pattern of isolation by distance, which suggests limited gene flow. In contrast, there was no clear pattern of differentiation in mitochondrial haplotypes between Atlantic versant populations in México and Costa Rica. Instead, Barrientos-Villalobos and Espinosa de los Monteros (2008) documented deeply divergent mitochondrial haplotypes from the lower Río Usumacinta and the Petén region of Guatemala that formed a sister clade to all other tropical gar samples, including individuals from the Atlantic and Pacific versants (Fig. 1B), which led to speculation about a putatively undescribed species. However, analyses of nuclear markers in our study did not corroborate the cryptic species interpretation. Although we were not able to obtain samples from Guatemala's Petén region, there was no evidence for a divergent population in sites associated with the Río Usumacinta. Hence, the divergent mitochondrial haplotypes likely represent ancient polymorphisms that have been maintained within populations of tropical gar, rather than distinct evolutionary lineages. Future analyses will nonetheless require a more comprehensive coverage of the species' range to obtain a complete picture of genetically distinct lineages present within the species, including the Ríos Papaloapan and Coatzacoalcos in México, the Petén region in Guatemala, and occurrences around Lago Xolotlán and Lago Cocibolca in Nicaragua.
One of the key findings of our study was that the aquaculture population (located on the Pacific versant of México) was genetically derived from populations on the Atlantic versant, which belong to a different genetic cluster than more proximate populations in rivers that flow westward. This affinity to Atlantic versant populations was confirmed through the neighbor-joining tree, principal component analysis, and the analysis of population structure, suggesting that brood stock for use in aquaculture has been transported across biogeographic boundaries. Follow-up inquiries with the operators of the facility have confirmed that the brood stock at the aquaculture cooperative in La Palma was originally derived from fish produced in a facility from Tabasco (Arias-Rodriguez, personal communication). Transfer of fish across biogeographic boundaries is potentially problematic if we consider the aquacultural practices of the region.
Aquaculture cooperatives breed tropical gar in their facilities and then sell juveniles to landowners that stock them in ponds for rearing to adult size and eventual harvesting. However, individual fish can escape into natural systems, especially during flooding periods that are common in this region and trigger spawning in A. tropicus. This scenario has been documented in other species, especially cichlids of the genus Oreochromis, which are commonly used in aquaculture in the same region and have become invasive nuisance species as a consequence (Esselman et al. 2013). Interestingly, mixing of natural stocks may have occurred already. Two of three individuals genotyped from a natural population on México's Pacific versant (P7) exhibit a comparatively high putative ancestry from the Atlantic clusters. Whether this finding is a consequence of human activity remains to be tested, especially using higher sample sizes from surrounding populations. At least in theory, the observed pattern may also be a consequence natural gene flow mediated through stream capture in the intersecting river systems of the upper Río Grijalva (although tropical gar does not typically occur in those habitats), or it may reflect ancestral genetic variation that is shared and has been maintained across geographic regions.
The sharing of brood stock among hatcheries could significantly impact the population genetic structure of tropical gar in the future, as new facilities producing gar are continuously appearing. Currently, there are at least six known facilities producing tropical gar in the Mexican state of Tabasco, two in Campeche, one in Veracruz (all associated with water bodies connected to the Gulf of Mexico), and two in Chiapas (both associated with Pacific drainages; Arias-Rodriguez, personal communication). One facility for which production numbers are available (at the Universidad Juárez Autónoma de Tabasco) sells about 200,000 individuals annually. The current practice is that brood stocks are recruited from laboratory-reared individuals and widely shared across facilities, because acclimating wild-caught fish to hatchery conditions is difficult and time consuming. Wild-caught tropical gar have to be fed with live fish, and transitioning them to artificial feeds is not always successful. In addition, wild-caught fish will only spawn in the rainy season even in captivity, which limits productivity. Using breeders recruited from laboratory-born individuals circumnavigates these problems as they readily accept artificial feeds and can be spawned year-round.
Moving forward, potential ecological and evolutionary differences among genetic clusters should be considered in the management of tropical gar populations and aquaculture practices. Aquaculture can be regulated at the local, regional, and international levels to avoid the trafficking and mixing of brood stock across biogeographic boundaries and maintain the genetic structure documented in our study. Such regulations would help to prevent introgression among deeply divergent genetic clusters that could represent evolutionary significant units with different fisheries and conservation management needs. However, for such regulations to be effective, there are several pressing questions that will need to be addressed through future research. From an aquaculture perspective, developing methods for the more effective recruitment of wild-caught individuals as breeders in captive populations should have high priority. Establishment of captive breeders from local populations will reduce the need of sharing brood stock across regions, and it will open opportunities to investigate how different stocks perform in captive care, with concomitant insights about potential differences between native and non-native bloodlines that could impact natural populations. Reliance on local brood stock would also facilitate the frequent turnover of breeding individuals to prevent hatchery adaptation, as documented in other freshwater fishes (Frankham 2008;Christie et al. 2012Christie et al. , 2014)). From a fisheries and conservation management perspective, it will be critical to investigate whether the different genetic clusters documented here also differ in aspects of their phenotypes, especially in terms of physiological, behavioral, and life history traits that could mediate local adaptation. These data will be critical to identify potential evolutionary significant units and designate conservation priorities. Future work will also need to investigate the extent and nature of potential introgression using genomic tools to better understand the risks associated with aquacultural releases. In the meantime, tropical gar farmers can be informed that introgression from non-local populations has the potential to reduce survival and reproductive rates, with concomitant negative effects for natural populations (Araki et al. 2008). Proactive fisheries and aquaculture management will hopefully help to counteract recent trends of population decline and local extirpation of tropical gar (Villa 1982;Mora Jamett et al. 1997;Aguilera et al. 2002;Siebe et al. 2005), such that future generations might continue to enjoy traditional pejelagarto dishes served in the region.