Strong genetic differentiation in tropical seagrass Enhalus acoroides (Hydrocharitaceae) at the Indo-Malay Archipelago revealed by microsatellite DNA

Study area and sampling

A total of 202 E. acoroides samples from seven localities from Java and Sumatra were collected in 2014 (Fig. 1, Table 1). At each location, 18–42 individuals were taken in a zigzag pattern along a line transect. To avoid collection of the same genet or clone, only one shoot was collected within a diameter of minimum of 5 m (DA Willette, pers. comm., 2013). Collected shoots were rinsed with fresh water to remove epiphytic algae. A young leaf from each shoot was desiccated with silica gel and preserved at room temperature until use.

DNA extraction and amplification

Silica gel-dried leaves (5 cm in length) from each shoot were ground using a mortar and pestle. Genomic DNA was extracted using DNeasy plant mini kit (Qiagen^®, Hilden, Germany) following the manufacturer’s protocol. Eight microsatellite loci (Eaco_001, Eaco_009, Eaco_019, Eaco_050, Eaco_051, Eaco_052, Eaco_054, Eaco_055) developed by Nakajima et al. (2012) were used to score individual genotypes (Table 2). Forward primers were each labeled with a fluorescent dye 6-FAM, NED, PET or VIC (Applied Biosystems^®, Foster City, CA, USA) (Table 2).

Polymerase Chain Reaction (PCR) was performed into two ways. (i) five loci (Eaco_001, Eaco_009, Eaco_019, Eaco_051, Eaco_054) were amplified using the Qiagen multiplex PCR Kit in a 10 µl total reaction containing 3 µl ddH₂O, 5 µl of 2× PCR Master Mix, 1 µl of 0.2 µM primer mix and 1 µl template DNA. PCR cycling was carried out for 5 min at 95 °C, followed by 35 cycles of 30 s at 95 °C, 1.5 min at 57 °C and 30 s at 72 °C with an extension of 30 min at 60 °C in the final cycle. (ii)three loci (Eaco_050, Eaco_052, Eaco_055) were amplified in a 20 µl total reaction containing 7.8 µlddH₂O, 2 µl of 1× PCR gold buffer, 2 µl MgCl₂ (2.5 mM), 2 µl dNTPs (200 µm), 1.5 µl of each 0.75 µm primer, 1U taq polymerase (Applied Biosystems^®, Foster City, CA, USA) and 3 µl DNA template. PCR cycling was carried out for 15 min at 95 °C, followed by 32 cycles of 30 s at 94 °C, 1.5 min at 58 °C and 60 s at 72 °C with an extension of 30 min at 60 °C in the final cycle. All PCRs were performed on a 2720 Thermal Cycler (Applied Biosystems^®, Foster City, CA, USA). GeneScan™ 500 LIZ^® (Applied Biosystems^®, Foster City, CA, USA) was used as an internal lane standard and PCR products were sent to UC Berkeley DNA Sequencing Facility, USA, where they were run on an ABI 3130 Xl automated DNA sequencer (Applied Biosystems^®, Foster City, CA, USA). Individual genotypes were scored using Geneious ver. 7.0.6 (Kearse et al., 2012).

Data analysis

Genetic diversity and Hardy-Weinberg equilibrium

Two loci (Eaco_052 and Eaco_050) were discarded after initial examination because of quality control issues. Eaco_052 had only a single allele in all populations (monomorphic), while Eaco_050 was not successfully amplified in most samples. Thus, only six loci (Eaco_001, Eaco_009, Eaco_019, Eaco_051, Eaco_054, and Eaco_055) were used for further analysis (Table 3). All samples were successfully amplified for all six loci, except a single sample from ANS for which Eaco_051 could not be amplified.

A total of 89 alleles were detected across these six microsatellite loci, ranging from one allele at the locus Eaco_019 in BM, KJ, BK, and ANS populations to 13 alleles at locus Eaco_054 in the BM population. The mean number of alleles per locus ranged from 1.57 to 9.57, and the mean number of alleles per population ranged from 4.33 to 7.00 (Table 2). The BM population had the highest average number of alleles, while the lowest was found in the KJ population. Observed heterozygosities ranged from 0.434 to 0.615 and expected (H_E) from 0.458 to 0.605, respectively. The TD population had the highest H_E value (0.605), while the ANS showed the lowest H_E (0.458).

Deviation from HWE is evidence that simple genetic models may not be appropriate (Waples & Allendorf, 2015). Hardy-Weinberg equilibrium tests revealed that five loci (Eaco_001 in TD and NK, Eaco_009 in KJ and ANS, Eaco_051 in KJ, Eaco_054 in TD, NK, ANS and Eaco_055 in ANS) deviated significantly from Hardy-Weinberg expectations (P < 0.05) prior to a multiple test correction with a false discovery rate. After correction for multiple tests, only two loci significantly deviated from Hardy-Weinberg equilibrium (Eaco_001 and Eaco_054).

Microchecker detected homozygote excess at locus Eaco_054 in TD and ANS and at locus Eaco_001 in TD and NK, likely a result of null alleles. Null alleles may also explain the significant deviation from Hardy-Weinberg expectations at these loci (Eaco_054 and Eaco_001), although inbreeding due to population substructure (Wahlund effect) and natural selection are also possible. Null alleles were potentially implicated in 16 out of 42 locus-population combinations, with an estimated frequency of <0.2 (Table 4). Null alleles may cause an overestimation of genetic differentiation (Chapuis & Estoup, 2007), but the effect is not considerable when their frequency is lower than 0.2 (Carlsson, 2008). To test the effects of null alleles to our datasets, we conducted pairwise F_ST tests by including or excluding locus Eaco_054 (in which three out of seven locations were deviated from HWE after FDR correction); removing Eaco_054 did not change the results of pairwise comparisons, therefore all loci were included in further analysis. Individual genotypes for these loci are reported in Data S1.

Population structure

Pairwise F_ST values ranged from 0.127 to 0.359. F_ST were significant between all pairs of samples with P < 0.001 (Table 5). The largest genetic differentiation was found between samples from NK and ANS (F_ST = 0.359), while the smallest pairwise difference was found between TD and PR and between BM and BK (F_ST = 0.127).

The Neighbor Joining method based on D_A distance identified three major clusters (Fig. 2). Cluster 1 (NK) was genetically distinct from all other populations. Cluster 2 consists of two populations in Java (TD and PR) and two populations in Sumatra (BK and BM). The highest bootstrap support was found between TD and PR. In contrast, two populations, KJ and ANS showed discordance between genetic and geographic distance (cluster 3).

The ΔK test in Structure (Pritchard, Stephens & Donnelly, 2000) indicated the maximum value of ΔK at K = 3 with a secondary peak at K = 6 and K = 7 (Fig. S1). At K = 3, the Structure analysis mirrors the groupings seen in Neighbor Joining tree, with the same three major clusters (NK; KJ, ANS; and TD, PR, BM, BK) (Fig. 2). Isolation-by-distance (IBD) revealed significant correlation between genetic differentiation and geographic distance across all pairs of samples, with P = 0.003 (Fig. S2). Due to the existence of hierarchical structuring, a partial Mantel test was conducted in which population NK was excluded, given the scale of its genetic and geographic distance from all other sites. A partial Mantel performed on the remaining Java and Natuna Sea locations indicated that IBD is present across these sites (p = 0.03, Fig. 3)

Genetic diversity

Genetic diversity is closely related to survivorship, resistance, and resilience of any individual or population to disturbance (Williams, 2001; Hughes & Stachowicz, 2004). Low genetic diversity can reduce the fitness of individuals, and lead to population extinction (Spielman, Brook & Frankham, 2004). In recent years, genetic diversity plays a central role in predicting the ability of individuals or populations to survive environmental change (Hughes & Stachowicz, 2004), and global climate change (Ehlers, Worm & Reusch, 2008).

Our study revealed that the highest genetic diversity was observed in the central part of Indonesia (BK and PR), while the northernmost population (ANS) possessed lower observed heterozygosity. Sample collection at BK and PR was conducted in undisturbed area where the seagrasses growth in a dense and continues bed (ING Putra & HH Madduppa, pers. comm., 2014). Seagrasses in undisturbed areas tend to be more genetically diverse (Larkin et al., 2006). Although ANS showed the lowest observed heterozygosity, this may not related to their local environment because sampling location was also conducted in undisturbed area. Microchecker results indicated that ANS had an excess of homozygosity, a potential indicator of inbreeding. Inbreeding mostly occurs in clonal organisms and reduces diversity by reducing the effective frequency of recombination throughout the genome (Charlesworth, 2003).

Previous studies of E. acoroides in Lembongan (Bali) and Waigeo (Papua), found that observed heterozygosities were 0.436–0.582, respectively (Pharmawati et al., 2015). Likewise, using nine microsatellite loci, Nakajima et al. (2014) found observed heterozygosities of E. acoroides ranging from 0.165–0.575 at three locations (Japan, China, and Philippines). Observed heterozygosities of E. acoroides varied from 0.100 to 0.567 in China (Gao et al., 2012), highlighting that our observed values are quite similar to those reported from other sites. The differences in reported heterozygosity ranges between studies could arise due to various factors such as natural selection, variation in the mutational rate, and the number and nature of the loci used.

Seagrasses in the Indo-Malay Archipelago, as in other parts in the world, remain subject to overexploitation with conservation efforts largely marginalized. Although the potential decrease of genetic diversity in seagrasses is quite high due to meadow fragmentation (Vermaat et al., 2004) and extensive clonal growth (Procaccini, Olsen & Reusch, 2007), E. acoroides appears to retain considerable genetic variation through sexual reproduction. Sexual reproduction contributes significantly to maintaining genetic variation in clonal organisms (Rollon, Cayabyab & Fortes, 2001). The available life history information also suggests that most seagrass species are predominantly outcrossing (Reusch, 2001).

Population genetic structure

Pairwise F_ST showed significant genetic differentiation between all pairs of populations. Moreover, Isolation-by-distance (IBD) showed a significant correlation between genetic differentiation and geographic distance, even when the most distant and differentiated population (NK) is excluded. These results indicate that there is low gene flow between populations, possibly due to limited dispersal potential of E. acoroides. For example, Lacap et al. (2002) showed that pollen dispersal of E. acoroides is limited over a small spatial scale, usually less than 10 km. The low gene flow discovered in the present study and from previous studies (e.g., Nakajima et al., 2014) are consistent with direct observations made on fruit and seed floating time, where median flotation periods were 7 days for fruit and 0.5 h for seeds, during which time they can cover distances of 41 km and 0.1 km, respectively (Lacap et al., 2002). After buoyancy is lost for seeds, further seed dispersal probably continues for 2 to 5 days, but the maximum distance was only 204 cm (Lacap et al., 2002).

Limited dispersal distances are also indicated in other seagrass species such as Thalassia hemprichii (Lacap et al., 2002) and Thalassia testudinum (Davis, Childers & Kuhn, 1999). Likewise, Hosokawa et al. (2015) found that seeds of Zostera marina do not disperse far from the parent bed. Thus, such local genetic differentiation may be common among seagrass species, although Zostera noltii shows a large geographic area of apparent genetic homogeneity (Coyer et al., 2004).

Phylogeographic relationships among E. acoroides suggest its partitioning into three major groups, consistent with Bayesian clustering analysis (Fig. 2). Clustering output indicated secondary peaks at K = 6 and K = 7. These results indicate the presence of further substructure similar to that found in F_ST analysis. At K = 6, all populations are separated except for populations in TD and PR, while at K = 7, all populations seem to be genetically distinct. NK, the most geographically distant population, showed distinct genetic structure from the six other sites, and represents the western Indonesian population (Indian Ocean population). Differentiation of the western Indonesian (Indian population) is common in other organisms and is believed to be largely a result of Pleistocene vicariance (Barber et al., 2000; Crandall et al., 2008; Yasuda et al., 2009). Low sea levels during the Pleistocene would result in the emergence of the Sunda Shelf, which presumably isolates E. acoroides by preventing dispersal. Studies on marine plants such as mangroves, Bruguiera gymnorrhiza (Urashi et al., 2013), genus Ceriops (Huang et al., 2008), and the seagrass H. ovalis (Nguyen et al., 2014) all showed the Sunda Shelf as being a geographical barrier between the Indian and Pacific Ocean populations. Further, similar genetic patterns consistent with Pleistocene isolation have also been reported for marine animals such as in false clownfish Amphiprion ocellaris (Nelson et al., 2000), mantis shrimp Haptosquilla pulchella (Barber et al., 2002), the coral reef fish Caesio cuning (Ackiss et al., 2013), the giant clam Tridacna crocea (DeBoer et al., 2014), and five commercially important pelagic fishes (Auxis thazard, Euthynnus affinis, Katsuwonus pelamis, Rastrelliger kanagurta, and Scomberomorus commerson) (Jackson et al., 2014) respectively.

The geographical partitioning of genetic groups observed in the present study may be a result of various events in the past, but the effect of recent physical oceanography cannot be ignored (DeBoer et al., 2014). The importance of currents in shaping the genetic structure have been reported on both large and small geographical scales (Barber et al., 2002; Galarza et al., 2009; Yasuda et al., 2009; Zhan et al., 2009; Wee et al., 2014). In this region, the South Java Current (SJC) is a unidirectional current, flowing down the west coast of Sumatra (Fig. 1). This current probably acts as a barrier to the dispersal of E. acoroides between NK and the other sites.

The genetic structuring of cluster 2 in this study mirrors the genetic structure found in the genus Hippocampus (Lourie, Green & Vincent, 2005) and Nerita albicilla (Crandall et al., 2008) respectively. However, no phylogeographic break was found across Java Sea, which is in contrast to previous studies carried out in stomatopods (Barber et al., 2002; Barber, Erdmann & Palumbi, 2006). It is possible that E. acoroides found appropriate habitat across the Seribu, Bangka, and Batam Islands and used these islands as stepping stones for long distance dispersal. Alternatively, chance may play a role such that each species responds to history in a slight different manner, resulting in a variety of genetic patterns (Bowen et al., 2016). If this were the case, we would expect a shared pattern of genetic structure to emerge if enough species were sampled, as has been reported recently in Hawaii (Toonen et al., 2011; Selkoe et al., 2014; Selkoe et al., 2016).

The alternating monsoon in Indonesia also likely plays an important role in driving dispersal in the region, and is a possible confounding factor for estimating gene flow. Alternating monsoons generate two Seasonally Reversing Currents (SRC) across the Java, Natuna, and South China Seas (Fig. 1). These large storms and reversing currents could also create a complexity of local oceanographic features. Although, previous research did not find evidence of any relationship between genetic patterns and seasonal currents in the Sunda Shelf region for clownfish (Nelson et al., 2000), our results presented here indicates that SRC may play an important role in the dispersal of E. acoroides.

Finally, although we found significant isolation-by-distance overall, two sites (KJ and ANS) showed high genetic similarity despite being geographically distant from one another. This genetic pattern may also have been influenced by sea level changes during the Pleistocene. When the sea level was at 30 m Below Present Level (BPL), land bridges connecting Sumatra and Borneo via Bangka Belitung seem to have been lost (Voris, 2000). Thus, both the Java and Natuna Seas became connected and this connection may have facilitated gene flow between ANS and KJ. At more recent timescales, oceanographic processes such as SRC may also create and maintain genetic structure between KJ and ANS. Seagrass populations in the Natuna Sea could be originating from a different source. Kool et al. (2011) suggest that areas in the Natuna Sea are supplied by larval transport from the South China Sea. Other studies showed that populations in the Natuna Sea may also receive propagules from the Andaman Sea via the Malaka Strait following the ice retreat after LGM (Timm & Kochzius, 2008; Wee et al., 2014). Due to observed differentiation between populations in BM and ANS, we suspect these populations might originate from different areas. The existence of biophysical barriers separating western and eastern part of the Natuna Sea might also responsible for the population differentiation in this area (Treml et al., 2015). However, in order to better understand the origins of the populations in BM and ANS, additional sampling sites from the Malaka Strait and the South China Sea are necessary.

Implications for restoration and management

Seagrass beds provide important habitat for a wide variety of organisms in marine ecosystems (Hemminga & Duarte, 2000). However, they are declining globally due to human activity and climate change (Duarte, 2002; Waycott et al., 2009). Restoration of seagrass beds is underway in many sites worldwide (Van Katwijk et al., 1998; Park & Lee, 2010). Recently, it has been found that genetic diversity plays an important role for seagrass resilience and restoration (Ehlers, Worm & Reusch, 2008; Reynolds, McGlathery & Waycott, 2012).

Here we show there is limited connectivity and significant genetic differentiation between all sampling sites. This study also revealed there are three different cluster based on genetic distance and Bayesian clustering. Therefore, we recommend resource managers treat populations across our sampling region as three Management Units (MUs) to maintain the unique genetic characteristics of each region to support seagrass conservation goals. Any restoration activities conducted using local stocks (the same structure cluster) would be beneficial because their genetic diversities have been assessed. Population with higher genetic diversity increases the fitness and reproductive success of the transplants (Williams, 2001). Genetic diversity also contributed to the resistance of seagrass beds to various disturbances (Hughes & Stachowicz, 2004); thus, high diversity creates more resilient ecosystems, and provides more ecosystem services (Reynolds, McGlathery & Waycott, 2012).

Introducing new genotypes through the transplanting seed or plants to other sites should be avoided. Introduction of new alleles that are distinct from native populations would lead to lower survivorship if not locally adapted (Williams, 2001). When migrants are translocated or transferred from one population to another, they may mate with native population and produce offspring that are less fit (outbreeding depression). Although, this has not been documented or recorded in seagrasses so far, this event has been found in wild populations of scarlet gilia (Ipomopsis aggregata) (Waser, Price & Shaw, 2000).