Genetic structure and historical and contemporary gene flow of Astyanaxmexicanus in the Gulf of Mexico slope: a microsatellite-based analysis

Sampling

The present study included samples of Astyanax spp. from at least seven of the main hydrographic basins of the Gulf of Mexico slope in Mexico: The Bravo, San Fernando, Soto la Marina, Pánuco, Tuxpan, Papaloapan, and Grijalva-Usumacinta River basins (Table 1; Fig. 1). Fish specimens were collected using electrofishing equipment and trawl nets. Tissue samples (fin clips) were preserved in 96% ethanol for subsequent DNA extraction. A total of 469 individuals from 16 localities, corresponding to 4 cave and 12 surface locations, were included. Collection of the samples was carried out with the approval of the Mexican government (DGOPA.05003.181010-5003; DGOPA.00570.288108-0291; DAPA/2/130409/0961, 230401-613-03).

DNA extraction, PCR and genotyping

Total DNA was extracted from fin tissue using a DNeasy blood and tissue kit from Qiagen. Ten polymorphic microsatellite loci developed for the genus Astyanax (Strecker, 2003; Protas et al., 2006) were studied (Table S1) through amplification in 15 µl reactions, following the cycling conditions shown in Table S2. Individual genotyping was obtained using 6% denaturing polyacrylamide electrophoresis gels, run in a vertical electrophoresis unit at 1,700 V, 40A and 45 W for 2.5 h. DNA fragments were visualized by silver staining (Bassam & Caetano-Anollés, 1993). Allelic size was obtained with gels using the Sequentix Gel Analyzer (http://www.sequentix.de/).

Statistical analyses

In order to detect the presence of null alleles at every locus and in each location, the Brookfield 1 Index (Brookfield, 1996) and the Expectation-Maximization algorithm (Dempster, Laird & Rubin, 1977) were applied using MICRO-CHECKER version 2.2.3 (Oosterhout, Wills & Hutchinson, 2004) and FreeNA (Chapuis & Estoup, 2007), respectively. Linkage disequilibrium and significant deviation from Hardy-Weinberg (HW) expectations using an exact test (10,000 dememorization steps, 5,000 batches of 10,000 iterations per batch) were determined with GENPOP (Raymond & Rousset, 2001). Statistical significance values were corrected using the sequential Bonferroni method (Rice, 1989), with a significance value of 5%.

Genetic variability within populations (as estimated by STRUCTURE) was evaluated by the number of alleles (N), number of alleles per locus (Na), effective number of alleles (Nae), observed heterozygosity (H_O), and the fixation index (F_IS), using GENETIX 4.05 (Belkhir et al., 1999-2004). Unbiased expected (He) heterozygosity was calculated in GenAlex 6.5 (Peakall & Smouse, 2012). Allelic richness (A_R) was calculated in the function allele.rich in PopGeReport. v3.0.4 (Adamack & Gruber, 2014), in R (R Core Development Team, 2017), which uses the rarefaction approach to deal with differences in sample sizes.

Clustering and genetic differentiation

The fixation index F_ST, based on Weir & Cockerham (1984), was estimated with GENETIX 4.05 (Belkhir et al., 1999-2004) to describe genetic differentiation based on the variance in allele frequencies among locations. Bonferroni corrections (adjusted P = 0.05) were performed by multiple assessment (Rice, 1989). In addition, F_ST were estimated using and not using ENA algorithm as implemented in FreeNA.

A Discriminant Analysis of Principal Components (DAPC) was performed using the ADEGENET package (Jombart, 2008) in R Core Development Team, (2017), through cross-validation of the POPPR package (Jombart, Devillard & Balloux, 2010) with the xvalDapc function. This function chooses the number of axes to retain by testing different numbers of PCs, and the DAPC quality is subsequently evaluated by cross-validation. The DAPC is carried out on the set comprised of 90% of the observations, and then used to predict the groups of the remaining 10% of observations. The number of PCs associated with the lowest Mean Squared Error is then retained for the DAPC. This method combines the advantages of principal component analysis, by assuming neither Hardy-Weinberg equilibrium nor linkage disequilibrium (Jombart, 2008), and from discriminant analysis, by which it attempts to summarize the genetic differentiation between groups, while omitting variation within groups (Jombart, Devillard & Balloux, 2010). Two independent analyses were carried out, one using all 16 locations and the other using only the 12 surface locations, while excluding the four cave locations.

Analysis of molecular variance (AMOVA) was conducted for four different partitions of the 16 locations using ARLEQUIN 3.5 (Excoffier & Lischer, 2010). The following four grouping models were employed under different criteria to maximize variance between the groups and thus reveal the genetic structure contained in the data: (A) All 16 locations; (B) phylogenetic criteria considering two subclades of A. mexicanus [(CC)+(SF, GV, TR, AL, ML, PCH, LCA, SLP, V)], subclades of A. aeneus [(CAT), (TE, TA, RTZ)], and the clade of A. hubbsi (SAB, TIN); (C) hydrographic basin criteria considering the seven independent river basins included in the present study [Bravo (CC), San Fernando (SF), Soto la Marina (GV, TR), Pánuco (AL, ML, PCH, SAB, TIN, LCA, SLP), Túxpan (V), Papaloapan (CAT), Grijalva-Usumacinta (TE, TA, RTZ)]; (D) Localities were grouped according to their degree of genetic differentiation as measured by the F_ST, which resulted in the grouping of those locations that presented non-significant differentiation levels [(CC), (SF), (GV, TR), (AL, LCA, SLP), (V), (PCH), (ML), (SAB, TIN), (CAT), (TE, TA, RTZ)]. Ten thousand permutations were used in each of these analyses.

A Bayesian model-based clustering method was also performed using STRUCTURE 2.3.3 software (Pritchard, Stephens & Donnelly, 2000; Pitchard, 2007) in order to evaluate population structure and assign individuals to genetic clusters. This approach was conducted using the admixture model, which excluded information on the location of origin and assumed independence among loci and non-informative priors. Values were tested for K = 1 to K = 16, including all locations, and for K = 1 to K = 12, without cave locations, that is, including only surface locations. The mean and variance of the log likelihoods for each K were calculated from 10 independent runs of 100,000 iterations in order to determine the highest posterior probability K. To estimate the true number of clusters represented along the sampled range of Astyanax spp. (Evanno, Regnaut & Goudet, 2005), ΔK statistics were calculated by http://taylor0.biology.ucla.edu/struct_harvest/ (Earl, 2011). Individual and mean population membership coefficients of ancestry in our inferred demes were presented graphically in DISTRUCT, version 1.1 (Rosenberg, 2004).

Isolation by distance and historical contemporary gene flow

The pattern of isolation-by-distance (IBD) among locations was assessed by performing a Mantel test, based on pairwise F_ST and geographic distances between distinct locations, using the Isolde algorithm with the on-line version of GENEPOP (http://genepop.curtin.edu.au/). Four distinct Mantel tests were performed as follows: (1) using all 16 locations; (2) using only the surface populations; (3) using only the locations located north of the TMVB and (4) using only the locations south of the TMVB. A total of 10,000 permutations were required to estimate the 95% upper tail probability of the matrix correlation coefficients.

Historical gene flow (for approximately the four previous N_e generations (Beerli & Felsenstein, 2001) was estimated using the program MIGRATE-N 3.6.11 (Beerli & Felsenstein, 2001). Mutation-scaled migration rates (M = m/µ, where m = migration rate and µ= mutation rate) and mutation-scaled size theta (Φ = 4N_eµ, where N_e = effective population size, and µ= mutation rate) were estimated between locations. Different migration model parameter distributions were used between every location detected by MIGRATE. Microsatellite mutation was modeled as a continuous Brownian process and the mutation rate was set to constant for all loci. A static heating scheme was used with four chains and temperatures (1.0, 1.5, 3.0 and 1000000), for which the initial Φ and m values were generated from the F_ST option. Exponential priors were placed for both Φ (bounded between 0 and 50.0) and M (bounded between 0 and 100.0). The analyses were conducted using 10 replicates of a single long Markov chain, with 100,000 steps recorded for every 100 generations and the first 500,000 trees per run discarded as burn-in. Three patterns of migration were estimated in A. mexicanus from the Gulf of Mexico slope: among surface populations, among cave populations, and from surface to cave. Gene flow from cave to surface has been regarded as a nonexistent event, as reflected in its negligible estimate (1.5 ×10⁻⁶ to 10.2 ×10⁻⁴, Fumey et al., 2018).

Recent gene flow (in the past few generations) was estimated using BAYESASS 3.04 (Wilson & Rannala, 2003). This program estimates the posterior probability of an individual’s migratory history and thus allows estimation of the rate and direction of recent dispersal. The Markov Chain Monte Carlo method was run for 10,000,000 iterations. Delta values (i.e., maximum parameter change per iteration) were left as default (Beerli & Felsenstein, 2001).

Population structure

Pairwise F_ST values revealed high, moderate, and low divergences (Table 2). High and significant values were observed among the cave locations (0.77–0.85), except in the comparison between the SAB and TIN caves (0.27). Moderate and significant values were recorded for most of the pairwise comparisons of surface locations, regardless of geographic distance between them (0.11–0.30), while low but significant values were found for the comparisons of surface locations of the Pánuco basin (AL, LCA, SLP) with those of the Papaloapan (CAT) and Grijalva-Usumacinta (TE, TA, RTZ) basins (0.05–0.10). The comparisons made within the Soto la Marina (GV vs TR) and Grijalva-Usumacinta (TE vs TA vs RTZ) basins presented the lowest and non-significant values (0.01–0.02). Pairwise F_ST values using and not using ENA did not present significant differences (Table S6).

An initial DAPC, with all 16 locations, showed the formation of four main clusters (Figs. 2A and 2B) in which 15 discriminant functions were retained. The first two functions accounted for 93.3% of the total variance. Cluster a corresponded to the PCH cave; cluster b included the SAB and TIN caves; cluster c was formed by the ML cave, and cluster d corresponded to all the surface locations. A second DAPC, with only the 12 surface locations, revealed the formation of seven clusters (Figs. 2C and 2D) in which 11 discriminant functions were retained. The first two functions accounted for 94% of the total variance. Cluster e was formed by the CC location; cluster f by the SF location; cluster g included the GV and TR locations; cluster h was formed by the AL, LCA, and SLP locations; cluster i included the V location; cluster j included the CAT location; and finally, cluster k included the remaining surface locations, TE, TA, and RTZ.

The AMOVA analysis indicated that model D of 10 populations showed the highest inter-group variance (F_CT = 0.2449, P = 0.00) and minimal intra-group variance (F_SC = 0.0455, P = 0.00) (Table 3). STRUCTURE also revealed 10 genetic groups. Three cave populations were assigned the greatest probability of ancestry (ML: 99%, PCH: 99%, SAB-TIN: 97%), as were the following epigean populations corresponding to distinct river basins: Bravo (CC: 95%), Soto la Marina (GV+TR: 92%), and Tuxpan (V: 91%). Other populations that also were homogeneous, but with lower assignment probability, were the Papaloapan (CAT: 85%), Grijalva-Usumacinta (TE+TA+RTZ: 84%), and Pánuco (AL+LCA+SLP: 70%) basins (Fig. 1). The Pánuco and Grijalva-Usumacinta basins presented a genetic admixture between them (Fig. 1), which was reflected in both analyses (DAPC and F_ST values) (Fig. 2; Table 2). Another genetic admixture was recorded between the populations of the Soto la Marina and Papaloapan basins (Fig. 1).

In summary, the results of F_ST, DAPC, AMOVA, and STRUCTURE agree that the number of homogeneous genetic groups along the Atlantic slope is ten, consisting of three cave populations, and seven surface populations.

Isolation by distance, effective sample size, and gene flow

In the Mantel tests for both the 16 (all) and 12 (surface only) locations, the results showed a non-significant correlation (R² = 0.07, P = 0.908; R² = 0.067, P = 0.125, respectively). A fourth Mantel test, using all surface locations except for CAT (11 locations), also recorded a significant correlation between genetic and geographic distances (R² = 0.633, P = 0.0; Fig. S1B).

The smallest values of the effective population sizes (Ne) were presented by caves locations (average values from 1,153.334 ML to 3,260 for TIN; Fig. 3). For surface locations, the lowest Ne values were recorded in northern localities (average values of 5,316.667 for CC and 10,816.667 for SF; Fig. 3). Remaining surface locations presented higher Ne values (average values from 16,677.78 for LCA to 22,450.00 for TE; Fig. 3).

In general, the historical migration rates (m) among locations and populations show low values (below m = 0.005), except for the highest genetic flow values recorded among distinct basins. Of particular note is the TR location from the Soto la Marina basin that presented a symmetric flow with the CC population, and the same TR that presented an asymmetric flow toward the cave population PCH. All three of these populations are assigned to A. mexicanus (Fig. 4). Intermediate gene flow values also were recorded among distinct basins (Fig. 4). It should be noted that in the lowest gene flow there were no differences within basins, between basins, or between caves (Fig. 4).

Considerable recent gene flow using BAYESASS 3.04 was detected only unidirectionally among the following locations within the same cave or basin system that formed a single population: for cave locations, from TIN to SAB, and for surface locations, from GV to TR in the Soto la Marina River basin, from SLP and LCA to AL in the Pánuco River basin, and from TA to TE in the Grijalva-Usumacinta River basin (Table 4).

Characterizing population genetic variability

The ten populations varied in terms of the polymorphism of the microsatellite loci. The number of alleles per locus (N_a) ranged from 1 to 23 (average = 8.30), number of effective alleles (N_ae) ranged from 1 to 11.8 (average = 4.0), allelic richness (A_R) ranged from 1 to 13.9 (average = 5.8), observed heterozygosity (Ho) presented values from 0 to 0.923 (average = 0.485), expected heterozygosity (uHe) ranged from 0 to 0.921 (average = 0.556), and the fixation index (F_IS) ranged from −0.292 to 0.665 (average = 0.145). As expected, the cave populations (ML, PCH, SAB/TIN) presented the lowest values of genetic variability and the highest values of the fixation index, while the surface populations showed higher variability and lower values of the fixation index (Fig. S2 and Table S3).

Confirming multiple independent cave invasions

Until a few years ago, the widely accepted evolutionary scenario regarding the origin of A. mexicanus cavefish was that some of these populations were ancient relict populations (i.e., hundreds of thousands of years to millions of years) whose surface relatives were extinct, while other cave populations were considered more recent, and therefore closer to current surface fish populations. This scenario involved two types of cave populations derived from an ancient invasion and a recent invasion (revision in Gross, 2012). However, recently Fumey et al. (2018) found that the origin of the blind cave fish is much more recent (<20,000 years), thus refuting the hypothesis of the two ancient invasions. In another study based on independent approach, the findings also refuted the origin theory of two ancient invasions, although they supported at least four independent origins for the cavefish populations (Coghill et al., 2014). Bradic et al. (2012), while recognizing the origin of the cave populations through the two ancient invasions, also found four independent invasion events within the Sierra del Abra (Pachón, Chica cave, Micos, and the six adjacent central caves [Yerbaniz, Japones, Arroyo, Tinajas, Curva, and Toro]) which correspond to the oldest invasion wave, plus a fifth invasion within the Sierra de Guatemala. The isolation and differentiation pattern of the cave populations included in the present study, ML, PCH, and SAB+TIN caves (Figs. 1; 2A and 2B; Table 2), indicate three independent origins, thus, corroborating those findings pointing to multiple independent subterranean invasions. Greater efforts in field research and modern genomic techniques will surely be required to accurately determine the invasion times of this highly dispersive species into caves.

The case of the SAB+TIN cave fish, while they were not separated by the STRUCTURE and DAPC analyses (Figs. 1; Figs. 2A and 2B), both presented considerable and significant differentiation coefficients (F_ST = 0.27, Table 2). Despite the physical connectivity between these caves, they are not adjacent and are in fact located at different elevations, with the Sabinos cave located at a higher elevation (239 masl) than the Tinajas cave (166 masl) (Elliott, 2016). Such inconsistency in the level of structure (Figs. 1; Figs. 2A and 2B) could be the result of two possible scenarios regarding how the SAB+TIN complex was formed. The first (e1) considers that a cave population within “Sistema de los Sabinos” experienced a partial fragmentation, implying the occurrence of an isolation stage with low connectivity as indicated by the historical gene flow (Fig. 4). The second (e2) scenario consists of independent isolation events in each of these caves derived from a common surface ancestor, which would imply two independent cave invasions in “Sistema de los Sabinos”, with subsequent secondary contact. These scenarios could both reflect the very dynamic karst system in northeastern Mexico, associated with the formation of new caves, and/or incidents of separation or fusion of these systems (Strecker, Hausdorf & Wilkens, 2012) (see below for more discussion).

Dispersal pattern of Astyanax mexicanus throughout the Gulf of Mexico slope

In concordance with previous studies (Bradic et al., 2012; Fumey et al., 2018), findings obtained herein showed very low gene flow in both intra-basin locations and populations, and inter-basin populations. However, there were six particular cases that recorded the highest migration rates, which will be discussed next (Fig. 4). Recent findings based on analysis of the whole genome, support the occurrence of considerable historical and contemporary gene flow between cave and surface populations (Herman et al., 2018).