Comparative population structure of two dominant species, Shinkaia crosnieri (Munidopsidae: Shinkaia) and Bathymodiolus platifrons (Mytilidae: Bathymodiolus), inhabiting both deep‐sea vent and cold seep inferred from mitochondrial multi‐genes

Abstract Deep‐sea hydrothermal vents and cold seeps, limited environments without sunlight, are two types of extreme habitat for marine organisms. The differences between vents and cold seeps may facilitate genetic isolation and produce population heterogeneity. However, information on such chemosynthetic fauna taxa is rare, especially regarding the population diversity of species inhabiting both vents and cold seeps. In this study, three mitochondrial DNA fragments (the cytochrome c oxidase submit I (COI), cytochrome b gene (Cytb), and 16S) were concatenated as a mitochondrial concatenated dataset (MCD) to examine the genetic diversity, population structure, and demographic history of Shinkaia crosnieri and Bathymodiolus platifrons. The genetic diversity differences between vent and seep populations were statistically significant for S. crosnieri but not for B. platifrons. S. crosnieri showed less gene flow and higher levels of genetic differentiation between the vent and seep populations than B. platifrons. In addition, the results suggest that all the B. platifrons populations, but only the S. crosnieri vent populations, passed through a recent expansion or bottleneck. Therefore, different population distribution patterns for the two dominant species were detected; a pattern of population differentiation for S. crosnieri and a homogeneity pattern for B. platifrons. These different population distribution patterns were related to both extrinsic restrictive factors and intrinsic factors. Based on the fact that the two species were collected in almost identical or adjacent sampling sites, we speculated that the primary factors underlying the differences in the population distribution patterns were intrinsic. The historical demographics, dispersal ability, and the tolerance level of environmental heterogeneity are most likely responsible for the different distribution patterns.


Introduction
Deep-sea hydrothermal vents, first discovered on the Galapagos Rift in the eastern Pacific Ocean in 1977 (Corliss et al. 1979), are typically located in mid-ocean ridges and back-arc spreading centers. Vents have extremely high temperatures, high pressure levels, and high levels of toxins (Yang et al. 2013). Cold seeps, situated in subduction zones with ecosystems similar to vents, are areas where chemically modified fluids derived from hydrocarbon reservoirs, methane hydrates, pore waters in sediments, and sites of organic enrichment (such as whale skeletons) are released into the ocean (Van Dover et al. 2002). A few decades ago, sunlight was considered the essential element for life because photosynthesis provides primary production for the biosphere. However, numerous unique "chemosynthetic fauna", which depend primarily on energy supplied by the chemosynthesis of bacterial endosymbionts, are found in the limited environments without sunlight (Tirmizi and Javed 1980;Kenk and Wilson 1985;Hessler and Lonsdale 1991;Kikuchi and Ohta 1995;Watsuji et al. 2010). Knowledge regarding the biological communities around vents and cold seeps led to a profound change in our perception of deep-sea life (Van Dover et al. 2002).
Within deep-sea chemosynthetic ecosystems, particularly hydrothermal vents, studies of population structure and gene flow has recently received much attention for many species, such as giant tubeworms (Hurtado et al. 2004;Coykendall et al. 2011), mussels (Won et al. 2003), swarming shrimp (Teixeira et al. 2012), and provannid gastropods (Kojima et al. 2001;Teixeira et al. 2012). Population diversity studies are also reported for deep-sea chemosynthetic communities at the methane seeps in the Gulf of Mexico and eastern Atlantic (Carney et al. 2006;Olu et al. 2010). From these studies, we are beginning to understand the extent of genetic connectivity among deep-sea populations associated with chemosynthetic ecosystems. However, few studies have been conducted to examine the genetic diversity and population structure for species inhabiting both vents and cold seeps, especially for the dominant species. As we know, relying on chemoautotrophic bacteria for nutritional support, all known species of macroorganisms in hydrothermal vents and cold seeps are highly specialized for a symbiotic lifestyle (Cavanaugh 1983;Fialamedioni and Lepennec 1988;Fisher 1990). The differences between the ecosystems of vents and cold seeps may facilitate the genetic isolation of species in both environments. Therefore, it is particularly important to determine the extent of population differentiation and genetic isolation for these species.
Bathymodiolus (Mytilidae: Bathymodiolus) mussels, which are the dominant species in vents and cold seeps, have been described beyond 18 species. Only three Bathymodiolus species in Japanese waters (B. japonicus, B. platifrons, and B. aduloides) are capable of inhabiting both vents and cold seeps. In particular, the same haplotypes have been found to be shared by B. platifrons from distant localities and highly differentiated environments (Miyazaki et al. 2004). Shinkaia crosnieri (Munidopsidae: Shinkaia) is also considered the dominant member of the fauna inhabiting many deep-sea hydrothermal vents and minority cold seeps, including Hatoma Knoll, the Formosa Ridge off southwest Taiwan, and the Okinawa Trough (Chevaldonne and Olu 1996;Tsuchida and Fujikura 2000;Watabe and Miyake 2000;Martin and Haney 2005;Tsuchida et al. 2007). Previous studies of B. platifrons from two seeps and one vent revealed no genetic differentiation between the mussel populations of either species from cold seeps versus vents and suggested high adaptability to deep-sea chemosynthetic environments (Miyazaki et al. 2004). However, the distribution and population structure of S. crosnieri have been reported only at hydrothermal vents in the southern Okinawa Trough (Tsuchida et al. 2003). The two dominant species (B. platifrons and S. crosnieri) are ideal models as they exhibit various population distances and inhabit different habitat types.
In this study, three mitochondrial DNA fragments (the cytochrome c oxidase submit I (COI), cytochrome b gene (Cytb), and 16S) were employed as a mitochondrial concatenated dataset (MCD) to evaluate the genetic diversity, genetic structure, and demographic history of B. platifrons and S. crosnieri. The current study also sought to explore the adaptability of the two dominant species to different deep-sea chemosynthetic environments.

Materials and Methods
The experiments were performed in accordance with the recommendations of the Ethics Committee of the Institute of Hydrobiology, Chinese Academy of Sciences. The policies were enacted according to the Chinese Association for Laboratory Animal Sciences and the Institutional Animal Care and Use Committee (IACUC) protocols.

Sampling, identification, and DNA extraction
A total of 64 specimens of S. crosnieri from three vents and one cold seep habitat, and 52 specimens of B. platifrons from two vents and one cold seep habitat were collected by the manned submersible "Jiaolong" during June 2013 and the ROV "Faxian" during April 2014. All the sample information and photographs used in this study are provided in Table 1 and the sampling sites are mapped in Fig. 1. All the specimens were frozen and preserved at À80°C or in 100% ethanol and deposited in the Institute of Oceanology, Chinese Academy of Science.
The species-level morphological identification was performed based on the original morphological descriptions, locality data, and additional information of S. crosnieri (Baba and Williams 1998;Chan et al. 2000) and of B. platifrons (Hashimoto and Okutani 1994;Barry et al. 2002).
The total genomic DNA of S. crosnieri was extracted from a small piece of abdominal muscle tissue using the E.Z.N.A. â Tissue DNA Kit (OMEGA, Wuhan, China) following the manufacturer's instructions. Several bivalve species have an unusual mode of mtDNA inheritance known as doubly uniparental inheritance (DUI) (Zouros et al. 1994a;Zouros et al. 1994b;Garrido-Ramos et al. 1998;Plazzi and Passamonti 2010). DUI species have two highly differentiated mitochondrial genomes, one of which follows the standard mode of maternal inheritance (known as type F) and the other is transmitted through sperm and found almost entirely in the male gonad (known as type M) (Mizi et al. 2005). Phylogenetic and population analyses require comparisons between orthologous sequences and M-or F-type genes under DUI that are not orthologous sequences (Plazzi and Passamonti 2010). Therefore, the total genomic DNA of B. platifrons was obtained from a small section of adductor or foot muscle, which carries very little M-type mtDNA in DUI species (Grant and Bowen 1998). Then, the extracted DNA was suspended in distilled water and stored at À20°C.

PCR amplification and sequencing
Three mitochondrial DNA fragments, COI, Cytb, and 16S were amplified using the primers listed in Table S1 (Palumbi et al. 1991;Folmer et al. 1994;Merritt et al. 1998;Baco-Taylor 2002;Passamonti 2007;Radulovici et al. 2009). The final reaction volume of 30 lL contained 21.125 lL of sterilized ultrapure water, 3.0 lL of 103 PCR buffer (including MgCl 2 ), 1.5 lL of each primer (10 mmol/L), 1.5 lL of dNTPs (2.5 mmol/L each), 0.375 lL of Taq DNA polymerase (2.5 U/lL, TaKaRa Bio, Shanghai, China), and 1.0 lL of DNA template (50-100 ng/lL). The cycling parameters were 94°C for 5 min; 32 cycles of 94°C for 40 sec, 43-48°C for 30 sec and 72°C for 1 min, and a final elongation step at 72°C for 10 min. The PCR products were visualized on 1.2% lowmelting agarose gels stained with ethidium bromide. Then, the products were purified and sequenced on an ABI3730 XL sequencing system using the primers described in Table S1.

Molecular data analysis
The sequences obtained in both directions were checked by the sequence peak figure and then assembled based on the contigs using the DNASTAR Lasergene package (DNASTAR, Inc., Madison, WI, USA). The sequences were aligned and trimmed to the same length using the software package MEGA 5.0 (Tamura et al. 2011). Then, the MCD was used directly for the subsequent population analyses. The average genetic distances within and between populations were estimated according to the Kimura 2-parameter (Kimura 1980) model in MEGA 5.0 (Tamura et al. 2011). Genetic diversity is reflected in the measures of nucleotide diversity (p) and haplotype diversity (h) (Nei 1987), and the values for each population were calculated using DnaSP 5.10 (Librado and Rozas 2009). The maximum likelihoods (ML) for phylogenetic analyses were assembled in PhyML 3.0 (Guindon and Gascuel 2003) with 1000 replicates, and the most appropriate model of DNA substitution, which is TIM2+I+G for S. crosnieri and HKY+I for B. platifrons, was identified by ModelTest 3.7 (Posada and Crandall 1998). Alvinocaris longirostris (GenBank: AB821296) and Geothelphusa dehaani (GenBank: AB187570) were set as outgroups for S. crosnieri. Mytilus californianus (GenBank: NC_015993) and Mytilus edulis (GenBank: NC_006161) were set as outgroups for B. platifrons.
Median-joining network (MJN) analysis was performed with Network 4.6 (Bandelt et al. 1999) to depict the relationships among all the haplotypes. The pairwise genetic divergences between the populations were estimated using F-statistics (F ST ) with 10,000 permutations, based on the distance method in Arlequin 3.5 (Excoffier and Lischer 2010). Analysis of the molecular variance analysis (AMOVA) was used to analyze the hierarchal population structure in Arlequin.
The population demography (e.g., bottlenecks or expansions) for the vent and seep populations were examined using two different approaches. First, the demographic history was investigated by comparing the mismatch distributions in each habitat-type sample with those expected in stationary and expanding populations using DnaSP 5.10 (Librado and Rozas 2009). In addition, we tested the goodness-of-fit of the actual distributions with the expected distributions using a model of population expansion by calculating the Harpending's raggedness index (r) (Harpending 1994) and by assessing significance with 1000 permutations. Second, Tajima's D (Tajima 1989) was also applied to seek evidence of demographical expansions in Arlequin 3.5 (Excoffier and Lischer 2010). In addition, gene flow (N m ) was evaluated from GammaSt (Nei 1982), and F ST (Slatkin 1985;Hudson et al. 1992) was evaluated using DnaSP 5.10 (Librado and Rozas 2009). The geographic distances between populations were estimated using Google Earth 4.3 based on the longitude and latitude.

Sequence information
All three mitochondrial gene fragments were successfully amplified for 64 S. crosnieri specimens and 52 B. platifrons specimens. Meanwhile, COI (GenBank accession nos. KR003111-KR003157), Cytb (GenBank accession nos. S. crosnieri: genetic diversity and population structure The number of haplotypes (H), the haplotype diversity (h), and nucleotide diversity (p) for each population of S. crosnieri are presented in Table 2. All populations showed high haplotype diversity and there was no shared haplotypes except one each for the vent and cold seep populations. All populations also showed high haplotype diversity, and several shared haplotypes were detected among the vent and cold seep populations. A comparison of the genetic diversity showed that the haplotype and nucleotide diversities were lowest in the cold seep population (h = 0.972, p = 0.0062) and highest in the vent WTS1 population (p = 0.0062).
Based on the Kimura 2-parameter model, the overall average genetic distance among individuals was 0.0112 for S. crosnieri. The mean genetic distances within and between each of the populations are shown in Table 3. The overall mean intrapopulation genetic distance of the WTS1 population (0.0062) was the largest among the four populations and was smallest for the WKC population (0.0032). The genetic distances among WTS1, WTS2, and NIS were very small, but the genetic distances between WKC and all the vent populations were higher.
The median-joining network analysis had a weblike topology, with many singletons connected through multiple nodes, indicating high genetic variability for S. crosnieri ( Fig. 2A). It was found that haplotypes from the vent and seep populations could be separated and did not share common haplotypes. Moreover, two obvious clades (Vent and Cold Seep) were also identified from the ML tree (Fig. 4A).
For population structure analysis, all samples were divided into vent and seep groups. An AMOVA analysis revealed that most of the molecular variance of S. crosnieri occurred among groups (79.70%), whereas the variance among the intragroup populations (0.35%) and within populations (19.95%) was relatively small (F ST = 0.801, P < 0.001), indicating a high level of geographical population structure (Table 4). In addition, all pairwise F ST comparisons between the vent and seep populations were statistically significant (P < 0.001). The values of Nm based GammaSt (from 1.08 to 29.93) and F ST (from 0.27 to 43.98) indicated the fluctuation of gene flow levels among the four S. crosnieri populations (Table S2). The level of gene flow between the vent and seep populations was very limited, which was consistent with the statistically significant F ST values among them.
The shape of the mismatch distribution of S. crosnieri was approximately unimodal for the whole vent dataset (Fig. 3A). Nevertheless, the distribution was obviously Harpending's raggedness index; Significant: * = P < 0.05; ** = P < 0.01; *** = P < 0.001. ragged and multimodal for the seep WKC population (Fig. 3B) and the whole dataset (Fig. 3C). In addition, none of the values for Harpending's raggedness index were significant, and the Tajima's D values were statistically significant for only the whole vents dataset (P < 0.01) ( Table 2).

B. platifrons: genetic diversity and population structure
The values of H, h, and p for each population of B. platifrons are presented in Table 2. All populations of B. platifrons showed high haplotype diversity and several shared haplotypes were detected among the vent and cold seep populations. A comparison of the genetic diversity of B. platifrons detected a similar level of haplotype and nucleotide diversities among all populations. The overall average genetic distance among individuals was 0.0021 for B. platifrons. The mean genetic distances within and between each of the populations are shown in Table 3. The overall mean intrapopulation genetic distance of the WKC population (0.0021) was the largest among the three populations, and that of NIS1 (0.0018) was the smallest. Therefore, the genetic distances among all populations were very small.
The median-joining network analysis had only several singletons connected through multiple nodes for B. platifrons (Fig. 2B). We found that individuals of B. platifrons from the vent and seep populations shared haplotypes. In particular, three common haplotypes were shared by the seep WKC and vent NIS1 populations. One common haplotype was shared by the vent NIS1 and NIS2 populations. In addition, the ML tree also showed that some individuals from the vent and seep populations occupied the same clades (Fig. 4B).
For population-structure analysis, all samples were divided into vent and seep groups. An AMOVA analysis revealed that most of the molecular variance of B. platifrons occurred within populations (89.65%), whereas  All samples were divided into vent and seep groups. For S. crosnieri, the vent group included the WTS1, WTS2, NIS populations, and the seep group included the WKC population. For B. platifrons, the vent group included NIS1, NIS2 populations, and the seep group included the WKC population.
variance among populations within groups (7.35%) and among groups (3.00%) was relatively small (3.00%) (F ST = 0.103, P < 0.001), indicating a low level of geographical structure (Table 4). In addition, significant F ST was found only between WKC and the vent NIS2 population (P < 0.001) ( Table 3). The values of Nm-based GammaSt (from 4.78 to 11.17) and F ST (from 2.37 to 9.96) indicated a slight stabilization of the gene flow levels among the three B. platifrons populations (Table S2). The gene flow levels among the vent populations were obviously large. Especially for vent NIS2 and the seep WKC populations of B. platifrons, the gene flow was rather weak, which was also consistent with the statistically significant F ST values between them. All the mismatch distribution shape of B. platifrons was obviously unimodal (Fig. 3D, E, F). In addition, none of the values for Harpending's raggedness index was significant, and all the Tajima's D values were statistically significant (P < 0.01) ( Table 2).

Discussion
Comparative genetic diversity and population structure Based on the haplotype and nucleotide diversities, all populations of S. crosnieri and B. platifrons showed high haplotype diversities but low nucleotide diversities, which are normally associated with an expansion model (Grant and Bowen 1998). All the vent populations showed higher genetic variability than the seep population for S. crosnieri, whereas a similar level of the genetic variability was detected among B. platifrons populations.
The AMOVA analysis and pairwise F ST values indicated that high levels of population structure and genetic differentiation occurred between the vent and seep populations for S. crosnieri, but low levels of population structure and genetic differentiation occurred among the B. platifrons populations (although the distance between vent and seep are above 983.8 km), except NIS2-WKC, which was probably related to small sample size. For S. crosnieri, the results were also supported by a weblike statistical haplotype network and ML tree analyses. The lack of interpopulation contact might be related to its ecological habits (vent and seep) and the long distance of the distribution areas (the distance between the vent and seep are more than 983.8 km). This phenomenon could be proven by the lack of gene flow between vent and seep populations but higher gene flow levels among vent populations, especially for WTS1-WTS2. However, for B. platifrons, the statistical haplotypes network and ML tree analysis showed that the haplotype distribution did not mirror the geographical origin of the populations because several haplotypes were shared among individuals from different habitat populations. In addition, gene flow was also detected among the populations. The results suggested that extrinsic restrictive factors for the distributions were not related to environmental types (vents vs. seeps) Fujikura et al. 2000).
The shape of the mismatch distribution of pairwise differences was ragged or multimodal for the populations at stationary demographic equilibrium, but it is typically smooth or unimodal for populations that have passed through a recent expansion or bottleneck (Rogers and Harpending 1992;Schneider and Excoffier 1999). The smooth and approximately unimodal shape of the mismatch distribution for the whole vent dataset of S. crosnieri populations showed that vent populations are likely to have passed through a recent expansion or bottleneck, but cold seep WKC S. crosnieri populations did not because of the multimodal shape of the mismatch distribution. The results demonstrated that S. crosnieri populations had different demographic histories in vent and cold seep environments. Significant D values may be due to factors such as population expansion and bottleneck (ArisBrosou and Excoffier 1996). The significant negative D values for the entire vents dataset also could support the results. Nevertheless, for B. platifrons, both the obviously unimodal shape and the statistically significant negative D values showed that all populations probably passed through a recent expansion or bottleneck. Generally, species living in patchy, fragmented, ephemeral habitats are expected to harbor less genetic diversity (Tunnicliffe and Juniper 1990) and are likely to have encountered an expansion or bottleneck. Moreover, cold seeps are relatively stable, whereas some vents persist for only a few decades (Tunnicliffe and Juniper 1990). In this study, the S. crosnieri vent populations have passed through an expansion or bottleneck but the cold seep populations have not, indicating that the vent habitats were more fluctuant than the cold seep habitats.
A previous study for Bathymodiolus mussels has also demonstrated extensive intraspecific genetic exchanges between vent and seep sites across thousands of kilometers (Van Dover et al. 2002;Miyazaki et al. 2004;Iwasaki et al. 2006). Our analysis also indicated that genetic exchanges occurred among three populations and was largest between the seep WKC and vent NIS1 populations. Based on these results, a homogeneity pattern was fit for the B. platifrons populations. The absence of differentiation between the seep and vent mussels demonstrated that the primary intrinsic factor can most likely be ascribed to the different physiological tolerance of the mussels to pressure Fujikura et al. 2000).

Comparison of different population distribution patterns
In our analyses, different population distribution patterns for the two dominant species were detected. However, the different population patterns were related to extrinsic restrictive factors, such as environmental heterogeneity, geological histories, and oceanic currents, as well as intrinsic factors, such as historical demographics, dispersal ability, and physiology (Iwasaki et al. 2006). The extrinsic restrictive factors for the two species were almost identical because the sampling sites were identical or adjacent. Thus, intrinsic factors are the primary factors underlying the different population distribution patterns.
First, the two species belong to different kingdoms and present different historical demographics. Differences in historical demography always have a great contribution to the different population distribution patterns. Therefore, it probably was the primary intrinsic factor. Second, larvae are the primary dispersal vector, especially during the planktotrophic developmental stage (Lutz et al. 1986;Iwasaki et al. 2006;Thaler et al. 2014). Munidopsids appear to produce very large eggs and brooded larvae with yolk sacs, consistent with lecithotrophy, and possibly a limited dispersal capability compared to the planktotrophic larvae with high dispersal capability for the bathymodioline mussels (Won et al. 2003;Miyazaki et al. 2004). At the scale of 1000 km, one could therefore expect to see differentiation between the seep and vent sites for S. crosnieri, as has been shown at a similar scale for the tubeworm Riftia pachyptila in the East Pacific (Coykendall et al. 2011).
Finally, the levels of genetic differentiation were also likely to be related to physiology. Miyazaki et al. suggested high adaptability of these Bathymodiolus species to deepsea chemosynthetic environments (Miyazaki et al. 2004). Our analysis also showed that B. platifrons had high adaptability to different habitats. Therefore, a possible explanation was that the higher degree of tolerance to environmental heterogeneity of B. platifrons, in comparison to S. crosnieri, should be responsible for their different population distribution patterns. However, this proposed explanation requires further proof of physiology tests.
This study provides the first documentation of the different population genetic structure and distribution patterns for S. crosnieri and B. platifrons inhabiting both deep-sea vents and cold seep. Nevertheless, caveats are necessary with respect to the limited sample locations and small sample sizes. More sample sites and molecular markers, especially nuclear loci, should be employed in further research.

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Table S1. PCR primers used in the present study. Table S2. Gene flow (N m : below the diagonal) was evaluated from the values of GammaSt and F ST for the Shinkaia crosnieri and Bathymodiolus platifrons populations.