The Yangtze Diluting Water and the Ocean Current System Caused an Obviously Genetic Divergence of Meretrix meretrix in China

Population genetic structure of marine organisms has been partially explained by biophysical factors, such as the biology, ecology and behavior of the species as well as hydrographical barriers to dispersal. Abundant and complex ocean currents system was served as the barriers or the accelerators for gene flow because the marine environment was lacked obvious space limitations. The gene flow was enhanced or reduced the genetic divergence of different local samples with a combination of other factors (Gilg and Hilbish, 2003; Weersing and Toonen, 2007). The life histories of marine bivalve contain a pelagic larval duration (PLD). The sessile adults would product vast planktonic larvae which were passively dispersed in the ocean by currents. Ocean currents and the river diluting water were the important factors that could influence the dispersal of larvae (White et al., 2010; Panhwar et al., 2018). The pelagic larval duration and ocean currents would determine the length of diffusion distance, and affect the gene flow between different populations (Weersing and Toonen, 2007). Without swimming organs, the larvae had to spread by following the currents until juvenile stage. In marginal sea of China, several coastal current and the river diluting water form a complex and diverse system. The coastal current will change direction seasonally. Therefore, the larvae will spread to different places seasonally and change the gene exchange between different populations which result gene flow. However, the coastal current was not the single factor that affects the gene flow of bivalve. With the abundant river system in the coastline of China, a great deal of fresh water flows into the ocean every year and cut the coastal current. As the third longest river in the world, the Yangtze River injects huge fresh water into ocean every year in summer. It was a significant factor that would influence the larvae disperse in the ocean (Su and Yuan, 2005; Ni et al., 2015).

The hard clam, Meretrix meretrix L., belongs to the Meretrix, with a widely distribution in the coastline and estuarine areas of Indo-Pacific region (Tang et al., 2006). M. meretrix is a commonly species and seafood in these areas. It is also an economically species for local people to develop aquaculture and processing industries (Zhang et al., 2005). The larva of M. meretrix is adapted to large scale of salinity and temperature so that became a popular species for farming (Lin and Xu, 1997). Since 1990s, due to reduce dependence on wild seeding and the increased demand for aquaculture seeding, artificial breeding became the main source of seeding for M. meretrix in China (Lin et al. 2002). The larval stage of M. meretrix is six days (Shao et al., 2017). With the shorten planktic larval stage and the sessile living of adults, M. meretrix is the better biology material to analysis the genetic structure under the influence of ocean current and river diluting water.

In this study, we focused on the coastal water of China to analysis the genetic variation of the M. meretrix. To investigate the impacts of the Yangtze diluting water and other coastal ocean currents on genetic variance of intertidal species with limited larval dispersal capability, two mitochondrial DNA markers (COI and Cytb) were employed to analysis seven location samples of M. meretrix from north to south in the coast of China. Results from this study may provide some suggests for fisheries managements policymaking to using the germplasm resources and maintain genetic resources.

We collected 134 specimens of Meretrix meretrix from seven locations, Qingdao (QD), Yancheng (YC), Nantong (NT), Zhangzhou (ZZ), Zhuhai (ZH), Beihai (BH) and Haikou (HK). Location details of each samples of M. meretrix were shown in Table I and Figure 1. To avoid confusion specimens with other possibly coexisting clam species, all the samples were identified by experienced fisheries researchers. The adductor muscle was taken out and saved in 100% alcohol before the DNA extraction.

The improved salting out method (Aljanabi and Martinez, 1997) was utilized to extract DNA. DNA quality was estimated by using 1.5% agarose gels electrophoresis and NanoDrop 2000 Spectrophotometer (Thermo Scientific) was used to evaluate DNA concentration in absorption at 260/280 nm. Then, DNA was diluted to a final concentration of 60-80 ng/μl in TE buffer and stored at –20°C for PCR.

N, number of individual; n, number of haplotype; h, haplotype diversity; π, nucleotide diversity.

The complete mitochondrial genome data of M. meretrix (NC_013188) (He et al., 2011) were utilized for the gene sequences search and primers design of COI and Cytb using the software Primer-premer v6.0 (Singh et al., 1998). The primers were as followed: COI-F: 5’-TTATGATAGAACAAGTAAACG-3’ and COI-R: 5’-AGAATAGCATAAATCATAGGC-3’; Cytb-F: 5’-TCTCTATAAGCCATTCTACAAAAAG-3’ and Cytb-R: 5’-TTAATGAGTCAAGATGTAGT-3’.

Each polymerase chain reaction (PCR) was performed in 50 μl volumes, containing 40-70 ng of template DNA, 6 ρmol (Cytb: 8 ρmol) of each forward and reverse primers, 25μl CW0716 2×Taq MasterMix (Cwbiotech Co. Ltd., Peking, China). Applied Biosystems Veriti 96-Well Thermal Cycler (Applied Biosystems, Inc., Foster City, CA, USA) was employed to perform the PCR amplification, consisting of initial denaturation at 95°C for 3 min, the course of reaction was 35 cycles at 95°C for 30s, 54°C (Cytb: 48°C) for 30s, 72°C for 1min, and final elongation at 72°C for 10 min. The products were placed in the electrophoresis on 2% agarose gels. Then, all of them were sent to Bi-directional sequence (Hangzhou TSINGKE Biotechnology Co. Ltd., Hangzhou, China).

The DNA fragments were assembled by the software BioEdit v7.2.5 (Hall, 1999). Each fragment was BLAST with database of the mitochondrial DNA sequence of M. meretrix online (https://www.ncbi.nlm.nih.gov/) to identify the specie of samples and the fragment position. Then, all fragments were aligned by the program ClustalX v2.0 (Larkin et al., 2007). The program jModel Test was utilized to estimate the nucleotide substitution model of these genes (Darriba et al., 2012). The genetic parameters and haplotype discrimination were identified by the software DnaSP v6.0 (Rozas et al., 2017). To further analysis, seven samples were divided into two groups based on geographic position. Group 1 (North) was QD, YC and NT, group 2 (South) was ZZ, ZH, BH and HK. The program Arlequin v3.5 (Excoffier and Lischer, 2010) was used to estimate the pairwise FST values, performed the AMOVA and test the neutrality test based on D value in Tajima (1989) and Fs value in Fu (1997). The haplotype networks were constructed by the Network v5.1 (Bandelt et al., 1999). The program MrBayes was used to construct the Bayesian tree of haplotypes (Huelsenbeck and Ronquist, 2001). The software BEAST v2.0 was used to estimate the population evolution history (Bouckaert et al., 2014). The program Tracer v1.6 was used to plot the Bayesian skyline plot (Rambaut and Drummond, 2007).

672 bp of the COI gene sequence (partly) was analyzed in 122 individuals, and 33 unique haplotypes were detected (GenBank accession numbers: MG698839-MG698871). The optimal nucleotide model was GTR model. The sample size was ranged from 13 (ZZ) to 22 (NT). The number of haplotypes ranged from 7 (ZZ and HK) to 12 (NT). The haplotype diversity ranged from 0.909 (HK) to 0.948 (NT), averaged 0.929 and totaled 0.968. The nucleotide diversity was from 0.0033 (QD) to 0.0110 (ZH), averaged 0.929 and totaled 0.0170 (Table I).

685 bp partial Cytb gene sequence was obtained and aligned with 134 individuals, which detected 33 unique haplotypes (GenBank accession numbers: MG698872-MG698904). The optimal nucleotide model was GTR model. The sample size ranged from 14 (ZH) to 23 (NT). Haplotype numbers ranged from 7 (ZH) to 14 (ZZ). The haplotype diversity varied from 0.879 (ZH) to 0.925 (QD), averaged 0.909 and total 0.959. The nuclear diversity from 0.0048 (NT) to 0.0105 (HK), averaged 0.0072 and total 0.0129 (Table I).

In COI gene, the pairwise FST values ranged from -0.04 to 0.85 (Fig. 2). The result showed that the genetic differentiation between north and south groups was significant. Samples from same sea area were not genetically different. AMOVA results (Table II) revealed that78.65% of variation was attributed among groups, 20.53% was with populations and 0.82% was among

populations with groups. In Φ-statistic, ΦST = 0.80 (P = 0.001), ΦCT = 0.79 (P = 0.030) and ΦSC = 0.04 (P = 0.048) (Table II). The haplotypes network showed a clear division of two groups, consistent with geographical distribution of Group 1 and Group 2 (Fig. 3). The clustering in the UPGMA tree of seven samples also suggested that two main clusters were detected (Fig. 4). The Bayesian tree of haplotypes of COI gene was clustered into two lineages (north and south) and the clustering of haplotypes was conformed to the grouping of samples (Fig. 5).

For Cytb gene, almost FST values were significantly from -0.04 to 0.75 (Fig. 2). AMOVA showed that 60.35% of variation was inter-group and 38.21% of variation was intra-group. Φ-statistic parameters showed that ΦST =0.62 (P = 0.001), ΦCT = 0.60 (P = 0.028) and ΦSC = 0.04 (P = 0.042). The haplotype network was similar to the network of COI gene that north (QD, YC and NT) and south (ZZ, ZH, BH and HK) groups were obviously divided into two parts. Furthermore, UPGMA tree of Cytb gene indicated that seven samples were divided into two subsets (Fig. 4). In the Bayesian phylogenetic trees of haplotypes (Fig. 5), the tree of Cytb gene showed that the haplotypes clustered into two groups. In the Bayesian trees, the clustering of haplotypes was in line with geography distribution.

In the neutrality test, the Tajima’s D and Fu’s Fs values were not significant (Table III). To comparison, we also constructed the Bayesian skyline plot to estimate the population evolution of M. meretrix using the COI and Cytb genes (Fig. 6).

The COI divergence rates was ranged from 0.7%/Ma (Marko, 2002) to 2.4%/Ma (Hellberg and Vacquier, 1999). Because of the lack of the rate data of Cytb gene in mollusk, we used the standard rate from 1.5%/Ma to 2%/Ma (Nei and Koehn, 1983). The Bayesian skyline plots showed that the sample size stayed stable. It was non-significantly showed that seven samples did not have a recent population expansion event.

The genetic parameters and haplotype sizes were in a high level in all seven lacation samples of M. meretrix in the coast of China. This trait was also found in others economical bivalve in China, such as Macridiscus multifarious (Ye et al., 2015), Meretrix lamarckii (Teng et al., 2015), etc. The shared haplotypes of COI and Cytb were widely detected within the groups. It was suggested that samples had a frequently genetic exchange within the group. The high level of genetic diversity may be caused by its short planktonic larvae stage (6 days). As the important part of its life cycle, the planktonic larvae stage was the main way for genetic exchange between different regions. The spawning period was range from May to September (Zhao and Yang, 2014). The ocean currents were actively that promoted the gamete to spread in a no-barriers ocean area. With the influence of the ocean currents, shorten larvae stage also reduced the integration time of different samples to establish a stable location sample. The neutrality test also suggested the samples were in a stable.

The pairwise FST values, haplotypes networks, UPGMA trees of samples and Bayesian trees of haplotypes were provided insights that seven samples of M. meretrix were clustered into two groups. It was obvious that north group was significantly differentiation with south group. Each group had its own private haplotypes and lacked share haplotypes (Fig. 1). Within the groups, it was no significantly divergence. The samples in the same sea area had frequently gene exchange which reduced the gene divergence. Therefore, it was obviously observed that a genetic divergence was happened between the two groups. The clustering results were conformed to the geographic distribution of seven samples.

To explain the results, it was necessary to find out which factors would influence of the genetic of bivalve. In the ocean, stereoscopic spaces lacked effective barriers to block communication between animals. In the case of

bivalve, the gene exchange of different populations only could happen in their planktonic larval stage because the adults were fixed. Therefore, the planktonic period duration and the ocean currents were the important elements that caused the genetic divergence between the populations. The larvae were easily affected by the ocean currents in the long floating. The planktonic period duration was determined the floating distance of larvae and the gene flow.

In this study, in terms of the planktonic period duration, it was shorter than most bivalve (Wang et al., 2008). It was not enough for larvae to across the Yangtze River diluting water and follows the coastal water from the Yellow Sea to the South China Sea (Liu, 2013). Shorten larvae stage did not support a long-distance dispersal, such the report of our previous study of Gomphina aequilatera (Ye et al., 2018). Oppositely, a long larvae stage would help the larvae to cross the long-distance, such the report of Crassostrea sikamea in Hu et al. (2018).

In terms of the ocean currents, the larvae could not follow the ocean currents to distinct places from north to south, because of the Yangtze diluting water (YDC) cut off the currents. The Yangtze diluting water and other currents meet around the estuary which cut down the spread of gametes and reduced the gene flow among groups which eventually resulted in the differentiation between samples and groups (Su and Yuan, 2005). Then, ocean current frontal system also reduced the exchange of larvae to from SBCC to MZCC, and in South China Sea, the ocean current system (GDCC) also formed complex circulation. It meant that the gene flow would not be able to across the coastal water and exchange the gene between north and south. This phenomenon also found in other studies that the genetic variation of bivalve was affected by the ocean currents (Liu et al., 2011; Ni et al., 2012, 2015).

Therefore, the north group could not spread their larvae from north to south. The south group also difficultly got away from the circulation system and exchanges the gametes with the north group. Lacking the gametes exchange between these two group, the genetic variation would be estimated an obviously divergence among groups.

In neutrality test, the Tajima’s D and Fu’s Fs values non-significantly supported the samples of M. meretrix happened a population expansion event. The Bayesian skyline plots showed that the samples of M. meretrix were in stable. It did not appear that any clams were adversely affected by glacial cycles. The analysis showed that the wild sample of M. meretrix was not impacted by the natural factors. The Bayesian skyline plot shown the stable sample size was formed since 0.5 Ka. The stable sample size promoted the gene flow and caused a high genetic diversity.

There were several reasons that threat to the wildlife of M. meretrix. The artificial seeding could not fully satisfy the aquaculture requirement, and the wild seeding was the important supplementary sources. The seeding plants would prefer to utilize the wild seeding for farmers which could reduce the cost of productions, which was possibly that increase the impact on wild populations. The mussel disease was also an assignable factor that decreased the wild population size (Murgarella et al., 2016). The natural habitats of M. meretrix were destroyed by the human activities. As a main bivalve for aquaculture, with these negative factors, it was necessary to investigate the current situation of the natural resources. Some research had done on the M. meretrix in small areas or with little location samples. In China, He et al. (2008) studied five wild and cultured samples in the Yellow Sea by the AFLP markers. The results suggested that these five samples showed significant differentiation. Xue et al. (2006) studied three wild samples with enzymatic in the South China Sea. There were one samples was divergence with others. Less or none researches investigated the genetic relationship of M. meretrix samples between north and south in China. As a maintain aquaculture species, it was necessary to provide insight for using to larval rearing and genetic resource protecting.

Revealing the genetic variation of an aquaculture species, it was necessary to regularly detect the basic genetics parameters of its main wild breeding source. It could provide insight of selecting breeding samples and protect the wild genetic resource. In this study, with distributed Chinese coastline, we investigated that seven samples of M. meretrix were detected an obviously genetic divergence. Based on the geographic distance, seven samples could be easily divided into two groups. Inside the group, the samples had a high-level genetic diversity and lower genetic differentiation. For the fisheries management, it was necessary to establish the genetic protection unit for M. meretrix in the North and South China Sea. Keeping the genetic divergence of the M. meretrix between north and south was a best way to maintain the genetic resource. A species with genetic divergence was useful for the fisheries researchers to solve the disease of hard clam. It was also helpful to select the individuals with good traits in different groups for hybridization.

This work was financially supported by the Natural Science Foundation of Zhejiang Province (Grant No. LQ18D060004), and the 2017 Scientific Research Startup Foundation of Zhejiang Ocean University (Funding No. 12245090318).

The authors declare no conflict of interest.