Population structure of Bathymodiolus manusensis, a deep-sea hydrothermal vent-dependent mussel from Manus Basin, Papua New Guinea

Sample collection and DNA extraction

As part of a larger study looking at multi-species biodiversity and population structure within Manus Basin, Bathymodiolus manusensis were collected from two hydrothermal vent sites (Solwara 8 and South Su; Fig. 1) during the M/V Nor Sky research campaign (June–July 2008; Chief Scientist: S Smith) using an ST200 ROV modified for biological sampling. For COI analyses, B. manusensis were analyzed from two discrete sulfide mounds at each site (Table 1), with 10–43 individuals per mound (Table 2). For microsatellite analyses, up to 142 individuals per locus were analyzed. All sampling was undertaken with the permission of the government of Papua New Guinea and did not involve endangered or protected species.

Mantle tissue was dissected from each individual and preserved in 95% ethanol prior to DNA extraction. Genomic DNA was isolated using a standard Chelex-Proteinase-K extraction (10–30 mg tissue digested with 120 µg Proteinase K (Bioline, Taunton, MA, USA) in 600 µl 10% Chelex-100 resin (Bio-Rad, Hercules, CA, USA)) overnight at 60 °C, heated to 100 °C for 15 min, and centrifuged at 10,000 rpm for 5 min; (Walsh, Metzger & Higuchi, 1991). Extracted DNA was stored at 4 °C until amplification and archived at −20 °C.

COI sequencing and analysis

Bathymodiolus manusensis mitochondrial COI fragments were amplified using the following reaction conditions: 10–100 ng of DNA template was combined with 2 µL 10× PCR buffer (200 mM Tris, pH 8.8; 500 mM KCl; 0.1% Triton X-100; 0.2 mg/ml BSA), 2 mM MgCl₂, 0.2 mM dNTPs, 0.5 µM LCOI1490 and 0.5 µM HCOI2198 primers (Folmer et al., 1994), and 1 unit of Taq polymerase in a 20 µL reaction with the following PCR protocol: initial melting temperature of 94 °C for 240 s; 35 cycles of 94 °C for 15 s, 48 °C for 15 s, 72 °C for 30 s; and a final extension of 72 °C for 300 s. Products were stored at 4 °C until purification.

Fourteen µl of PCR product was incubated with 0.2 µl 10× ExoAP buffer (500 mM Bis-Tris, 10 mM MgCl2, 1 mM ZnSO4), 0.05 µl Antarctic Phosphatase (New England Biolabs, Ipswich, MA, USA), 0.05 µl Exonuclease I (New England Biolabs, Ipswich, MA, USA) at 37 °C for 60 min followed by 85 °C for 15 min to remove unincorporated nucleotides. Sequencing reactions were executed with BigDye Terminator v3 reactions (Applied Biosystems, Foster City, CA, USA). AMPure magnetic beads (Agencourt; Morrisville, NC, USA) were used to remove excess dye, products were analyzed on an ABI 3730xl DNA Analyzer (Applied Biosystems International, Rotkreuz, Switzerland), and chromatograms were edited using CodonCode Aligner (version 3.7.1; CodonCode Corporation, Dedham, MA, USA). Consensus sequences were compared against the NCBI GenBank database to confirm identity when available (Benson, 1997) and sequence alignments were constructed using the MUSCLE alignment algorithm (Edgar, 2004) implemented in CodonCode Aligner. Representative sequences of dominant haplotypes were deposited in GenBank (accession numbers KF498731 –KF498847). Full COI sequences for each individual are provided as FASTA files (Data S1).

Standard summary statistics, including number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), and Fu’s F_S were calculated using DnaSP version 5.10.01 (Librado & Rozas, 2009). To detect potential cryptic species, maximum-parsimony phylograms of aligned mitochondrial sequences were assembled in MEGA version 5 (10,000 replicates; Tamura 3-parameter substitution model determined by Mega 5: Find Best-Fit Substitution Model; Tamura et al., 2011). To visualize potential population structure, statistical-parsimony networks were assembled in TCS version 1.21 (default settings; Clement, Posada & Crandall, 2000). To detect population structure, Arlequin version 3.5.1.2 (Excoffier, Laval & Schneider, 2005) was used to estimate pairwise φ_ST. Sequential Bonferroni was used in all appropriate comparisons to correct for multiple tests (Rice, 1989).

Additional Bathymodiolus manusensis samples were identified in NCBI GenBank (accession numbers KU597590.1 through KU597592.1 from PACMANUS in Manus Basin (Assié et al., 2016); AB101431.1 through AB101434.1 from PACMANUS in Manus Basin (Miyazaki et al., 2004); AB257539.1, AB257541.1, and AB257543.1 from Lau Basin (Miyazaki et al., 2010); and AB255739.1, AB255740.1, AB255741.1, and AB255742.1 from offshore New Zealand (Miyazaki et al., 2010)). A Neighbor-Joining tree (10,000 bootstrap replicates) was assembled from three sequenced individuals from Solwara 8, three individuals from South Su, as well as seven individuals from PACMANUS, three individuals from Lau Basin, and four individuals from offshore New Zealand GenBank sequences using MEGA 7 (Kumar, Stecher & Tamura, 2016).

Microsatellite genotyping and statistical analyses

Eight microsatellite markers (Bm17, Bm22, Bm23, Bm53, Bm63, Bm76, Bm81, Bm83) were amplified from Bathymodiolus manusensis in Manus Basin following methods reported in (Schultz et al., 2010). To test whether these markers provided sufficient power to evaluate the null hypothesis of genetic homogeneity, models of the dataset were implemented in POWSIM (Settings based on observed allele distributions, Supplement 1; Ryman & Plam, 2006). Full microsatellite genotypes for each individual are provided as GENPOP files (Data S2).

Standard summary statistics, including divergence from expected Hardy-Weinberg Equilibrium (HWE) and allelic richness were assessed using GENEPOP (default settings; version 4.0; Rousset, 2008) and Microsatellite Analyzer (version 4.05; Dieringer & Schlötterer, 2003), respectively. Permutation tests were used to determine significant variation in allelic richness (F-stat; default settings; version 2.9.3.2; Goudet, 1995). MicroChecker (version 2.2.3; 1,000 randomizations; Van Oosterhout et al., 2004) was used to detect the potential presence of null alleles, stutter, and large allele dropout. To test for the potential influence of selection, loci were screened using LOSITAN (25,000 simulations; IA and SMM; Antao et al., 2008; Beaumont & Nichols, 1996).

Pairwise genetic differentiation (F_ST) between aggregations, sites, and basins was analyzed using Microsatellite Analyzer. Alpha levels were adjusted via Sequential Bonferroni to correct for multiple tests (Rice, 1989). Structure version 2.3.3 (admixture model, sampling locations as prior distributions; Pritchard, Stephens & Donnelly, 2000) was used to visualize potential population structure. Analyses were conducted with a 1,000,000 step burn-in, 10,000,000 repetitions, and three replicates per level from K = 1–7. Effective population size was estimated based on microsatellite linkage-disequilibrium using LDNe (default parameters; Waples & Do, 2008).

Population structure of Bathymodiolus manusensis in Manus Basin

Bathymodiolus manusensis form a single, coherent population between Solwara 8 and South Su in Manus Basin, Papua New Guinea. No genetic differentiation was detected at any spatial scale using either mitochondrial COI or nuclear microsatellite markers. Despite this apparent lack of population structure, B. manusensis is absent from Solwara 1, a site that occurs between Solwara 8 and South Su and that is within 2.5 km of South Su. Further, Solwara 1 shares many vent-dependent and vent-associated species with Solwara 8 and South Su (Coffey Natural Systems, 2008; Erickson, Macko & Van Dover, 2009; Thaler et al., 2011; Thaler et al., 2014; Plouviez et al., 2013).

The relatively homogeneous distribution of both COI haplotype and microsatellite markers for Bathymodiolus manusensis within Manus Basin is consistent with high gene flow between Solwara 8 and South Su. Similar levels of gene flow were observed in Ifremeria nautilei (Thaler et al., 2011) and Chorocaris sp. 2 (Thaler et al., 2014), although in both cases, the species were also found at Solwara 1. A significant, negative Fu’s F_S is consistent with a recent, rapid expansion, a pattern also observed in other species examined from these sites (Thaler et al., 2011; Plouviez et al., 2013; Thaler et al., 2014). COI haplotype diversity is higher at South Su (0.59 compared to Solwara 8’s 0.45; Table 2), however, there is no consistent pattern of microsatellite richness between the two sites, nor are unique alleles consistently identified at one site over the other. A potential alternate explanation could be that both sites were recently colonized by the same cohort, and though currently isolated, have not had enough time for significant differentiation to accumulate via genetic drift.

In more than ten years of exploration and environmental observations, consisting of at least four research campaigns, neither Bathymodiolus manusensis nor any other mussel in the genus Bathymodiolus has been observed at Solwara 1 (W Saleu, pers. obs., 2014). Visual surveys of the seafloor suggest that adequate substrate (hard basalt surrounding low temperature venting fissures) exists within the Solwara 1 site for B. manusensis to settle (A Thaler, pers. obs., 2008), although the fluid chemistry that might influence mussel recruitment has not been characterized for Manus vents. In a previous study, we identified a similar pattern of presence/absence among populations of Munidopsis lauensis at Solwara 8 and South Su (Thaler et al., 2014). One population of M. lauensis was found at Solwara 8 and South Su, but absent at Solwara 1, while a second population was restricted to samples from Solwara 1 (Thaler et al., 2014). We hypothesized that sweepstakes effects related to the survival and settlement of recruits at vent sites in Manus Basin was responsible for the observed population structure of M. lauensis and that time series sampling would reveal a stochastic, dynamic distribution of these populations throughout the basin (Thaler et al., 2014).

When compared with putative Bathymodiolus manusensis samples from other sites in Manus Basin as well as Lau Basin and offshore New Zealand, there is a similarly high affinity between Solwara 8, South Su, and PACMANUS samples, suggesting a larger, well-mixed Manus population. B. manusensis sequences from outside of Manus Basin were basal to all Manus samples, suggesting that the population within Manus Basin is younger than those from surrounding regions, and that there is a greater retention of propagules within Manus Basin, resulting in less gene flow and increased isolation within Manus Basin. Similar patterns of isolation within Manus Basin were observed in Ifremeria nautilei and Chorocaris sp. 2 (Thaler et al., 2011; Thaler et al., 2014).

That Bathymodiolus manusensis shares the same pattern of presence/absence with one population of Munidopsis lauensis suggests that the apparent exclusion of certain species or populations from Solwara 1 may be the result of a consistent, species- and population-dependent, dispersal barrier, rather than stochastic recruitment events. Other “leaky” dispersal barriers have been observed for hydrothermal vent populations across the equatorial East Pacific Rise (Plouviez et al., 2009; Plouviez et al., 2010; Vrijenhoek, 2010), but those sites were separated by thousands of kilometers, whereas the Manus Basin sites are 2.5–40 km apart. To determine if there is a barrier restricting some, but not all, species (or populations) from recruiting to Solwara 1, we need to sample additional species to identify consistent patterns across multiple taxa and sample the same species at additional time points to establish if observed patterns are temporally stable. The alternative hypothesis that species and populations are adapted to particular environmental conditions that are not always present at a site remains plausible (and not mutually exclusive), especially given well-documented evidence for such circumstances in Alviniconcha species in Lau Basin (Beinart et al., 2012) and the lack of fluid chemistry data from evident and putative mussel habitats at Manus Basin vents.

Implications for management strategies

The limited distribution of Bathymodiolus manusensis and of a Munidopsis lauensis population within Manus Basin underscores the potential complexity of connectivity and habitat availability within vent ecosystems and the value of comprehensive environmental baselines prior to the initiation of an extractive regime (Collins, Kennedy & Van Dover, 2012; Collins et al., 2013; Thaler et al., 2014; Boschen et al., 2016). It is possible that multiple mining events in Manus Basin could affect source–sink dynamics of B. manusensis and other taxa, resulting in regime shifts in vent communities of Manus Basin as has been noted in other marine ecosystems (Scheffer & Carpenter, 2003). A similar phenomenon was observed in Moorea coral reef communities, where persistent disturbance, caused, in this case, by invasive crown-of-thorn starfish resulted in permanent changes in community structure as opportunistic recruits occupied newly exposed ecologic niches (Berumen & Pratchett, 2006).

The potential for regime shifts, where species not present at the disturbance site but occurring at neighboring sites establish a foothold following anthropogenic impacts, creates a challenge for environmental management and mitigation programs. At the very least, there is a need to understand if such a regime shift constitutes a significant adverse impact that should trigger a management response. Managers need to understand the extent of local variation in population and community structure to anticipate cumulative impacts and ecological consequences of regime changes following disturbance.