Genetic Di ﬀ erentiation of an Endangered Megalobrama terminalis Population in the Heilong River within the Genus Megalobrama

: Megalobrama terminalis , which inhabits the Sino-Russian Heilong-Amur River Basin, has decreased critically since the 1960s. It has been listed in the Red Book of Endangered Fish Species by the Russian Federation in 2004. To guide the utilization and conservation programs of M. terminalis in the Heilong River (MTH), 3.1 kb of mitochondrial DNA (mtDNA) concatenated sequences and sequence-related ampliﬁed polymorphism (SRAP) markers (15 primer combinations) were applied to explore the genetic divergence and population di ﬀ erentiation of MTH within the genus Megalobrama . Clear genetic divergence between MTH and six other populations of the genus Megalobrama was found by haplotype network (mtDNA) and principal component (SRAP) analyses. Moreover, the STRUCTURE analysis based on SRAP data showed that MTH could be assigned to a particular cluster, whereas conspeciﬁc M. terminalis in the Qiantang River and Jinsha River Reservoir belonged to the same cluster. Analysis of molecular variance (AMOVA) and F st statistics for the mtDNA and SRAP data revealed signiﬁcant genetic variance and di ﬀ erentiation among all detected populations. Taken together, the results suggest that MTH has a strong genetic di ﬀ erentiation from other populations within the genus Megalobrama , which contributes to e ﬀ ective utilization in artiﬁcial cultivation and breeding of MTH. Furthermore, these results also provide a scientiﬁc basis for the management of MTH as a separate conservation unit.


Introduction
The genus Megalobrama is a member of Cultrinae in Cyprinidae fish. This genus is a kind of middle-large, economically valuable fish and is famous for its tasty meat quality [1]. The genus Megalobrama in China is mainly composed of M. amblycephala, M. terminalis, M. pellegrini and M. hoffmanni [2,3]. Among these species, M. amblycephala and M. terminalis formed a certain aquaculture scale. The M. amblycephala has two new selective breeding varieties (Pujiang No. 1 and Huahai No. 1), and its production reached 700,000 tons in 2014 [4]. In contrast, no new breed of M. terminalis has been systematically selected, and germplasm resources for aquaculture originated from the Qiantang River in Zhejiang Province and the Jinsha River Reservoir of Hongan County in Hubei Province [5,6]. According to statistics from 2016, the aquaculture area of M. terminalis in the Qiantang River (MTQ) has exceeded 600 ha in the areas of Hangzhou, Jinhua and other places in Zhejiang Province [1]. Since the success of artificial reproduction in 1988, M. terminalis in the Jinsha River Reservoir (MTJ) has been extended to the Hubei, Henan, Anhui, Sichuan, Guangdong and Fujian Provinces [6].
Remarkably, there is another native wild population of M. terminalis in the Heilong River (the Amur River in Russia), which is the northernmost, highest-latitude river in China. Although M. terminalis was historically widespread in bodies of water in the Heilong River Basin, such as the Heilong River, the Nen River, the Songhua River, Jingpo Lake and Khanka Lake [7] (Supplementary Materials Figure S1), its habitats have been shrinking year by year since the 1960s with the effects of overfishing and environmental changes [8]. Currently, only dozens of individuals can be captured per year in the Fuyuan section (134 • 28 E, 48 • 37 N) of the Heilong River ( Figure S1). M. terminalis in the Amur River has been included in the Red Book of endangered species of the Russian Federation [9]. M. terminalis in the Heilong River (MTH) is highly popular among local consumers because of its large size and delicious taste. The potential value of these excellent properties determines that MTH will be an ideal commercially exploited species [3,7]. Due to the difficulty of wild fish collection and the high mortality in the course of temporary rearing and transport, however, systematic research on population genetics for MTH is highly limited. Recently, the phylogenetic relationships among seven populations from four species in the genus Megalobrama, including M. terminalis in the Heilong River (MTH), were clarified at the mitochondrial genome level [10]. Moreover, the artificial reproduction of MTH was achieved under pond domestication conditions [11]. However, the genetic divergence and differentiation at a population level between MTH and other populations from M. terminalis and other species within the genus Megalobrama have been not determined, which limits the effective utilization and protection of MTH.
Regarding fish species of the genus Megalobrama, genetic variance and differentiation have been well-evaluated in different natural and selected populations from M. amblycephala using morphometrics, isozymes, mitochondrial DNA (mtDNA), random amplified polymorphic DNA (RAPD) and sequence-related amplified polymorphism (SRAP) markers [12][13][14][15]. However, existing data only reflect the morphological differences between the Qiantang River and the Jinsha River reservoir, as well as between the Heilong River and the Liangzi Lake, in M. terminalis [8,16]. At the interspecies level, Li et al. completed a preliminary comparison of genetic variations among M. terminalis in the Qiantang River, M. amblycephala and M. hoffmanni using morphometrics, isozymes and RAPD markers [17]. It can be seen from existing reports that the investigations on intra-and interspecies relationships involved in M. terminalis are highly limited.
The amplified target site of SRAP is in the coding regions of the genes, and it is a highly reproducible and informative marker system. SRAP markers have been widely used to estimate the genetic diversity and structure, and varietal identification in plants [18,19] and their discriminating power is stronger than those of simple sequence repeats (SSR), inter-simple sequence repeats (ISSR) and RAPD [20]. In addition, this technique has been extended to studies of population genetics in other living organisms [21], including aquatic animals [15,[22][23][24][25][26][27][28]. To understand the genetic relationships, especially the level of genetic differentiation of MTH with MTQ, MTJ and other representative geographic samples within the genus Megalobrama, containing M. amblycephala from the Liangzi Lake (MAL) and the Yi River (MAY) and M. Pellegrini from the Longxi River (MPL) and M. hoffmanni from the Xi River (MHX), we performed this survey based on mtDNA and SRAP markers.

Sampling and DNA Extraction
Two hundred and ten individuals, including 7 wild populations (30 per population) of 4 species in the genus Megalobrama, were sampled by means of drift-net fishing. Sampling sites are shown in Figure 1. The numbers of individuals used for mitochondrial and SRAP analysis are shown in Table S1. Genomic DNA from fin tissues was isolated using the phenol/chloroform method and preserved in 1× TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0). The quality and quantity of the DNA were detected by agarose gel electrophoresis and absorbance at 260 nm and 280 nm. DNA samples were quantified 200 ng/µL and diluted to 50 ng/µL for PCR. All animal experiments were conducted in accordance with the guidelines and approval of the Animal Research and Ethics Committees of Heilongjiang River Fisheries Research Institute.
All sequences amplified by PCR were edited and aligned using the CLUSTAL X program, as implemented in MEGA 7 (Temple University, Philadelphia, PA, USA) [29] and, subsequently, aligned manually. Analysis of population variation, such as the number of haplotypes, haplotype diversity (h) and nucleotide diversity (π), were calculated by DNAsp5.0 software (Universitat de Barcelona, Barcelona, Spain) [30]. A haplotype network was constructed using the median-joining method in NETWORK 5.0.1.1 software (Fluxus Technology Ltd., Clare, Suffolk, UK) to analyze and visualize the relationships among DNA sequences. The software ARLEQUIN v3.5.2 (University of Berne, Berne, Switzerland) [31] was used to complete the analysis of molecular (AMOVA) for detecting genetic variation levels among populations and calculate pairwise Fst values for estimating the extent of genetic differentiation among populations. Significance was tested using 1000 permutations. All accession numbers of the mtDNA fragments are MN958398-MN958523.
All sequences amplified by PCR were edited and aligned using the CLUSTAL X program, as implemented in MEGA 7 (Temple University, Philadelphia, PA, USA) [29] and, subsequently, aligned manually. Analysis of population variation, such as the number of haplotypes, haplotype diversity (h) and nucleotide diversity (π), were calculated by DNAsp5.0 software (Universitat de Barcelona, Barcelona, Spain) [30]. A haplotype network was constructed using the median-joining method in NETWORK 5.0.1.1 software (Fluxus Technology Ltd., Clare, Suffolk, UK) to analyze and visualize the relationships among DNA sequences. The software ARLEQUIN v3.5.2 (University of Berne, Berne, Switzerland) [31] was used to complete the analysis of molecular (AMOVA) for detecting genetic variation levels among populations and calculate pairwise Fst values for estimating the extent of genetic differentiation among populations. Significance was tested using 1000 permutations. All accession numbers of the mtDNA fragments are MN958398-MN958523.

SRAP Primer, PCR Amplification and Analysis
The primers required for SRAP refer to the primer sequences published by Li et al. [21]; 15 primer combinations (Table S3) with clear amplification bands, strong reaction stability and high reproducibility were identified from 240 combinations. All 210 individuals from 7 populations (30 per population) were genotyped using 15 primer combinations selected. The reaction system of 20 µL contained 7-µL H 2 O, 1-µL 10 × PCR buffer, 0.8-µL 1.5-mmol/L MgCl 2 , 0.2-µL 10-mmol/L dNTPs, 0.5 µL of each primer (10 µM) and 0.5 µL of 50-ng/µL genomic DNA. The PCR conditions were as follows: pre-denaturation for 5 min at 94 • C followed by the first cycle for 5 cycles of 94 • C for 1 min, 35 • C for 1 min and 72 • C for 1 min, followed by the second cycle for 35 cycles of 94 • C for 1 min, 50 • C for 1 min and 72 • C for 5 min, followed by an extension of 10 min at 72 • C. The PCR product was detected by 1.5% agarose gel electrophoresis and was further loaded to a 10% polyacrylamide gel electrophoresis with reference to the pBR 322/Msp I marker (TIANGEN, Beijing, China). The silver staining method [32] was used to display the bands. The electropherogram was analyzed by Gel-pro analyzer 4.5 software (Media Cybernetics, Inc., Rockville, MD, USA). In most of the gels, in addition to the markers to indicate the size of the bands, we also clicked a sample from a certain population used for the marker or control between gels. Moreover, the samples were from several populations in the same gel, which were contrasted with each other to reduce the influence of differences between gels on the result judgment ( Figure S2). The principle for determining bands was the stability of amplification. The judgment of weak bands needed to take into account the overall (multiple populations) amplification situations. If some bands were weak but others were clear in the same pair of primers, it should be caused by the different amplification efficiency among samples; such bands were recorded in the statistics. Whereas bands in different populations were weak, indicating unstable amplification; such bands were not included. Undecidable bands were regarded as missing data. Different and reproducible SRAP bands were scored as present (1) or absent (0), and ambiguous bands or missing data were recorded with number "9".
The polymorphism information content (PIC) of each SRAP marker was calculated by the following formula: , where f is the frequency at which the band appeared in all samples [33]. Principal component analysis (PCoA) was performed based on genetic distance using the GenAIEx 6.5 package (Australian National University, Canberra, Australia) [34] to examine the positional relationship of each group on the two-dimensional clustering map of the main coordinates. Genetic structure analysis was carried out based on the Bayesian model using Structure 2.3.4 software (Stanford University, Palo Alto, CA, USA) [35,36]. The number of subgroups was set to 1~7 with 10 repetitions. The length of the burn period for each run was set to 10,000, and the number of iterations for Markov chain Monte Carlo (MCMC) after burning was set to 100,000. The results of the calculations were analyzed online using Structure Harvester program (University of California, Santa Cruz, CA, USA), and the most likely K values were determined based on LnP(K) or ∆K values. The analysis of molecular variance (AMOVA) and the pairwise Fst values with 1000 permutations were performed using ARLEQUIN v3.5.2 software [31].

Genetic Variability
The effective sequences of the D-loop, ND2 and CYTB from each specimen were aligned and processed separately. The concatenated sequence length for haplotype analysis was 3116 bp in samples from M. terminalis and M. pellegrini, including a D-loop of 930 bp, ND2 of 1045 bp and CYTB of 1141 bp (Table S2) Table 1. Fifteen SRAP primer combinations generated 221 stable and reproducible bands, of which 214 (97%) were polymorphic ( Table 2). The amplification profiles of primer combinations 1f5r, 2f9r, 4f20r and 7f20r were shown in Figure S2. The most polymorphism was from the 7f20r combination, which produced 22 polymorphic bands. The PIC of each combination ranged from 0.262 to 0.412, with an average of 0.311 (Table 2). Based on the detection values of four SRAP diversity parameters, seven populations can be divided into two distinct groups in which data for different parameters were similar within each group. Three  Table 3).

Population Genetic Structure and Differentiation
Median-joining network analysis showed that MTH formed a single haplogroup, whereas there were only three haplotypes in MTJ, and one of them was clustered with the haplogroup of MTQ. Two populations for MAL and MAY constituted a haplogroup in which there were two shared haplotypes ( Figure 2). Diversity 2020, 12, x 6 of 12

Population Genetic Structure and Differentiation
Median-joining network analysis showed that MTH formed a single haplogroup, whereas there were only three haplotypes in MTJ, and one of them was clustered with the haplogroup of MTQ. Two populations for MAL and MAY constituted a haplogroup in which there were two shared haplotypes ( Figure 2).  In the principal component analysis (PCoA), the first and second axes explained 30.79% and 47.80% of the variation, respectively. The single scattered point clusters corresponding to MTH, MPL and MHX had almost no overlap with other populations, while the points corresponding to MTQ and MTJ and between MAL and MAY had a certain degree of overlap ( Figure 3).
In the principal component analysis (PCoA), the first and second axes explained 30.79% and 47.80% of the variation, respectively. The single scattered point clusters corresponding to MTH, MPL and MHX had almost no overlap with other populations, while the points corresponding to MTQ and MTJ and between MAL and MAY had a certain degree of overlap ( Figure 3).  To characterize the genetic variance among populations in four species of the genus Megalobrama, four groups were considered based on each species and compared on the basis of hierarchical AMOVA with both mtDNA and SRAP markers (Table 4). Significant genetic differentiation was detected among the four groups (p < 0.01), among the seven populations within the groups (p < 0.01) and within the populations (p < 0.01). In addition, each of the seven populations showed highly significant genetic differentiation (p < 0.01), except between MAL and MAY (p < 0.05), by pairwise Fst based on mtDNA and SRAP genotypes ( Figure 5).   To characterize the genetic variance among populations in four species of the genus Megalobrama, four groups were considered based on each species and compared on the basis of hierarchical AMOVA with both mtDNA and SRAP markers (Table 4). Significant genetic differentiation was detected among the four groups (p < 0.01), among the seven populations within the groups (p < 0.01) and within the populations (p < 0.01). In addition, each of the seven populations showed highly significant genetic differentiation (p < 0.01), except between MAL and MAY (p < 0.05), by pairwise Fst based on mtDNA and SRAP genotypes ( Figure 5). To characterize the genetic variance among populations in four species of the genus Megalobrama, four groups were considered based on each species and compared on the basis of hierarchical AMOVA with both mtDNA and SRAP markers (Table 4). Significant genetic differentiation was detected among the four groups (p < 0.01), among the seven populations within the groups (p < 0.01) and within the populations (p < 0.01). In addition, each of the seven populations showed highly significant genetic differentiation (p < 0.01), except between MAL and MAY (p < 0.05), by pairwise Fst based on mtDNA and SRAP genotypes ( Figure 5).

Discussion
Our recent research clarified the characterization of the mitochondrial genome of MTH in the genus Megalobrama [10]. Based on the analysis of the genomic-level variation sites, the concatenated fragment of D-loop, ND2 and CYTB was selected for the present study. MTH displayed a low level of genetic variability in the mtDNA sequences (π = 0.00045) compared with the other two M. terminalis populations (MTQ, π = 0. 0033 and MTJ, π = 0.0036), while MTH was at the same level as MPL (π = 0.00073), MAL (π = 0.00023) and MAY (π = 0.0002). This level of variation of MPL was highly consistent with a previous report (π = 0.00077) based on the D-loop fragment [38], providing evidence of the reliability of this data. MTH need at least five years to reach sexual maturity, and it is easy to catch on account of body shape [7,9]. These characteristics of MTH makes it more susceptible to environmental and anthropogenic factors; we therefore speculate that the low diversity of MTH may be due to the bottleneck effect. Since samples were not collected before the population declined, the exact mechanism needs further study.
Moreover, SRAP markers have also been used in genetic variability analyses at the nuclear genome level. Among these markers in the MAL, we obtained 82 polymorphic bands, which is similar to the result (88 polymorphic bands) reported by Ji et al. based on 13 pairs of primers in M.

Discussion
Our recent research clarified the characterization of the mitochondrial genome of MTH in the genus Megalobrama [10]. Based on the analysis of the genomic-level variation sites, the concatenated fragment of D-loop, ND2 and CYTB was selected for the present study. MTH displayed a low level of genetic variability in the mtDNA sequences (π = 0.00045) compared with the other two M. terminalis populations (MTQ, π = 0. 0033 and MTJ, π = 0.0036), while MTH was at the same level as MPL (π = 0.00073), MAL (π = 0.00023) and MAY (π = 0.0002). This level of variation of MPL was highly consistent with a previous report (π = 0.00077) based on the D-loop fragment [38], providing evidence of the reliability of this data. MTH need at least five years to reach sexual maturity, and it is easy to catch on account of body shape [7,9]. These characteristics of MTH makes it more susceptible to environmental and anthropogenic factors; we therefore speculate that the low diversity of MTH may be due to the bottleneck effect. Since samples were not collected before the population declined, the exact mechanism needs further study.
Moreover, SRAP markers have also been used in genetic variability analyses at the nuclear genome level. Among these markers in the MAL, we obtained 82 polymorphic bands, which is similar to the result (88 polymorphic bands) reported by Ji et al. based on 13 pairs of primers in M. amblycephala of the Liangzi Lake [15], indicating that SRAP has desirable repeatability and stability. On the other hand, based on results of the SRAP analysis, the genetic variability level of MTH was comparable to MTQ and MTJ, which was inconsistent with the low variation level of MTH obtained by the mtDNA fragments. In many species, including fish, discordant results between nuclear diversity and mtDNA diversity are common [39][40][41][42], which may be related to inheritance patterns, mutation rates and differences in marker characteristics between the mitochondrial and nuclear genomes. Another probable reason is that the effective population size of MTH is sufficiently small, so that the combined effects of genetic drift produced a greater impact on mtDNA diversity than on nuclear DNA diversity [43].
Although phylogenetic research clarified that MTH shared a common ancestor with MTQ, MTJ and MPL, and provided evidence that MPL should be regarded as one of the geographical populations of M. terminalis, the exact genetic differentiation between MTH and other populations in the genus Megalobrama remains to be determined at the population level. In this study, information from the haplotype (D-loop, ND2 and CYTB) network and PCoA (SRAP markers) analysis suggested that there was a strong genetic divergence of MTH with other closely related populations. Furthermore, the STRUCTURE analysis showed that MTH could be assigned to a particular cluster, whereas MTQ and MTJ belonged to the same cluster. These data demonstrated that MTH underwent long-term reproductive isolation resulting from habitat fragmentation and barriers to gene flow. On the other hand, this finding also reflected that MTH maintains a pure gene pool without any gene pollution from farmed M. terminalis or M. amblycephala. Notably, the hierarchical AMOVA and pairwise Fst values results for the mtDNA and SRAP data revealed significant genetic variation and differentiation among all populations, and it showed higher pairwise Fst values between MTH and MPL than MTH with MTQ and MTJ. These findings indicate the genetic differentiation level of MTH from other populations, and further provide evidence to support the species taxonomy of MTH [10]. It is not difficult to conclude that MTH has generated a substantial population structure.
The Heilong River Basin is currently one of the most well-preserved bodies of waters for fish diversity in China [44]. Fish species distributed in the Heilong River usually have the characteristics of large size, strong stress resistance and good meat quality [45]. As the only naturally distributed fish species in the genus Megalobrama north of the Yellow River, MTH possesses notably high economic, scientific and ecological value. In particular, MTH has the potential to be an excellent hybrid parent, because long-term segregation in unique habitats may enable it to evolve key genes required for the breeding of fish species in the genus Megalobrama. To the best of our knowledge, this report is the first to describe the genetic divergence and population differentiation between MTH and other populations in the genus Megalobrama. These data provide an important scientific basis for the further artificial cultivation, breeding application and scientific protection of this endangered species in the future.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/10/404/s1: Figure S1: Historical distribution of M. terminalis in the Heilong River Basin. HLR, the Heilong River (Fuyuan); AR, the Amur River; SHR, the Songhua River; NR, the Nen River; JPL, Jingpo Lake and KL, Khanka Lake. A red dot represents that wild M. terminalis in the location has disappeared. A blue dot indicates that wild M. terminalis in the location is still present. Figure Table output following the method of Evanno et al. [37]. Yellow highlight indicates the largest value in the Delta K column. Table S1: Sampling sites and number of individuals used for mtDNA and SRAP analysis. Table S2: Information on fragments for amplification and haplotype analysis in the D-loop, ND2 and CYTB genes. Reference mtDNA sequence of Megalobrama terminalis in the Heilong River (MTH) (accession no. MH289765).