Population genetics of swamp eel in the Yangtze River: comparative analyses between mitochondrial and microsatellite data provide novel insights

Ethics statement and Sample collection

Procedures involving animals and their care were approved by the Animal Care and Use Committee of Anhui Academy of Agricultural Sciences under approval number 201003076. Field experiments were approved by Fisheries Bureau of Anhui (project number: FB/AH 2017-10).

A total of 180 individuals were collected using net from six sites along the Anhui basin of the Yangtze River (Table S1). Sampling sites from tributaries of the Yangtze River included Dang Tu (DT), Fan Chang (FC) and Huai Ning (HN); sampling sites from the main stem of the Yangtze River included Wu Wei (WW), Gui Chi (GC) and Wang Jiang (WJ) (Fig. 1).

DNA extraction and Marker genotyping

Total genomic DNA was extracted from muscle tissue using a standard phenol/chloroform procedure via proteinase K digestion (Sambrook, Fritsch & Maniatis, 1989), and then kept at −20 °C for PCR amplification.

The mitochondrial cytochrome c oxidase subunit I (COI) gene and cytochrome b (Cyt b) gene were chosen. Two pairs of primers were designed here for the amplification (Table S2). PCR were conducted in 50 µL reaction mixtures containing 200 ng template DNA, 5 µL 10 × buffer (TaKaRa, Dalian, China), 4 µL MgCl₂ (2.5 mol/L), 3 µL dNTP (2.5 m mol/L), 2 µL of each primer (5 µmol/L), and 1 U Taq DNA polymerase (5 U/ µL, TaKaRa). PCR conditions were as follows: initial denaturation (95 °C, 1 min), then 35 cycles of denaturation (94 °C, 50 s), primer annealing (55 °C, 45 s), and elongation (72 °C, 1 min) and a final extension (72 °C, 10 min). All fragments were sequenced in both directions with an ABI3730 automated sequencer (Invitrogen Biotechnology Co., Ltd, USA). Then these two gene sequences were combined for subsequent analyses.

Eight unlinked polymorphic microsatellite loci were selected from previous studies (Table S1) (Lei et al., 2012; Li et al., 2007; Zhuo et al., 2011). PCR were conducted in 50 µL reaction mixtures containing 200 ng template DNA, 5 µL 10 × buffer (TaKaRa, Dalian, China), 4 µL MgCl₂ (2.5 mol/L), 2.5 µL dNTP (2.5 m mol/L), 2 µL of each primer (5 µmol/L), and 1 U Taq DNA polymerase (25 U/ µL, TaKaRa). PCR conditions were as follows: initial denaturation (95 °C, 5 min), then 32 cycles of denaturation (94 °C, 30 s), primer annealing (57 °C, 60 s), and elongation (72 °C, 90 s) and a final extension (72 °C, 5 min). Genotypes were detected by ABIPRISM 3730.

Mitochondrial sequence

Sequences were assembled by DNASTAR Lasergene package. Subsequent homologous alignment was performed by Mafft v.7 online program (https://mafft.cbrc.jp/alignment/software/) (Katoh, Rozewicki & Yamada, 2019).

Haplotype and nucleotide diversity were estimated using DNAsp V.6 (Rozas et al., 2017). A haplotype network was constructed using Median Joining (MJ) in NETWORK v.5.0 (Bandelt, Forster & Röhl, 1999). Analysis of Molecular Variance (AMOVA) was performed using Arlequin v.3.11 (Excoffier, Laval & Schneider, 2005). Genetic variation within and among sampling sites was assessed. Pairwise F_st was estimated in order to evaluate the levels of population differentiation (Slatkin & Barton, 1989) and the P values were corrected using multiple testing.

The demographic history was explored using three approaches, e.g., neutrality tests, mismatch distribution and Bayesian Skyline Plots (BSP) analyses. Tajima’s D (Tajima, 1989) and Fu’s F_s (Fu, 1997) values were calculated using DNAsp V.6.Mismatch distribution analyses were performed using Arlequin v.3.11 (Rogers & Harpending, 1992). The expansion time was calculated by the τ value with the equation τ=2 µt, where µrepresents the nucleotide mutation rate and t represents the estimated expansion time. BSP analysis was performed using Beast v1.10.4 (Suchard et al., 2018) under an uncorrelated relaxed clock mode for 5 × 10⁷ generations.

Microsatellites data

The results of 8 microsatellites loci were read using GeneMarker (Holland & Parson, 2011) and reformatted using Convert v.1.31 (Glaubitz, 2004). Hardy-Weinberg equilibrium (HWE) tests were performed using Popgene v 1.32 (Yeh et al., 1997). Expected and observed heterozygosity were calculated with Arlequin v.3.11 (Rogers & Harpending, 1992). AMOVA was implemented with Arlequin. Pairwise F_st was computed based on Slatkin’s method (Slatkin & Barton, 1989). The geographical and genetic distance between sample sites was measured by GPS and Popgene v 1.32, respectively. The correlation between geographical and genetic distance was analyzed using Pcord v 5 (Grandin, 2006).

Population structure was estimated using an MCMC (Markov Chain Monte Carlo) algorithm as implemented in Structure v.2.3.3 (Hubisz et al., 2009). The number of clusters (K) was calculated under 1 × 10⁶ iterations with 10 replications and the optimal number of K was deduced by Structure Harvester Web v.0.6.94 (Evanno, Regnaut & Goudet, 2005; Earl, 2012).

Mitochondrial genes

A total of 1752 bp of mitochondrial sequence (COI 665 bp, Accession number: MN097948–MN098127; Cyt b 1087 bp, Accession number: MN098128–MN098307) were obtained for analyses. The contents of the bases A, T, G and C were 24.6%, 29.3%, 14.6% and 31.5% respectively, which showed obvious anti-G bias (Saccone et al., 1999).

The 180 mitochondrial sequences corresponded to 86 distinct haplotypes (Table 1). All haplotypes were divided into four clades based on MJ method (Fig. 2). Clade A was the largest one which contained samples from five sampling sites. Haplotype 5 (H-5) had the largest number of shared individuals, and its central placement in the network suggests that this is the ancestral haplotype. The other three clades were separated by 13, 47, 24 mutational steps, respectively. Clade B and C only contained samples from the HN and FC sampling sites, respectively. Clade D mainly consisted of DT samples.

Haplotype diversity ranged from 0.6620 to 0.9793 and nucleotide diversity ranged from 0.0017 to 0.0148 based on mitochondrial sequence (Table 1). The results of AMOVA showed that genetic variation among sampling sites (71.23%, P < 0.001) were much higher than the variation within the sampling sites (28.77%, P < 0.001) (Table 2A). Subsequent F_st values further confirmed this result. Strong gene flow was detected between the main stem sampling sites (F_st = 0.0242 between GC and WW, F_st = 0.0286 between WJ and WW, F_st = 0.0305 between WJ and GC, see Table 3A). And high differentiation was revealed between the tributary sampling sites (F_st = 0.3069 − 0.9431) (Table 3A). Fu’s F_s and Tajima’s D tests of main stem samples were significant (P < 0.01) but negative. No explicit expansion or decline were revealed for the tributary samples and the Fu’s F_s and Tajima’s D values except the Tajima’s D of FC were not significant (Table 4). Mismatch distribution analysis revealed similar results. The values of sum of squares deviations (SSD) for samples from main stem and DT were not significant (P > 0.05), indicating that sudden expansion could not be rejected (Table 4). The BSP analysis suggested the main stem samples had expanded roughly in 0.46 MYA (Fig. 3).

Microsatellite loci

The eight microsatellite loci amplified unambiguous and repeatable products in the size range expected. All loci were in Hardy-Weinberg equilibrium (P > 0.05). High genetic diversity was also supported by microsatellite data. Expected and observed heterozygosity for the six sampling sites were 0.8052 –0.8749 and 0.5958 –0.7917, respectively (Table 1).

Structure results suggested the highest posterior probability for K = 4 (Fig. S1). The ΔK method revealed four potential genetic clusters, aligning with the three tributaries and all the main stem sites together. Samples from main stem showed high levels of genetic admixture (Fig. 4).

AMOVA was performed using the microsatellite data and suggested that genetic variation was mainly within sampling sites (94.65%, P < 0.001), opposite to the results from mitochondrial data (Table 2B). F_st values suggested low levels of differentiation (F_st = 0.0005 − 0.0982) (Table 3B).

The correlation between genetic and geographic matrixes was assessed using Mantel test (Table S3). The results suggested a significant correlation between them (r = 0.8791, P = 0.004) (Fig. 5).

Population genetics of swamp eel

High levels of genetic diversity were found across sampling sites at both mitochondrial and microsatellite markers. The genetic diversity level of this study except FC sample site was higher than previous study in the same basin (Hd = 0.708 and Pi = 0.002 based on mitochondrial D-loop sequences) (Liang et al., 2016). The genetic diversity of FC samples was the lowest (Hd = 0.6620, Pi = 0.0017). The significant differentiation of three tributary sampling sites was revealed by the population genetic analyses (F_st >0.25). The haplotype network and structure results suggested the tributary samples formed three separate clades which contained site-specific lineages. Significant correlation between genetic and geographic distance was detected. Interestingly, strong gene flow was detected among the main stem sampling sites and the expansion of main stem samples was detected.

It is well known that the swamp eel is a burrowing fish whose fins are vestigial or absent (Nelson, Grande & Wilson, 2016). Compared with most fishes, the swimming ability of eel is weak. Thus, we were curious about the reasons for this long-distance gene flow among main stem sampling sites. The flow rate of main stem in Anhui basin range up to 1.0 m/s (Guo & Xia, 2007). Rapid flow provides the dynamic for the migration of swamp eel. The eggs and juvenile fishes can slip downstream to the farther places. However, due to the flat stream gradient and curved channel, the tributary flow becomes slower (Zhang, YT & Jiang, 2008) and long-distance migration is difficult for swamp eel. We propose that geographic isolation and the habitat patchiness caused by seasonal cutoff are the reasons for the differentiation between tributary samples. The tributary of FC site was much more isolated from the main stem. It connected to the main stem during the wet season and isolated during the dry season. Long-term isolation from the main stem may cause the low genetic diversity of FC samples.

Comparative analyses between mitochondrial and microsatellite data

Analyses of nuclear and mitochondrial markers revealed discordant population structure. Based on mitochondrial data, genetic variation was mainly found among sampling sites. Samples between the tributary sites were highly differentiated (F_st > 0.25) and represented three monophyletic clades. However, microsatellite analyses suggest that the majority of genetic variation is within these sampling sites; samples between the tributary sites were moderately differentiated (0.05 < F_st < 0.15). Mean F_st values among six sampling sites based on mitochondrial and microsatellite data were 0.548 and 0.055, respectively. The mean mitochondrial F_st value (0.548) was almost ten times higher than the F_st (0.055) estimated with microsatellite data.

Our study provided an interesting pattern of discordance between markers for population genetics. According to previous studies, sex-biased dispersal, genetic admixture and lineage sorting may be the potential reasons for the discordance caused by different molecular markers (Funk & Omland, 2003; Qu et al., 2012; Yang et al., 2016; Zarza, Reynoso & Emerson, 2011). Considering the sex reversal in this species, we inferred the unique life history of the swamp eel contributed to this discordance. Initially, the swamp eel is female and provides both mitochondrial and nuclear DNA to the population genetic pool. After spawning, the swamp eel becomes male. Male swamp eels are much bigger and stronger than females. Males could migrate farther and have a higher survival rates, which provides a potential of male sex-biased dispersal. Due to mitochondrial maternal inheritance, male swamp eels only provide nuclear DNA to the population genetic pool (Fig. 6). So male sex-biased dispersal may cause the differences in population structure between the markers. As mentioned above, male stage is much longer than female stage in the whole life of swamp eel that may cause different genetic frequencies between mitochondrial and nuclear data. The ten–fold difference between mitochondrial and nuclear F_st values also confirmed this hypothesis.