The swamp eel (Monopterus albus) is a typical sex reversal fish, belonging to the family Synbranchidae, which usually inhabits swamps, ponds and rice fields (Nelson, Grande & Wilson, 2016). Due to its high nutritional value and good taste, the swamp eel is used as a significant aquatic food in China. In 2018, 0.32 million tons of swamp eel were produced by Chinese aquaculture (data from The Ministry of Agriculture of China, 2019).
With the development of swamp eel farming, population genetic research of the wild eel became a hot topic, which has been used to guide the genetic breeding and swamp eel aquaculture. Several studies have used different markers including microsatellites, mitochondrial sequences and Inter-Simple Sequence Repeats (ISSR) to explore swamp eel population genetics (Lei et al., 2012; Li et al., 2013; Liang et al., 2016). Previous studies suggested a rapid decrease of wild swamp eel caused by the large use of pesticides and over-exploitation (Li et al., 2013; Liang et al., 2016). Cultured populations were genetically less diverse than wild populations. Due to farm escapes, the populations of wild swamp eel suffer from genetic homogenization and degeneration of genetic characters (Li et al., 2013). Thus, further research of population genetics and dynamics are necessary to protect the wild swamp eel.
Monopterus albus is a sex reversal fish with a unique life history (Liu, 1944; Mazzoni et al., 2018; Qu, 2018). It starts the reproductive cycle as a functional female. During 1 to 1.5 years of age, females obtain well-developed ova. After spawning, the swamp eel’s sex is reversed from female to male. The intersexes appear in the two-year-old age class. Then it lives as male without sex reversal (Liem, 1963). Generally, swamp eels can live to 6-8 years. Mitochondrial DNA is maternally inherited. For this reason, the inheritance patterns of protogynous hermaphrodites differ to that of species that do not practice sex reversal (Coscia et al., 2016). Therefore, it may be important to include both nuclear and mitochondrial markers to explore the population genetics of sex reversal fish.
This fine-scale study was performed across six sites along the Yangtze River including three sites on the main stem and three sites from tributaries. Subsequently, the genetic structure of these six sampling sites was explored, using sequence from two mitochondrial genes and eight microsatellite loci. Here, we show how different types of molecular markers can provide new insights regarding the population genetics of sex reversal fish.
Materials & Methods
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 MgCl2 (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 MgCl2 (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.
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 Fst 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 Fs (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 × 107 generations.
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 Fst 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 × 106 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).
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.
|Sampling sites||Mitochondrial genes||Microsatellite loci|
|Num. of haplotypes||Hd||Pi||Num. of alleles||He||Ho|
|WW (Main stem)||24||0.9793||0.0078||14||0.8638||0.7292|
|GC (Main stem)||15||0.8480||0.0021||13||0.8420||0.7417|
|WJ (Main stem)||19||0.9240||0.0023||14||0.8516||0.7802|
Notes. Hd represents haplotype diversity; Pi represents nucleotide diversity; He represents expected heterozygosity; Ho represents observed heterozygosity.
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 Fst values further confirmed this result. Strong gene flow was detected between the main stem sampling sites (Fst = 0.0242 between GC and WW, Fst = 0.0286 between WJ and WW, Fst = 0.0305 between WJ and GC, see Table 3A). And high differentiation was revealed between the tributary sampling sites (Fst = 0.3069 − 0.9431) (Table 3A). Fu’s Fs 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 Fs 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).
|Table 2a AMOVA of 6 sampling sites indicating the source of variation by mitochondrial data used.|
|Source of variation||df||Sum of squares||Percentage of variation||P value|
|Among sampling sites||5||1965.433||71.23||<0.001|
|Within sampling sites||174||908.533||28.77||<0.001|
|Table 2b AMOVA of 6 sampling sites indicating the source of variation by microsatellite loci used.|
|Source of variation||df||Sum of squares||Percentage of variation||P value|
|Among sampling sites||5||76.325||5.35||<0.001|
|Within sampling sites||352||672.105||94.65||<0.001|
|Table 3a Pairwise values of Fst (below diagonal) and P (above diagonal) between sampling sites estimated using mitochondrial data.|
|Table 3b Pairwise values of Fst (below diagonal) and P (above diagonal) between sampling sites estimated using 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). Fst values suggested low levels of differentiation (Fst = 0.0005 − 0.0982) (Table 3B).
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 (Fst >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 (Fst > 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 < Fst < 0.15). Mean Fst values among six sampling sites based on mitochondrial and microsatellite data were 0.548 and 0.055, respectively. The mean mitochondrial Fst value (0.548) was almost ten times higher than the Fst (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 Fst values also confirmed this hypothesis.
Our study used two data sets, mitochondrial DNA and microsatellites, to explore the demography, genetic variation and population structure of swamp eels. Compared with previous studies, high levels of genetic diversity suggest that swamp eels are an abundant resource in the Anhui basin and have potential commercial value. Samples from each tributary site in this study should be treated as an independent genetic unit. The unique sex reversal life history of the swamp eel may be a significant factor affecting the population genetic structure and may generate the discordance we found between different molecular markers. Our study provides novel insights regarding the population genetics of sex reversal fish.