Highly divergent mussel lineages in isolated Indonesian marine lakes

Sample collection

Collection was conducted under research permits #3246/FRP/SM/VII/2012, and #0094/FRP/SM/V/2009 provided by the Indonesian Institute of Sciences (LIPI) and the Indonesian State Ministry of Research and Technology (RISTEK) to LE Becking. The mussels were collected from three marine lakes in Berau, East Kalimantan, Indonesia (Fig. 1 and Table 1): Kakaban lake (Lake 1), Lake Haji Buang (Lake 2), Tanah Bamban lake (Lake 3). Mussels were abundant in the lakes displaying dense mussel-beds along the shorelines, covering rock or mangrove roots. Yet, these species of Brachidontes have previously not been found outside of lakes in the Indo-Pacific (Colin, 2009; Goto, Tamate & Hanzawa, 2011). After extensive searching along the coasts of the islands in Berau, we found small populations of mussels in two coastal mangrove systems with a mixture of mangroves species Rhizophora mucronata and Avecinnia alba, here referred to as Mangrove 1 (Maratua island) and Mangrove 2 (Samama island) (Fig. 1 and Table 1). The mussels from these two areas were burrowed in the sand. All specimens were conserved in 96% laboratory grade EtOH and stored at 4 °C. For a full description of the Berau marine lakes, see Tomascik & Mah (1994) and Becking et al. (2011). The Berau marine lakes have a tidal regime, which is typically delayed (ranging from 20 min to 4 h) and dampened (ranging from 10–80%) compared to the adjacent sea (Becking et al., 2011). Based on the level of tidal dampening and delay, marine lakes can be ranked by their degree of connectivity to the surrounding sea (Hamner & Hamner, 1998), i.e., Lake 1 was the most isolated lake, followed by Lake 2 and Lake 3. Additional COI sequences of Brachidontes spp. from Palau were obtained from GenBank (accession numbers: AB465561, AB465566, AB465569). Palau is located at approximately 1,900 km north-east from Berau, East Kalimantan. For a full description of the Palauan marine lakes, see Hamner & Hamner (1998) and Colin (2009).

DNA extraction, gene amplification and sequencing

Posterior adductor muscle tissue was sampled from Brachidontes spp. and incubated at room temperature for 24 h in a lysis buffer consisting of 250 mM EDTA, 5% SDS, and 50 mM Tris (pH = 8) (Holland, 1993). Total DNA was purified using DNeasy Blood and Tissue Kit (Qiagen) following the manufacturer’s protocol. COI is an informative marker for high resolution phylogeographic analysis (Lee & Foighil, 2004; Aguirre, Perez & Sirch, 2006; Smietanka et al., 2014). Partial COI was amplified using the CB1F/CO1R primers of Goto, Tamate & Hanzawa (2011) designed for the female type and the jgLCO/jgHCO primers from Geller et al. (2013) for marine invertebrates. Two individuals per lineage were amplified 28S and 18S, using the primers D23F/D6R (Park & Ò Foighil, 2000) and 22F/1789R (Medlin et al., 1988), respectively. These markers were used only on this subset of individuals to assess the phylogenetic context of the lineages within the Mytillidae and we did not expect to find intra-specific or intra-lineage variation with 28S and 18S (Medlin et al., 1988; Park & Ò Foighil, 2000). PCR reactions for all primer pairs were performed in 25 µL volumes containing 0.8 µL 25 mM MgCl₂, 3 µL dNTP’s (1 mM each), 2.5 µL 10× PCR Buffer (SpheroQ), 0.8 µL BSA (10 mg/mL), 1 µL of both primers (10 µM), 0.2 µL Taq polymerase (5 units/µL) and 1 µL template DNA. The PCR program of the COI primers consisted of an initial denaturation step of 94 °C for 3 m followed by 35 PCR cycles of 94 °C for 45s, 49 °C for 30s and 72 °C for 120s, with a final extension step of 72 °C for 5 m. The PCR program of the 18S and 28S primers was identical to the previously described except for the annealing step at 47 °C for 35s. Quality of PCR products was assessed using gel electrophoresis on 1% agarose gels. If PCR product quality was insufficient, a new PCR reaction was performed in a 25 µL volume containing PCR Beads (Illustra, GE Healthcare) using 21.6 µL ddH₂O, 0.2 µL of both primers (10 µM) and 3 µL template DNA, and an identical PCR cycling program. PCR products were purified and sequenced by Macrogen Inc. (Korea and The Netherlands).

Doubly uniparental inheritance (DUI)

Some species of the family Mytilidae display a type of mitochondrial inheritance called ‘doubly uniparental inheritance’ (DUI) (Zouros et al., 1994) in which two different mitochondrial lineages, a male (M-type) and female (F-type) lineage, are present in male individuals (Terranova et al., 2007; Goto, Tamate & Hanzawa, 2011). Male and female mitochondrial lineages can be highly diverged (Saavedra, Reyero & Zouros, 1997; Kenchington et al., 2002). We avoided obtaining M-type sequences by sampling muscle tissue, avoiding gonads, and using specific primers for F-type sequences. With an additional check, the distinction between F- and M-type sequences in our data was made by constructing a phylogeny. We obtained two M-type sequences (Genbank accession numbers KX688100 and KX688101). Subsequently, the analyses in the present paper were done only on the female lineages, which is standard in phylogeographic studies of mussels (e.g., Goto, Tamate & Hanzawa, 2011; Smietanka et al., 2014).

Genetic data analysis

The bivalve origin of the obtained sequences was verified through BLAST searches (http://blast.ncbi.nlm.nih.gov/Blast.cgi). From the blast results for the 28S and 18S sequences, the most similar sequences were selected and included in the subsequent phylogenetic analysis. Sequences were aligned and handled in CodonCode aligner v4.1.1 (CodonCode Corporation, USA) and DAMBE 5.2.15 (Xia & Xie, 2001; Xia, 2013). The best-fit DNA substitution model was selected by the Akaike Information Criterion deployed in jMODELTEST v. 2.1.6 (Darriba et al., 2012). These models (COI: Hasegawa-Kishino-Yano (Hasegawa, Kishino & Yano, 1985); 18S: Kimura 2-parameter (Kimura, 1980); 28S: Tamura-Nei (Tamura & Nei, 1993); all with a discrete gamma distribution) were used for subsequent maximum likelihood phylogeny inferences implemented in MEGA v6.06 (Tamura et al., 2011) using a heuristic search with 1,000 bootstrap replicates. The resulting phylogenetic trees were visualized in FigTree v1.4.0 (Rambaut, 2009). Separate maximum likelihood analyses were carried out for each marker. Net nucleotide distances between lineages were calculated in MEGA v6.06 using the models obtained with jMODELTEST. Net nucleotide divergence corrects for discrepancies between gene divergence and population divergence due to ancestral polymorphism in populations (Edwards & Beerli, 2000), since it subtracts the average within-group divergence from the observed between-group estimate. Estimates of genetic variation in samples pooled per location were obtained as haplotype diversities h (Nei, 1987), nucleotide diversities π (Tajima, 1983; Nei, 1987), and mean pairwise differences using Arlequin version v.3.5 (Excoffier & Lischer, 2010). In order to create a haplotype network, separate maximum likelihood tree analyses were carried out for each lineage (following the description above), and the resulting phylogenetic trees were used as input for Haploviewer (Salzburger, Ewing & Von Haeseler, 2011).

Two methods were used to test for signatures of recent population expansion. First, Tajima’s D tests of selective neutrality (Tajima, 1983) were carried out in Arlequin to compare the observed numbers of pairwise nucleotide differences between haplotypes in a sample with expectations under an infinite-sites model of sequence evolution, and under assumptions of selective neutrality and stable population size. Significance was tested by generating 10,000 random permutations. Second, mismatch distributions were calculated in Arlequin (Excoffier & Lischer, 2010) and DnaSP (Rozas et al., 2003) to test for signatures of demographic expansion and to test the null hypothesis of population growth. The observed distribution of pairwise differences between sequences was compared with a theoretical distribution, as expected under a sudden expansion model (Rogers & Harpending, 1992) computed in DnaSP (Rozas et al., 2003). To test the fit of the sudden-expansion-model, the sum of squared deviations (SSD) between the observed data and theoretical model was calculated in Arlequin. Harpending’s raggedness index (rg) was used to determine the smoothness of the observed mismatch distribution, which can be used to distinguish between expanded and stationary populations (Harpending et al., 1993). Expanding populations generate smooth and unimodal distributions, while more stationary populations the mismatch distribution becomes more ragged and erratic. The value of the raggedness index will be low and non-significant in expanding populations, while it is usually high and significant in stationary populations (Harpending, 1994; Harpending et al., 1993).

Morphometric analysis

Brachidontes spp. mussels were photographed in a standardized orientation for geometric morphometric analyses. In total 172 digital images were stored as Nikon RAW format (.nef) and converted to 3,008 × 2,000 pixel JPEG images using Photoshop 5.0 (Adobe). JPEG images were sampled into TPS files using the program tpsUTIL (Rohlf, 2010a). Shell outlines were used to capture variation in shell shape of Brachidontes spp. We used a sliding semi-landmark analysis, in which semi-landmarks are allowed to slide along the outline of a shell in order to find the position that optimally matches the positions of corresponding semi-landmarks in a consensus specimen (Bookstein, 1997; Adams, Rohlf & Slice, 2004). Shell outlines were drawn as curves and digitized as 68 semi-landmarks at equal distance using tpsDig (Rohlf, 2010b), using the beak of the mussel (umbo, see Fig. S1) as a standardized starting point for drawing an outline. A “sliders file” indicating sliding semi-landmarks was made using tpsUtil (Rohlf, 2010a). To standardize for size and orientation we used tpsRelw (Rohlf, 2010c) with Generalized Procrustes Superimposition (Rohlf, 1999). Residuals from the superimposition were analysed with the thin-plate spline interpolating function, producing principal warps, followed by relative warp (RW) analysis. The program TpsRelw was used to obtain centroid size (Bookstein, 1991) and RW scores for each individual. RW axes are analogous to the eigenvectors of principal component analysis, which combine the major patterns of shell shape variation in the data. Repeatability of RW axes was tested using regression analysis and a non-parametric analysis of similarity in PAST 2.11 (Hammer, Harper & Ryan, 2001) of RW scores extracted from 17 specimens of Brachidontes spp., which were independently photographed. RW axes were considered repeatable when they showed a non-significant and close to zero R-value in the analysis of similarity and a strong (r > 0.7) and significant correlation. Only repeatable relative warp axes were included in further analyses of shell shape variation. Correlations of RW scores with centroid size were tested, and if significantly correlated, we used residuals of the regression in analyses of shell shape. Significant differentiation between populations was tested using a non-parametric analysis of similarity (ANOSIM, 10,000 randomizations) (Clarke, 1993) based on Euclidian distance as implemented in PAST v2.11 (Hammer, Harper & Ryan, 2001).

Genetic analysis

A total of 103 COI sequences of 516 bp (aligned length) were collected, resulting in 36 F-type haplotypes (Genbank accession numbers KX346179 –KX346214). No stop codons, indels or double peaks were observed in these sequences. The COI sequences represented four major lineages that were strongly supported by maximum likelihood analyses and these lineages corresponded with locality (Fig. 2 and Table 2). Two of these four divergent lineages correspond to ‘lineage A’ and ‘lineage B’ from the marine lakes in Palau (Goto, Tamate & Hanzawa, 2011) which naming we adopt (Fig. 2). In addition, two new lineages were detected, which we have named C & D. In the present study lineage A was present in Lake 3 (Tanah Bamban Lake), while lineage B was present in lake 1 (Kakaban Lake), lineage C was present in both mangrove localities (Samama and Maratua), and lineage D was present in lake 2 (Haji Buang Lake). A subset of individuals of Lineage A–D were sequenced for a 976 bp long fragment of 28S (aligned length, Genbank accession number KX346215 –KX346217) and 18S (1,323 bp, KX346218 and KX346219) (Figs. S2 and S3) this resulted in identical sequences for Lineage A–C, which fall in the clade Brachidontes s.s. established by Trovant et al. (2015) and are most closely related to Brachidontes spp. from marine lakes in Palau. Though the lineages A & B occur in both Indonesia and Palau, not one haplotype was shared between the different marine lakes populations of Brachidontes spp. from Indonesia and Palau. The mussels from Lake 3 are 1.8% diverged from Lineage A from Palau, and the mussels from Lake 1 are 0.9–2.8% diverged from Lineage B from Palau (Table 2). Lineage D was the most divergent from the other three lineages: 67–74% corrected divergence in COI (Table 2), 15% in 28S, and 5% with 18S. Blast results of COI, 28S, and 18S of Lineage D did not show high similarity (all less than 90%) with any species in the database.

Each lake harbored only one lineage (Fig. 2), while no lineages were shared among lakes in East Kalimantan. In contrast, the two mangrove locations (Mangrove 1 & 2) shared the same lineage (C) and shared haplotypes. The nucleotide diversity, haplotype diversity and mean pairwise difference (Table 1) were lowest in the lakes (π: .1–0.2%; h: 0.3–0.7%; pairwise difference: 0.4–1.7) and highest in the open populations of the mangroves (π: 0.5–75%; h: 0.9; pairwise difference: 2.5–3.7).

Demographic analysis

In each lake the haplotype network was starlike shape with a common central haplotype different from most of the other haplotypes by one or a few nucleotide substitutions (Fig. 2). Mismatch distributions were clearly unimodal for the three lake populations (Fig. 3), typical of populations that have undergone a recent bottleneck in population size (e.g., during colonization events) followed by rapid population expansion (Slatkin & Hudson, 1991; Rogers & Harpending, 1992).

This pattern was congruent with the raggedness index, which was low and non-significant, and with the Tajima’s D statistic, which for lakes 1, 2, and 3 was negative and for lakes 1 and 2 differed significantly from expectations under a neutral model of evolution assuming constant population size (Table 1). By contrast, the mangroves did not have a single dominant haplotype and displayed a multimodal mismatch distribution (Fig. 3), which significantly deviated from an expansion model (Table 1). The mean number of pairwise differences among the mangroves is higher than that in the lakes and assuming a constant mutation rate, this suggest that if a population bottleneck has occurred in the mangroves, this would have happened longer ago in the mangroves than in the lakes (Fig. 3 and Table 1). A significant Tajima’s D could indicate a recent population expansion, but could also be a result of a selective sweep (Tajima, 1989a; Tajima, 1989b; Aris-Brosou & Excoffier, 1996; Peck & Congdon, 2004). Likewise, a unimodal mismatch distribution can also be the result of a selective sweep. Here we take a comparative approach utilizing the same markers in species with a very similar biology. Hence, we assume that differences in genetic signatures between samples are most likely due to differences in demographic histories, and not selection on the markers.

Morphometric analysis

The relative warp (RW) axes 1–2 were significantly repeatable, explaining 90.27% of the total observed variation in shell outline shape (Fig. 4). The residuals of the RW axes 1 and 2 were used to correct for the correlation with centroid size (RW1: R² = 0.4095, p = 0.0001, RW2: R² = 0.0285, p = 0.0296). Figure 4 shows a scatter plot of the individuals along the two first axes of variation in shell outline shape, indicating how 90.27% of total observed variation in shell outline is distributed along these axes. Variation in shell outline was found in the position of the umbo relative to the longitudinal axis of the shell (axis 1) and the shell length to width ratio (axis 2). Overall shell shape variation among populations was highly significant (One-Way Permanova, global F =21.68, p = 0.0001) and all pairwise comparisons between locations showed significant differences (Table 3). A brief morphological description of each of the lineages is provided in S5.

Divergent Brachidontes lineages in Indonesia and Palau

We detected four mitochondrial female-type lineages based on COI sequences in Brachidontes spp. mussels sampled in marine lakes and adjacent mangrove habitats in Indonesia, viz. lineage A, B, C and D (Fig. 2). Two major lineages (A and B) are closely related to the lineages A and B, respectively, reported by Goto, Tamate & Hanzawa (2011) for Brachidontes spp. from marine lakes in Palau. Considering the degree of genetic divergence between the four lineages (14–74%), and corresponding morphological differentiation, the lineages probably represent four distinct species. Furthermore, the genetic distances between lineages A–D are comparable to those reported between (cryptic) species of other mussels (Lee & Foighil, 2004; Lee & Foighil, 2005; Terranova et al., 2007; Cunha et al., 2011) Lineage A–C fall within the clade Brachidontes s.s. established by Trovant et al. (2015) and are closely related to Brachidontes spp. from marine lakes in Palau, as well as Brachidontes variabilis, and Brachidontes mutabilis (Fig. 2, Fig. S3). In contrast, Lineage D does not fit within the three proposed clades of Brachidontinae and probably represents a distinct genus.

The patterns of COI genetic differentiation of the Brachidontes spp. lineages were largely congruent with shell outline shape differentiation. There was, however, one exception: the two mangroves are genetically part of the same lineage and share haplotypes, yet morphometrically they differ significantly. Morphology can provide a proxy for underlying genetic variation, adaptive evolution, as well as phenotypic plasticity related to environmental regimes (e.g., Luttikhuizen, Drent & Baker, 2003; Dawson & Hamner, 2005; Mariani, Peijnenburg & Weetman, 2012; Burridge et al., 2015). It is likely that part of the morphometric variation in shell outline shapes of Brachidontes spp. mussels are the result of highly distinct environmental conditions at each locality (Table 1). Worldwide 31 species of Brachidontes have been described, but the phylogenetic positions of the different species within this genus remains unclear with several reports of the occurrence of cryptic species (Lee & Foighil, 2004; Lee & Foighil, 2005; Aguirre, Perez & Sirch, 2006). Further taxonomic review of these Indo-Pacific lineages is clearly needed.

Isolated populations in Indonesian marine lakes

In our study we found that each marine lake harbors a population of mytilid bivalves, but each lake is dominated by a different species that is morphologically and genetically distinct. It is remarkable that the lakes are at such short distance from each other (2–6 km), yet do not share lineages, while the two mangrove locations 20 km apart do. The overall high level of variation is most pronounced in the mangroves, which is interpreted as evidence of a larger population, more stable demographic history and/or more connected systems than the lake systems. Dawson & Hamner (2005) found that the genetic distance between marine lake and lagoon populations of the jellyfish Mastigias papua in Palau was highly correlated with the degree of physical isolation of the lake to the adjacent sea, and not the actual geographic proximity of populations to each other. Isolation of marine lake populations may be the result of strong physical isolating barriers as well as different selective environmental regimes in the lakes. The subterranean channels that connect each lake with the surrounding sea may provide a formidable dispersal barrier for propagules. Alternatively, a propagule may be able to enter but may not be able to survive due to the environmental regime within the lake or competition with resident founder lineages/species. The pattern that is almost consistently seen in a variety of taxa (jellyfish, fish, and bivalves) is that each lake harbors a single lineage per taxon (Dawson, 2005; Dawson & Hamner, 2005; Gotoh et al., 2009; Goto, Tamate & Hanzawa, 2011). An explanation for this pattern may be fierce inter-lineage competition, as would be expected for closely related species that occupy similar niches. This is further supported by the observation that lakes either contain mussels or oysters as the dominant bivalves, rarely both (Becking et al., 2011). It could be that priority effects (Alford & Wilbur, 1985; Wilbur & Alford, 1985) provide competitive advantage to whichever lineage colonizes a marine lake first The dominance of a lineage would be largely stochastic, likely dictated by the chance arrival of propagules (Gillespie & Baldwin, 2009): whichever mussel lineage arrives first is the one that fills that niche space. Paleoecological records (e.g., sediment cores with bivalve shells) could potentially verify which morphotypes have been dominant historically.

Each lake appears to show the genetic signature of a recent bottleneck in population size followed by an expansion such as would be expected in founder populations (Slatkin & Hudson, 1991; Rogers & Harpending, 1992). We assume that the floodwaters that filled the lakes during the Holocene sea level rise allowed for independent, chance colonization of lakes by propagules from the surrounding sea and that these propagules were the progenitors of the present-day populations (Dawson & Hamner, 2005; Dawson, 2006). The deep divergences that are observed between the lineages in the mussels are probably ancient lineages that have taken refuge in the lakes. We suggest, however, that the within lineage diversity that is unique to each lake may have resulted from in situ divergence during the 6,000–12,000 years that the lakes have existed. This scenario would imply relatively rapid evolutionary rates, but such rates are not uncommon for recently diverged taxa (e.g., Genner et al., 2007; Ho et al., 2011) and have been recorded generally in mtDNA of mussels (e.g., Lee & Foighil, 2004; Trovant et al., 2015). Combined effects of stochastic processes (e.g., founder effects), local adaptation and increased evolutionary rates could produce high levels of differentiation in small isolated populations such as in marine lake environments (Dawson & Hamner, 2005; Ho et al., 2011). The patterns of genetic variation found so far in marine lake populations of Mastigias papua (Dawson & Hamner, 2005) Brachidontes spp (Goto, Tamate & Hanzawa, 2011; this study), Suberites diversicolor (Becking et al., 2013), Sphaeramia orbicularis (Gotoh et al., 2009) are generally consistent with taxa evolving in isolation in peripatric environments, such as islands or satellite lakes of ancient rift lakes in Africa (e.g., Genner et al., 2007; Emerson & Gillespie, 2008; Chen & He, 2009; Rosindell & Phillimore, 2011). The role of marine lakes in supporting endemism in the IAA region may be through enhanced survival of endemics (relicts from earlier anchialine systems), with the possibility of local population divergence that may lead to speciation.

The islands of Samama, Kakaban, Maratua are part of the Berau Marine Protected Area that was formed in 2005 and Kakaban island with the large marine lake (Lake 1) is planned to be proposed for a UNESCO World Heritage Site (Gunawan & Visser, 2012; Vermaat, Estradivari & Becking, 2012). The implication of our results for marine conservation is that each lake is barely connected by effective dispersal and must be viewed as unique and isolated entities.