Ancient Geographical Barriers Drive Differentiation among Sonneratia caseolaris Populations and Recent Divergence from S. lanceolata

Glacial vicariance is thought to influence population dynamics and speciation of many marine organisms. Mangroves, a plant group inhabiting intertidal zones, were also profoundly influenced by Pleistocene glaciations. In this study, we investigated phylogeographic patterns of a widespread mangrove species Sonneratia caseolaris and a narrowly distributed, closely related species S. lanceolata to infer their divergence histories and related it to historical geological events. We sequenced two chloroplast fragments and five nuclear genes for one population of S. lanceolata and 12 populations of S. caseolaris across the Indo-West Pacific (IWP) region to evaluate genetic differentiation and divergence time among them. Phylogenetic analysis based on sequences of nuclear ribosomal internal transcribed spacer and a nuclear gene rpl9 for all Sonneratia species indicate that S. lanceolata individuals are nested within S. caseolaris. We found strong genetic structure among geographic regions (South China Sea, the Indian Ocean, and eastern Australia) inhabited by S. caseolaris. We estimated that divergence between the Indo-Malesia and Australasia populations occurred 4.035 million years ago (MYA), prior to the onset of Pleistocene. BARRIERS analysis suggested that complex geographic features in the IWP region had largely shaped the phylogeographic patterns of S. caseolaris. Furthermore, haplotype analyses provided convincing evidence for secondary contact of the South China Sea and the Indian Ocean lineages at the Indo-Pacific boundary. Demographic history inference under isolation and migration (IM) model detected substantial gene flow from the Sri Lanka populations to the populations in the Java Island. Moreover, multi-locus sequence analysis indicated that S. lanceolata was most closely related to the Indian Ocean populations of S. caseolaris and the divergence time between them was 2.057 MYA, coinciding with the onset of the Pleistocene glaciation. Our results suggest that geographic isolation driven by the Pleistocene ice age resulted in the recent origin of S. lanceolata.

Introduction

Pleistocene glaciations have profoundly shaped the biogeography of extant global flora and fauna (Hewitt, 2000, 2001). While tropical animals and plants were less influenced than their temperate counterparts by temperature fluctuations during this period, dramatic sea level changes that accompanied glaciation affected organisms inhabiting tropical marine or intertidal zone (Barber et al., 2000; Kittiwattanawong et al., 2001; Timm and Kochzius, 2008; Macieira et al., 2015). Indo-Australian Archipelago (IAA), which extends from the Malay Peninsula through Sumatra, Borneo and Java to Sulawesi and the island of New Guinea, has the richest biodiversity on the planet (Mora et al., 2003). Molecular phylogenetic and population genetic studies in many marine organisms have revealed strong genetic discontinuities between the bodies of water separated by IAA (McMillan and Palumbi, 1995; Williams and Benzie, 1998; Duda and Palumbi, 1999; Nelson et al., 2000; Kochzius et al., 2003; Froukh and Kochzius, 2008). Barber et al. (2000) proposed a marine Wallace’s line according to the genetic breaks between populations from north and south of Flores and Java Seas, and suggested that Pleistocene glacial isolation acted as a driver of intra-species genetic differentiation.

Mangroves are the dominant woody plants in the intertidal zones of the tropics (Tomlinson, 1986). Like other organisms inhabiting tropical marine or intertidal zones, they exhibit a strong association of genetic divergence with geographical breaks (Duke, 1995; Dodd and Rafii, 2002; Triest, 2008). Mangrove species are split into two main biogeographic regions of species diversity, the Atlantic-East Pacific (AEP) and the Indo-West Pacific (IWP; Tomlinson, 1986; Duke et al., 1998). Within the IWP region, there are two biodiversity centers, Indo-Malesia (Indian subcontinent and Southeast Asia) and Australasia (from Australia and New Guinea to west Polynesia) (Ricklefs and Latham, 1993), which may be evolutionarily distinct (Daru et al., 2013). Although the majority of species are shared between the two regions (Hogarth, 1999, 2001), there are 7–9 species exclusively distributed in Australasia, while 11 are exclusive to Indo-Malesia (Duke et al., 1998). Genetic breaks between populations from the two regions have been observed in many mangrove species, such as Lumnitzera littorea (Su et al., 2007), Bruguiera gymnorrhiza (Urashi et al., 2013), and Rhizophora mucronata (Wee et al., 2015). Understanding the linkage between intra- and inter-species genetic diversity on one hand, and historical geological events on the other hand should provide important insights into the evolutionary mechanisms of biodiversity generation in mangrove plants.

Sonneratia, a typical mangrove genus widely distributed across the IWP, comprises six species and three hybrid taxa (Duke, 1984, 1994; Wang and Chen, 2002; Zhou et al., 2005; Qiu et al., 2008). Sonneratia caseolaris is a wide-spread species, ranging from Sri Lanka through the Malay Peninsula to China and Australia (Tomlinson, 1986; Duke, 2006). It grows in upstream reaches of river-dominated estuaries with low levels of salinity and a high supply of fresh water (Duke and Jackes, 1987; Duke et al., 1998). Its sibling species, S. lanceolata, has a high similarity in morphology and ecological requirement with S. caseolaris (Duke, 2006). Both species have red petals, flat and expanded calyx lobes, distinct mucronate leaf apex and usually occur in the upstream of rivers (Duke, 1994, 2006). However, their geographical distributions do not overlap (Duke and Jackes, 1987). S. lanceolata is restricted to northwestern Australia, southern New Guinea and few locations in Indonesia. In contrast, S. caseolaris is found in the Indo-Malesia region, northeastern Australia and northern New Guinea (Duke, 2006). These two species thus provide an ideal system to study the impact of Pleistocene sea-level changes on speciation of mangrove plants. However, previous intrageneric classifications of Sonneratia based on either morphological characteristics (Ko, 1985, 1993; Wang and Chen, 2002) or molecular data didn’t take S. lanceolata into account (Shi et al., 2000; Zhou et al., 2005; Yang et al., 2015), thus, phylogenetic relationship between S. lanceolata and S. caseolaris remains unclear.

In this study, we first determined the phylogenetic position of S. lanceolata in Sonneratia using sequences of the nuclear ribosomal internal transcribed spacer (nrITS) and a nuclear gene rpl9. We then sequenced two chloroplast fragments and five nuclear genes in samples from one population of S. lanceolata and 12 populations of S. caseolaris to assess the genetic diversity within S. caseolaris and divergence between it and S. lanceolata. Our results provide new insights into population dynamics and inter-species divergence of mangrove plants under the influence of sea-level fluctuations during the Pleistocene epoch.

Materials and Methods

Plant Materials

We sampled 16 individuals from a single population of S. lanceolata in Northern Territory, Australia, and five to 20 individuals from 12 populations of S. caseolaris across the IWP region. These 12 populations include four populations around the South China Sea (SCS), two at the Indo-Pacific boundary, two from Java Island, Indonesia, three from the coasts of the Indian Ocean, and one from Queensland, Australia. We collected leaf tissue from each individual and stored it in a plastic bag with silica gel for DNA extraction. The details of sampling locations and sample sizes are presented in Table 1.

TABLE 1

TABLE 1. Locations and sample sizes of the 12 Sonneratia caseolaris populations and the single S. lanceolata population used in this study.

DNA Extraction, Polymerase Chain Reaction (PCR) Amplification, and Sequencing

Total genomic DNA of each individual was extracted using the CTAB method (Doyle and Doyle, 1987). To resolve the phylogenetic position of S. lanceolata within the genus Sonneratia, nrITS and one nuclear gene rpl9 were amplified from two randomly chosen individuals from each population of S. caseolaris and S. lanceolata by polymerase chain reaction (PCR) with LA Taq DNA polymerase (Takara Bio, Inc., Shiga, Japan). Additionally, two chloroplast fragments (trn L – trn F and trn V – trn M) and five nuclear genes (rpl9, cpi, ppi, phi, and nhx2) were amplified from each individual of S. lanceolata and S. caseolaris to assess the patterns of genetic diversity between and within the two species. The primers for the nuclear genes were from Zhou et al. (2007). PCR products were purified using 2% agarose gel electrophoresis and extracted using the StarPrep Gel Extraction Kit (GeneStar Biosolutions, Co., Ltd, Beijing, China). The purified products were sequenced using the Sanger method in an ABI 3730 DNA automated sequencer with BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster city, CA, USA) using the same sets of primers as were used for PCR amplification. All the primers used in this study were presented in Supplementary Table S1.

Phylogenetic Analysis

To identify the phylogenetic position of S. lanceolata within the genus, phylogenetic analyses of all six Sonneratia species and the outgroup Duabanga grandiflora were carried out using the combined sequence data of nrITS and the nuclear gene rpl9, as only these two fragments could be successfully amplified in D. grandiflora. The partition-homogeneity test suggests that these two fragments can be combined for phylogenetic analysis in the Sonneratia genus (p = 0.624; Farris et al., 1995). The sequences from four other Sonneratia species, S. alba, S. ovata, S. griffithii, and S. apetala, and the outgroup D. grandiflora were downloaded from the NCBI database. GenBank accessions of all the sequences used for phylogenetic tree construction were listed in Supplementary Table S2. The sequences of these two fragments were aligned in CLUSTALX (Thompson et al., 1997), followed by manual adjustments in SeqMan v. 7.10 (DNAStar, London, UK). The phylogeny was reconstructed using two methods: maximum parsimony (MP) and maximum likelihood (ML). The analyses were performed in Mega v. 6.06 (Tamura et al., 2013). The MP analysis was carried out using the Tree-Bisection-Recombination (TBR) algorithm (Nei and Kumar, 2000) on random trees with 100 random addition replicates. The Maxtree parameter was set to 500. Statistical support for nodes was evaluated by bootstrap support (BS) values from 1000 replications. D. grandiflora was used to root the phylogenetic trees. Consistency index (CI) and retention index (RI) were also employed to estimate homoplasy. Indels were ignored. For ML analysis, TrN + I was selected as the appropriate nucleotide substitution model from a set of 56 models using the Akaike information criterion, as implemented in Modeltest 3.7 (Posada and Buckley, 2004). ML trees were constructed with Nearest-Neighbor-Interchange (NNI) and the confidence measures of nodes were estimated by BS values from 1000 replications.

Sequence Data Analyses

The sequences of chloroplast fragments and nuclear genes were aligned and edited in SeqMan v. 7.10. The two chloroplast fragments were concatenated for further analyses. Haplotypes of the combined chloroplast fragments and each nuclear gene were phased and validated in DnaSP v. 5.10 (Librado and Rozas, 2009). For each population, the number of segregating sites (S) and haplotypes (H), haplotype diversity (Hd), nucleotide diversity (θ_π) and DNA polymorphism (θ_w, Watterson, 1977) were calculated for the combined chloroplast fragments and each of the five nuclear genes using DnaSP v. 5.10. Haplotype networks were constructed in NETWORK v. 4.6.1.1 (Fluxus Technology, Ltd, Suffolk, UK) using the median-joining algorithm (Bandelt et al., 1999).

To assess genetic divergence between S. caseolaris and S. lanceolata, as well as among populations of S. caseolaris, we performed a Bayesian clustering analysis, implemented in STRUCTURE v. 2.3.3 (Pritchard et al., 2000), using the polymorphic sites of the five nuclear genes from all the populations of these species we sampled. The likelihood for K (ranging from 1 to 8) clusters was computed with 20 replicates per K. Each run had a burn-in of 200,000 iterations, followed by 1,000,000 Markov chain Monte Carlo (MCMC) iterations based on the admixture model and correlated allele frequencies model (Falush et al., 2003). The most likely number of clusters was estimated using the ΔK statistic (Evanno et al., 2005). Graphics were produced using DISTRUCT v.1.1 (Rosenberg, 2004).

To assess among-population sequence divergence, pairwise Kimura two parameter distances and K_XY (the average number of nucleotide differences between populations) were calculated at each locus using Mega v. 6.06 and DnaSP v. 5.10, respectively. Multidimensional scaling (MDS) analysis was performed based on the pairwise Kimura two-parameter distances in R v. 3.1.3¹. For both species, only the individuals without any missing data were included in PCA. Starting from K_XY estimates, putative geographic barriers in the IWP region were inferred for all 13 populations using Monmonier’s algorithm in the software Barrier v. 2.2 (Manni et al., 2004).

We also estimated divergence time between S. lanceolata and S. caseolaris using the BEAST software v. 1.8.2 (Drummond et al., 2012). According to the results of Bayesian clustering and MDS analyses, the CWC population was selected to represent the SCS clade of S. caseolaris, while the two populations SRK and BSK were selected for the Indian Ocean clade. For each population, the major allele of each nuclear gene was selected and concatenated with alleles from other genes to represent the genetic composition of each clade. The topology of the prior tree was chosen based on the phylogenetic relationships described in Figure 1. The estimated divergence times of S. ovata, S. alba, and S. caseolaris were set as three calibration points under Gaussian priors with means of 7.0, 9.0, and 11.0 million years and a standard deviation of 2.0 million years, respectively, based on a separate genomic analysis of mangrove species.

FIGURE 1

FIGURE 1. Maximum parsimony (MP) and maximum likelihood (ML) phylogenetic trees of Sonneratia based on the combined sequences of nrITS and one nuclear gene rpl9. S. lanceolata samples are shown in bold. Numbers up and down the branches are bootstrap values (>50%) of the MP and ML analysis, respectively. The abbreviations of the sampling location after species names are listed in Additional file 6: Supplementary Table S4.

Demographic History Estimation Using the Isolation with Migration Model

To assess the level of interregional gene flow, we estimated pairwise migration rates (2NM), as well as effective population size (Ne) and divergence time (t), among populations from the SCS, Indian Ocean and the Indo-Pacific boundary under the isolation with migration (IM) model. We used data from the five nuclear genes and the IMa2 program (Hey and Nielsen, 2007; Hey, 2010). Two populations from Hainan, China (CWC and CBA) were selected to represent the SCS lineage and three populations from Sri Lanka (SRK) and Bangladesh (BSK and BPF) were chosen for the Indian Ocean lineage, according to the results of Bayesian clustering and MDS analysis (Figure 3A; Supplementary Figure S2). Individuals from the BSK and BPF populations were merged to increase sample size. The population from Sungaibuntu, Java, Indonesia (ISG) with putative signal of genetic admixture was included to examine whether there was substantial secondary contact at the Indo-Pacific boundary. According to our Bayesian clustering analysis, the ISG population showed a closer relationship to the Indian Ocean lineage than the SCS lineage. Thus, in our prior tree, the SCS and the Indian Ocean lineages split, and then the ISG population diverged from the Indian Ocean lineage at t2 (Figure 5A). The divergence time of the two SCS populations CWC and CBA was denoted as t0, while that between the populations from Sri Lanka (SRK) and Bangladesh (BSK and BPF) was represented by t1. For each gene, the longest non-recombining block was extracted from the aligned data set using the Perl script IMgc (Woerner et al., 2007) to avoid biased parameter estimation caused by recombination. The HKY substitution model was employed in this analysis. The priori upper bounds of demographic parameters were set to include most of their posterior distributions as determined by preliminary simulations. However, when there were no good attainable peaks, the maximum cut-off of the ancestral population size was set to 8, pairwise migration rates to 8 and divergence time to 5. Based on the results of preliminary trial runs, we ran 100 chains each with a burn-in of 1^∗10⁵ iterations and 5^∗10⁶ sampling steps with a geometric heating scheme (-hfg -ha0.96 -hb0.9). Three independent runs were carried out with different, random number seeds to test for convergence. The 5^∗10⁶ well-mixed samples were thinned to obtain 100,000 draws from marginal posterior distributions of each parameter. Posterior most likely values were estimated by finding posterior modes (Nielsen and Wakeley, 2001). Due to the lack of accurately established mutation rates for either the genus Sonneratia or the family Lythraceae, we used the substitution rates estimated by the BEAST software for the IM analysis.

Results

Phylogenetic Position of S. lanceolata

The combined sequences of nrITS and the nuclear gene rpl9 were 1,734 bp long after alignment, with 135 parsimoniously informative sites. The MP analysis produced one most parsimonious tree with tree length of 370, CI of 0.911 and RI of 0.933 (Figure 1). The ML analysis yielded the same topology as the MP method. The six Sonneratia species were divided into two clades with strong support, one of which grouped S. lanceolata and S. caseolaris (BS = 98% and 78% for MP and ML analysis, respectively) and the other contained the remaining four species (BS = 99% for both MP and ML analyses). Surprisingly, S. lanceolata was nested within the clade formed by S. caseolaris samples and was most closely related to the S. caseolaris individuals from Indo-Malesia region with high bootstrap support (BS = 97% and 96% for MP and ML analysis, respectively).

Polymorphisms within S. caseolaris

We sequenced 210 individuals from 12 populations of S. caseolaris. We observed two single-nucleotide polymorphisms (SNPs) and one one-bp indel in the combined chloroplast fragments. These sites generated four haplotypes across the 12 populations of S. caseolaris (Figure 2A). Only two populations (MSB and ISG) of S. caseolaris were polymorphic with extremely low levels of nucleotide diversity (θ_π) and polymorphisms (θ_w, Supplementary Table S3).

FIGURE 2

FIGURE 2. Geographic distribution of haplotypes and Median-Joining network for the combined chloroplast fragments (A) and five nuclear genes (B: rpl9, C: cpi, D: ppi, E: phi, and F: nhx2) in the 12 populations of Sonneratia caseolaris and the single population of S. lanceolata. Each haplotype was represented by a single circle whose size reflects haplotype frequency. Haplotypes with close relationships are marked with the same color. Each branch with more than one mutational step is labeled. Population abbreviations are defined in Table 1.

The aligned sequence length of the five nuclear genes ranged from 620 to 1,170 bp, containing 12 to 49 SNPs in the 12 populations. At the species level, S. caseolaris harbored relatively high levels of genetic diversity, with θ_π ranging from 3.660 to 8.850 per kb, and θ_w from 2.820 to 5.560 per kb (Supplementary Figure S1; Supplementary Table S3). At the population level, populations of S. caseolaris from its range margins (CWC, SRK, BPF, BSK, and AQS) exhibited lower levels of polymorphism than those from its distribution center (MSB, MKK, MSN, ISG, and ICJ, Supplementary Figure S1). MSB and ISG populations had the highest levels of nucleotide diversity ( = 3.556 and 3.558 per kb, respectively), while the population from Sri Lanka (SRK) had the lowest diversity ( = 0.054 per kb).

Population Structure within S. caseolaris Across the IWP Region

Haplotype network analyses at both the chloroplast fragments and the five nuclear genes showed strong population structure that tracked with geographical distribution of S. caseolaris (Figure 2). Four chlorotypes in the 12 populations fell into three major clusters corresponding to three geographic regions. The first cluster comprises populations around the SCS (CWC, CBA, TST, MSB, MKK, and MSN), the second is restricted to the Indian Ocean (SRK, BPF, BSK, and ICJ), and the third is found exclusively in Queensland, Australia (AQS, Figure 2A). Similar to the chloroplast fragments, haplotype networks of four of the five nuclear genes (rpl9, ppi, phi, and nhx2) showed three distinct groups restricted to the SCS, the Indian Ocean and eastern Australia (Figures 2B,D,E,F). Polymorphisms in the cpi gene revealed only two groups, with one predominant in the SCS and the other combining the Indian Ocean and Australian populations (Figure 2C). The haplotype networks, in addition to a MDS analysis of the genetic distance matrix, showed that the population from eastern Australia was most highly diverged from the SCS and the Indian Ocean (Figure 2; Supplementary Figure S2). This pattern is consistent with early separation of the Indo-Malesia from Australasia.

Bayesian clustering analysis using STRUCTURE suggests that the statistically optimal number of populations is K = 2 (Supplementary Figure S3). However, the results for K = 3 and 4 were also biologically meaningful and were included in this study. At K = 3 Bayesian clustering analysis generated population assignments consistent with haplotype network and MDS analysis, splitting samples from the SCS, the Indian Ocean and eastern Australia into their own groups (Figure 3A; Supplementary Figure S2). Using the Barrier software, strongest barriers, supported by four of the five genes, were identified between the populations from the SCS and the Indian Ocean, corresponding to the high genetic divergence of the SCS and the Indian Ocean lineages (Figure 3B). Within the Indian Ocean group, a relatively weak barrier marked a subdivision between Sri Lanka and Bangladesh, which was in agreement with the Bayesian clustering result with K set to 4. These results indicated that geographic features in the IWP region might have profoundly shaped population structure within S. caseolaris.

FIGURE 3

FIGURE 3. Likely genetic clusters and geographical barriers existing among the 12 populations of S. caseolaris and the single population of S. lanceolata. (A) Bayesian clustering analysis of the five nuclear genes implemented in STRUCTURE. To illustrate the hierarchical population structure across the IWP region, we show clustering with the best-fitting model of K = 2, as well as K = 3 and K = 4. (B) Putative geographical barriers identifed within the IWP regions. The bars show the output of Bayesian clustering analysis with K = 4. Population abbreviations are defined in Table 1.

Within S. caseolaris, the divergence time between the Indo-Malesia and Australasia populations was estimated to be 4.035 million years ago (MYA; 95% of the highest probability density interval (HPDI): 2.298–6.318 MYA, Figure 4). Within the Indo-Malesia region, the split between the SCS and the Indian Ocean lineages occurred at 2.561 MYA (95% HPDI: 1.375–3.981 MYA), while the Sri Lanka population diverged from Bangladesh (BSK and BPF) 1.286 MYA (95% HPDI: 0.522–2.174 MYA).

FIGURE 4

FIGURE 4. Divergence time of the six Sonneratia species, as well as the different lineages of S. caseolaris. The scale bar is 1.0 million years (MYA). The value and purple bar at each node indicate the estimated divergence time (MYA) with its 95% highest probability density interval (HPDI).

Genetic Admixture at the Indo-Pacific Boundary

We observed extensive admixture between the SCS and the Indian Ocean groups, as evidenced by the analyses of chloroplast fragments as well as nuclear genes, in the Sangaibuntu population (ISG), located near the Indo-Pacific boundary. Similar patterns were also observed in three other populations that are close to the Malacca Strait (MSB, MKK, and MSN). Bayesian clustering likewise showed signals of genetic admixture between the SCS and the Indian Ocean groups in the ISG, MSB, and MSN populations (Figure 3A). When K was set to 2, approximately 74.6% of the nucleotide variations present in the ISG population came from the Indian Ocean lineage, while 25.4% were from the SCS lineage. In the two Malaysian populations, the proportion of SNPs that belong to the SCS lineages was slightly higher (78.1 and 87.3% in MSB and MSN, respectively).

We then used the isolation with migration (IM) model to estimate demographic parameters of S. caseolaris populations from the SCS (CWC and CBA), the Indian Ocean (BSK, BPF and SRK) and the Indo-Pacific boundary (ISG) (Figure 5; Table 2). The effective sample size (ESS) was larger than 70,000, indicating a well-mixed chain. As a result of running the BEAST program on our data, we obtained an estimate of the mutation rate at 1.673 ^∗ 10^-9 substitutions per site per year (s/s/y; 95% HPDI: 1.219 ^∗ 10^-9–2.195 ^∗ 10^-9 s/s/y). Under this assumed mutation rate, effective size (Ne) of the two SCS populations (CWC and CBA) was estimated to be 107 (95% HPDI: 0–1,253) and 322 (95% HPDI: 0–64,092), respectively; Ne of the populations from the Indian Ocean (BSK/BPF and SRK) was 537 (95% HPDI: 179–2,111) and 36 (95% HPDI: 0–680), respectively (Figure 5C). Ne of the ISG population was estimated at 1,396 (95% HPDI: 394–4,187), much greater than that in other populations. The migration rate from the Indian Ocean population SRK to the population ISG was estimated to 0.119 (95% HPDI: 0.000–0.452), which was significantly greater than zero (LLR = 3.224, p < 0.05, Figure 5D).

FIGURE 5

FIGURE 5. Probability density plots of the demographic parameters estimated from the isolation with migration model for populations from the SCS (CWC and CBA), the Indian Ocean (SRK, BSK and BPF) and the Indo-Pacific boundary region (ISG). (A) A schematic of the isolation with migration model. Population abbreviations are defined in Table 1. (B) Probability density of the divergence time between CWC and CBA (T0), between the Sri Lanka population (SRK) and two Bangladesh populations (BSK and BPF, T1) and the formation time of the ISG population (T2). (C) Probability density of the effectively population sizes (Ne) of the five populations. (D) Probability density estimation of the migration rate from the CWC, CBA, BSK and BPF and SRK to ISG populations.

TABLE 2

TABLE 2. MLEs and the 95% highest probability density interval (HPDI) of all demographic parameters estimated using the isolation with migration model (IMa2).

Genetic Divergence between S. caseolaris and S. lanceolata

In contrast to S. caseolaris, genetic diversity of the S. lanceolata population (SLA) was extremely low in both chloroplast fragments and the five nuclear genes, with θ_π ranging from 0.000 to 0.110 and θ_w from 0.000 to 0.350 per kb (Supplementary Table S3).

The position of S. lanceolata within the clade that comprises the S. caseolaris samples suggested by our phylogenetic analysis was verified by the haplotype network reconstruction and Bayesian clustering. No polymorphism was detected in the chloroplast fragments sequenced in S. lanceolata, and this species differs from the Indian Ocean populations of S. caseolaris by only one mutation (Figure 2A). At the two of the five nuclear genes (ppi and nhx2), S. lanceolata is most closely related to the Indian Ocean populations of S. caseolaris (Figures 2D,F). At the gene phi, the haplotype of S. lanceolata forms its own clade that differs from the Indian Ocean S. caseolaris samples by 11 mutations (Figure 2E). In contrast, the remaining two genes, rpl9 and cpi, contained haplotypes that grouped with the Malaysia population (MSN) of S. caseolaris (Figures 2B,C).

For the optimal clustering (K = 2, Supplementary Figure S3) of the Bayesian clustering, S. lanceolata was assigned to the same cluster as the Indian Ocean and Australian populations of S. caseolaris, distinct from the SCS group of that species (Figure 3A). When K = 3, there was a division between the S. caseolaris populations from the Indian Ocean and Australia of S. caseolaris with the S. lanceolata population falling into the Indian Ocean population cluster. Increasing K to 4 results in S. lanceolata samples forming a new cluster with the two Bangladesh populations of S. caseolaris. These results suggest that S. lanceolata is closely related to the Indian Ocean lineage of S. caseolaris. Moreover, based on the pairwise K_XY values (Supplementary Figure S4; Supplementary Table S4), two putative barriers to gene flow between S. caseolaris and S. lanceolata were identified. One is located between northern Australia and Java Island, and the other is between the Arafura Sea and the Coral Sea. These two geographic barriers surround the current distribution region of S. lanceolata (Figure 3B). This suggests that the isolation of the Arafura Sea from the surrounding oceanic regions might have triggered the S. lanceolata speciation event. Divergence time estimates indicates that S. lanceolata diverged from the Indian Ocean lineage of S. caseolaris at approximately 2.057 MYA with 95% HPDI extending from 1.053 to 3.275 MYA (Figure 4).

Discussion

Geographic Isolation Drives Strong Genetic Differentiation among S. caseolaris Populations

Strong genetic differentiation among the S. caseolaris populations was associated with geographic barriers. We identified three distinct genetic lineages located in the SCS, Indian Ocean and Australia (Figures 2 and 3A; Supplementary Figure S2). We dated the split between the Indo-Malesia and the Australasia regions to 4.035 MYA (95% HPDI: 2.298–6.318 MYA, Figure 4), prior to the Pleistocene. This result suggests a critical role for pre-Pleistocene events in genetic divergence of the two main lineages of S. caseolaris.

We further estimated that the SCS lineage diverged from the Indian Ocean lineage 2.561 MYA (95% HPDI: 1.375–3.981 MYA), in the early Pleistocene. Perhaps sea level changes during this epoch drove the genetic differentiation between populations from these geographic areas. Dispersal of seeds via sea currents plays an important role in determining the extant distribution and genetic composition of mangrove species within and among different oceanic regions (Duke, 1995; Dodd and Rafii, 2002; Triest, 2008). During glacial periods, global sea levels dropped and exposed the bulk of the shallow seabed in the SCS, forming temporary land bridges connecting mainland Asia with the three Sunda Islands (Sumatra, Borneo and Java; Voris, 2000; Sathiamurthy and Voris, 2006). This landmass likely acted as a strong physical barrier halting sea-drifted gene flow between the SCS and the Indian Ocean populations of S. caseolaris during Pleistocene, consistent with the differentiation mechanism suggested for other mangrove species (Ge and Sun, 2001; Liao et al., 2007; Huang et al., 2008; Urashi et al., 2013; Wee et al., 2015) and numerous marine organisms (Barber et al., 2000; Nelson et al., 2000; Crandall et al., 2008; DeBoer et al., 2008; Froukh and Kochzius, 2008).

It is usually difficult to distinguish genetic admixture from incomplete lineage sorting, thus reducing the resolution of signals of gene flow (Ng et al., 2015). In this study, we provide novel evidence for genetic admixture of the SCS and Indian Ocean lineages at the Indo-Pacific boundary in S. caseolaris. Analyses of haplotype frequencies at the chloroplast fragments and four of the five nuclear genes showed that the ISG population harbors haplotypes from both the SCS and the Indian Ocean groups (Figures 2B–D,F). Haplotypes of nuclear genes fell into two highly divergent clusters separated from each other by 5–14 mutational steps and there were no intermediate haplotypes. This pattern is highly unlikely to have arisen within ISG in the face of ancestral polymorphism, strongly suggesting that previously isolated lineages came into contact and mixed in this population. This idea is consistent with our simulation under the isolation with migration (IM) model that identified substantial gene flow from the Sri Lanka population (SRK) into ISG (Figure 5D).

Although the basin of the Java Sea was exposed during glaciations, especially the Last Glacial Maximum (LGM), subsequent interglacial sea-level rises restored the sea corridors connecting the SCS and the Indian Ocean (Wyrtki, 1961). It is thus easily to imagine that seeds from both the SCS and the Indian Ocean could rapidly re-colonize the seashore of the Java Sea via drift on sea currents, resulting in secondary contact of the two previously isolated lineages in this area. Similar genetic admixture of SCS and Indian Ocean groups was observed in two mangrove species Rhizophora apiculata (Yahya et al., 2014) and Lumnitzera racemosa (Li et al., 2016), as well as the false clown anemonefish (Amphiprion ocellaris, Timm and Kochzius, 2008) and reef communities (Hoeksema, 2007).

Recent Origin of S. lanceolata Triggered by Glacial Isolation

S. lanceolata was first described as a new species by Blume (1851) and is highly similar to S. caseolaris in both morphology and ecology, suggesting a close relationship between these two species (Duke and Jackes, 1987; Duke et al., 1998). Coincidently, our phylogenetic analysis using the combined sequences of nrITS and the nuclear gene rpl9 indicates that S. lanceolata is nested within the S. caseolaris clade (Figure 1). The main differentiating character between these species is stamen color, which is white in S. lanceolata (Duke and Jackes, 1987; Duke, 1994, 2006). However, the stamen colors of S. caseolaris vary in different locations. For example, the stamens of S. caseolaris are purplish red in Australia and Bangladesh, while they are white in the upper part and red in the basal part in China and Thailand (S. Shi, personal observation). Bayesian clustering analysis and K_XY statistics showed that S. lanceolata individuals (collected in Australia) were genetically closer to the Indian Ocean group than either the SCS or Australia groups of S. caseolaris (Figure 3A; Supplementary Figure S4). The split between these two species was estimated to occur 2.057 MYA (95% HPDI: 1.053–3.275 MYA, Figure 4), during early Pleistocene.

Our results of genetic similarity and divergence time are consistent with an important role of geographic isolation driven by Pleistocene glaciations as a trigger of lineage break between S. lanceolata and S. caseolaris. At the boundary between Indo-Malesia and Australasia, over 25,000 volcanic ocean islands and the shallow seas surrounding them are particularly sensitive to sea level fluctuations (Briggs, 1987, 2005; Bellwood et al., 2004; Hoeksema, 2007). During Pleistocene glaciations, sea levels dropped as much as 116 m (Miller et al., 2005; Sathiamurthy and Voris, 2006), connecting the islands and forming a strong physical barrier that likely prevented genetic exchange between the two regions (Triest, 2008; Esselstyn et al., 2009; Lohman et al., 2011). Likewise, the Torres Strait is a shallow watercourse connecting the Arafura Sea and the Coral Sea (only 7–15 m deep; Harris, 1988) that spent approximately 70% of the time since the Pleistocene closed to sea traffic (Miller et al., 2005). These two barriers to gene flow likely enabled the split between S. lanceolata and S. caseolaris.

Similar genetic differentiation and speciation might have also occurred in other mangrove genera. Avicennia integra is a narrowly distributed species confined to the Northern Territory, Australia (Duke, 2006). In contrast, its sister species A. officinalis has a much wider range extending from South India through Malaysia and Indonesia to New Guinea and eastern Australia (Tomlinson, 1986). The divergence time between these two species was estimated at 2.5 MYA, corresponding to the onset of Pleistocene (Li, 2015). Another case was observed between two sibling species of Ceriops, C. tagal and C. australis. Morphologically, C. australis differs from C. tagal by the shorter terete and smooth hypocotyls and longer flowers, calyx lobes and petals (White, 1926; Ballment et al., 1988; Sheue et al., 2009). Geographically, C. australis is restricted to Australia, along the costal lines of North Territory and Queensland, and a part of New Guinea and Indonesia, while C. tagal widely occurs from southern India through Southeast Asia to Australia (Duke, 2006). However, C. australis is found more often than C. tagal in Australia (Sheue et al., 2009). These findings suggested that Pleistocene glacial vicariance may have been a general force contributing to recent inter-species divergence among numerous mangrove genera.

Author Contributions

SS and RZ designed the study. YY, FP, and CZ collected materials and performed experiments. YY, FP, JL, and SY analyzed and interpreted the data. YY, ND, FP, RZ, and SS wrote the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by grants from the National Natural Science Foundation of China (Grant Nos. 41130208, 91331202, and 31130069); Grant 33000-18811202 from the 985 Project; the Science Foundation of State Key Laboratory of Biocontrol (grant numbers SKLBC13A03) and the Chang Hungta Science Foundation of Sun Yat-Sen University.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer YS declared a shared affiliation, though no other collaboration, with several of the authors YY, FP, JL, SY, RZ, SS to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.

Acknowledgments

We wish to thank Dr. John Yong (Singapore University of Technology and Design) for his help in sample collection and Dr. Anthony Greenberg (Bayesic Research LLC) for his help in language editing.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01618

TABLE S1 | Gene IDs, gene descriptions, primer sequences and GenBank accession numbers of haplotypes for the two chloroplast fragments and five nuclear genes used in this study.

TABLE S2 | GenBank accession numbers of the sequences of the nuclear ribosomal internal transcribed spacer (nrITS) and the nuclear gene rpl9 for the six Sonneratia species and the outgroup D. grandiflora used in this study.

TABLE S3 | Statistics of the combined chloroplast fragments and the five nuclear genes for the populations of S. caseolaris and S. lanceolata. Population IDs are defined in Table 1. S, the number of segregating sites; h, the number of haplotypes; Hd, haplotype diversity; θπ, nucleotide diversity; θw, DNA polymorphism.

TABLE S4 | Means and standard variations (SD) of pairwise K_XY values of the five nuclear genes among the 13 populations of S. caseolaris and S. lanceolata. Means were listed in lower left matrix, while SD were listed in the upper right matrix. Population abbreviations are defined in Table 1.

Footnotes

1. ^ https://www.r-project.org/

References

Ballment, E. R., Tj Iii, S., and Stoddart, J. A. (1988). Sibling species in the mangrove genus Ceriops (Rhizophoraceae), detected using biochemical genetics*. Aust. Syst. Bot. 1, 391–397. doi: 10.1071/SB9880391

CrossRef Full Text | Google Scholar

Bandelt, H.-J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Barber, P. H., Palumbi, S. R., Erdmann, M. V., and Moosa, M. K. (2000). Biogeography: a marine Wallace’s line? Nature 406, 692–693. doi: 10.1038/35021135

PubMed Abstract | CrossRef Full Text | Google Scholar

Bellwood, D. R., van Herwerden, L., and Konow, N. (2004). Evolution and biogeography of marine angelfishes (Pisces: Pomacanthidae). Mol. Phylogenet. Evol. 33, 140–155. doi: 10.1016/j.ympev.2004.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Blume, C. L. (1851). Museum Botanicum Lugduno-Batavium: Sive, Stirpium Exoticarum Novarum vel Minus Cognitarum ex Vivis aut Siccis Brevis Expsitio et Descriptio. Leiden: EJ Brill.

Google Scholar

Briggs, J. C. (1987). Biogeography and Plate Tectonics. Amsterdam: Elsevier.

Google Scholar

Briggs, J. C. (2005). The marine East Indies: diversity and speciation. J. Biogeogr. 32, 1517–1522. doi: 10.1111/j.1365-2699.2005.01266.x

CrossRef Full Text | Google Scholar

Crandall, E. D., Jones, M. E., Munoz, M. M., Akinronbi, B., Erdmann, M. V., and Barber, P. H. (2008). Comparative phylogeography of two seastars and their ectosymbionts within the Coral Triangle. Mol. Ecol. 17, 5276–5290. doi: 10.1111/j.1365-294X.2008.03995.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Daru, B. H., Yessoufou, K., Mankga, L. T., and Davies, T. J. (2013). A global trend towards the loss of evolutionarily unique species in mangrove ecosystems. PLoS ONE 8:e66686. doi: 10.1371/journal.pone.0066686

PubMed Abstract | CrossRef Full Text | Google Scholar

DeBoer, T. S., Subia, M. D., Erdmann, M. V., Kovitvongsa, K., and Barber, P. H. (2008). Phylogeography and limited genetic connectivity in the endangered boring giant clam across the Coral Triangle. Conserv. Biol. 22, 1255–1266. doi: 10.1111/j.1523-1739.2008.00983.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dodd, R. S., and Rafii, Z. A. (2002). Evolutionary genetics of mangroves: continental drift to recent climate change. Trees 16, 80–86. doi: 10.1007/s00468-001-0142-6

CrossRef Full Text | Google Scholar

Doyle, J. J., and Doyle, J. L. (1987). A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 19, 11–15.

Google Scholar

Drummond, A. J., Suchard, M. A., Xie, D., and Rambaut, A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. doi: 10.1093/molbev/mss075

PubMed Abstract | CrossRef Full Text | Google Scholar

Duda, T. F. Jr., and Palumbi, S. R. (1999). Population structure of the black tiger prawn, Penaeus monodon, among western Indian Ocean and western Pacific populations. Mar. Biol. 134, 705–710. doi: 10.1007/s002270050586

CrossRef Full Text | Google Scholar

Duke, N. C. (1984). A mangrove hybrid, Sonneratia × gulngai (Sonneratiaceae) from north-eastern Australia. Austrobaileya 25, 103–105.

Google Scholar

Duke, N. C. (1994). A mangrove hybrid Sonneratia xurama (Sonneratiaceae) from northern Australia and southern New Guinea. Aust. Syst. Bot. 7, 521–526. doi: 10.1071/SB9940521

CrossRef Full Text | Google Scholar

Duke, N. C. (1995). Genetic diversity, distributional barriers and rafting continents—more thoughts on the evolution of mangroves. Hydrobiologia 295, 167–181. doi: 10.1007/BF00029124

CrossRef Full Text | Google Scholar

Duke, N. C. (2006). Australia’s Mangroves: The Authoritative Guide to Australia’s Mangrove Plants. Brisbane, QLD: University of Queensland.

Google Scholar

Duke, N. C., Ball, M. C., and Ellison, J. C. (1998). Factors influencing biodiversity and distributional gradients in mangroves. Glob. Ecol. Biogeogr. Lett. 7, 27–47. doi: 10.2307/2997695

CrossRef Full Text | Google Scholar

Duke, N. C., and Jackes, B. R. (1987). A systematic revision of the mangrove genus Sonneratia (Sonneratiaceae) in Australasia. Blumea 32, 277–302.

Google Scholar

Esselstyn, J. A., Timm, R. M., and Brown, R. M. (2009). Do geological or climatic processes drive speciation in dynamic archipelagos? The tempo and mode of diversification in Southeast Asian shrews. Evolution 63, 2595–2610. doi: 10.1111/j.1558-5646.2009.00743.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Falush, D., Stephens, M., and Pritchard, J. K. (2003). Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587.

PubMed Abstract | Google Scholar

Farris, J. S., Kallersjo, M., Kluge, A. G., and Bult, C. (1995). Constructing a significant test for incongruence. Syst. Biol. 44, 570–572. doi: 10.1093/sysbio/44.4.570

PubMed Abstract | CrossRef Full Text | Google Scholar

Froukh, T., and Kochzius, M. (2008). Species boundaries and evolutionary lineages in the blue green damselfishes Chromis viridis and Chromis atripectoralis (Pomacentridae). J. Fish Biol. 72, 451–457. doi: 10.1111/j.1095-8649.2007.01746.x

CrossRef Full Text | Google Scholar

Ge, X.-J., and Sun, M. (2001). Population genetic structure of Ceriops tagal (Rhizophoraceae) in Thailand and China. Wetl. Ecol. Manag. 9, 213–219. doi: 10.1023/A:1011156707160

CrossRef Full Text | Google Scholar

Harris, P. T. (1988). Sediments, bedforms and bedload transport pathways on the continental shelf adjacent to Torres Strait, Australia—Papua New Guinea. Cont. Shelf Res. 8, 979–1003. doi: 10.1016/0278-4343(88)90058-1

CrossRef Full Text | Google Scholar

Hewitt, G. (2000). The genetic legacy of the Quaternary ice ages. Nature 405, 907–913. doi: 10.1038/35016000

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewitt, G. M. (2001). Speciation, hybrid zones and phylogeography—or seeing genes in space and time. Mol. Ecol. 10, 537–549. doi: 10.1046/j.1365-294x.2001.01202.x

CrossRef Full Text | Google Scholar

Hey, J. (2010). Isolation with migration models for more than two populations. Mol. Biol. Evol. 27, 905–920. doi: 10.1093/molbev/msp296

PubMed Abstract | CrossRef Full Text | Google Scholar

Hey, J., and Nielsen, R. (2007). Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc. Natl. Acad. Sci. U.S.A. 104, 2785–2790. doi: 10.1073/pnas.0611164104

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoeksema, B. W. (2007). “Delineation of the Indo-Malayan centre of maximum marine biodiversity: the Coral Triangle,” in Biogeography, Time, and Place: Distributions, Barriers, and Islands, ed. W. Renema (Berlin: Springer), 117–178.

Google Scholar

Hogarth, P. J. (1999). The Biology of Mangroves. New York, NY: Oxford University Press.

Google Scholar

Hogarth, P. J. (2001). Mangrove ecosystems. Encycl. Biodivers. 3, 853–870. doi: 10.1006/rwbd.1999.0372

CrossRef Full Text | Google Scholar

Huang, Y., Tan, F., Su, G., Deng, S., He, H., and Shi, S. (2008). Population genetic structure of three tree species in the mangrove genus Ceriops (Rhizophoraceae) from the Indo West Pacific. Genetica 133, 47–56. doi: 10.1007/s10709-007-9182-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kittiwattanawong, K., Nugranad, J., and Srisawat, T. (2001). High genetic divergence of Tridacna squamosa living at the west and east coasts of Thailand. Phuket Mar. Biol. Cent. Spec. Publication 25, 343–347.

Google Scholar

Ko, W.-C. (1985). Notes on the genus Sonneratia (Sonneratiaceae) in Southeast Asia. Acta Phytotaxonomica Sin. 23, 311–314.

Google Scholar

Ko, W.-C. (1993). Notes on some Sonneratia (Sonneratiaceae) of China. J. Trop. Subtrop. Bot. 1, 11–13.

Google Scholar

Kochzius, M., Söller, R., Khalaf, M. A., and Blohm, D. (2003). Molecular phylogeny of the lionfish genera Dendrochirus and Pterois (Scorpaenidae, Pteroinae) based on mitochondrial DNA sequences. Mol. Phylogenet. Evol. 28, 396–403. doi: 10.1016/S1055-7903(02)00444-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Yang, Y., Chen, Q., Fang, L., He, Z., Guo, W., et al. (2016). Pronounced genetic differentiation and recent secondary contact in the mangrove tree Lumnitzera racemosa revealed by population genomic analyses. Sci. Rep. 6, 29486. doi: 10.1038/srep29486

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X. (2015). The Mechanism of Speciation for Mangroves and the Evolution of Avicennia in Indo-West Pacific. Doctor’s thesis, Sun Yat-sen University, Guangzhou.

Liao, P.-C., Havanond, S., and Huang, S. (2007). Phylogeography of Ceriops tagal (Rhizophoraceae) in Southeast Asia: the land barrier of the Malay Peninsula has caused population differentiation between the Indian Ocean and South China Sea. Conserv. Genet. 8, 89–98. doi: 10.1007/s10592-006-9151-8

CrossRef Full Text | Google Scholar

Librado, P., and Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi: 10.1093/bioinformatics/btp187

PubMed Abstract | CrossRef Full Text | Google Scholar

Lohman, D. J., de Bruyn, M., Page, T., von Rintelen, K., Hall, R., Ng, P. K. L., et al. (2011). Biogeography of the Indo-Australian archipelago. Annu. Rev. Ecol. Evol. Syst. 42, 205–226. doi: 10.1146/annurev-ecolsys-102710-145001

CrossRef Full Text | Google Scholar

Macieira, R. M., Simon, T., Pimentel, C. R., and Joyeux, J.-C. (2015). Isolation and speciation of tidepool fishes as a consequence of Quaternary sea-level fluctuations. Environ. Biol. Fishes 98, 385–393. doi: 10.1007/s10641-014-0269-0

CrossRef Full Text | Google Scholar

Manni, F., Guerard, E., and Heyer, E. (2004). Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmonier’s algorithm. Hum. Biol. 76, 173–190. doi: 10.1353/hub.2004.0034

PubMed Abstract | CrossRef Full Text | Google Scholar

McMillan, W. O., and Palumbi, S. R. (1995). Concordant evolutionary patterns among Indo-West Pacific butterflyfishes. Proc. R. Soc. Lond. B Biol. Sci. 260, 229–236. doi: 10.1098/rspb.1995.0085

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, K. G., Kominz, M. A., Browning, J. V., Wright, J. D., Mountain, G. S., Katz, M. E., et al. (2005). The Phanerozoic record of global sea-level change. Science 310, 1293–1298. doi: 10.1126/science.1116412

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, C., Chittaro, P. M., Sale, P. F., Kritzer, J. P., and Ludsin, S. A. (2003). Patterns and processes in reef fish diversity. Nature 421, 933–936. doi: 10.1038/nature01393

PubMed Abstract | CrossRef Full Text | Google Scholar

Nei, M., and Kumar, S. (2000). Molecular Evolution and Phylogenetics. New York, NY: Oxford University Press.

Google Scholar

Nelson, J. S., Hoddell, R. J., Chou, L. M., Chan, W. K., and Phang, V. P. E. (2000). Phylogeographic structure of false clownfish, Amphiprion ocellaris, explained by sea level changes on the Sunda shelf. Mar. Biol. 137, 727–736. doi: 10.1007/s002270000379

CrossRef Full Text | Google Scholar

Ng, W. L., Onishi, Y., Inomata, N., Teshima, K. M., Chan, H. T., Baba, S., et al. (2015). Closely related and sympatric but not all the same: genetic variation of Indo-West Pacific Rhizophora mangroves across the Malay Peninsula. Conserv. Genet. 16, 137–150. doi: 10.1007/s10592-014-0647-3

CrossRef Full Text | Google Scholar

Nielsen, R., and Wakeley, J. (2001). Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics 158, 885–896.

PubMed Abstract | Google Scholar

Posada, D., and Buckley, T. R. (2004). Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests. Syst. Biol. 53, 793–808. doi: 10.1080/10635150490522304

PubMed Abstract | CrossRef Full Text | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959.

Google Scholar

Qiu, L., Zhou, R.-C., Li, Y.-Q., Havanond, S., Jaengjai, C., and Shi, S.-H. (2008). Molecular evidence for natural hybridization between Sonneratia alba and S-griffithii. J. Syst. Evol. 46, 391–395.

Google Scholar

Ricklefs, R. E., and Latham, R. E. (1993). “Global patterns of diversity in mangrove floras,” in Species Diversity in Ecological Communities: Historical and Geographical Perspectives. Chicago, IL: University of Chicago Press, 215–229.

Google Scholar

Rosenberg, N. A. (2004). DISTRUCT: a program for the graphical display of population structure. Mol. Ecol. Notes 4, 137–138. doi: 10.1046/j.1471-8286.2003.00566.x

CrossRef Full Text | Google Scholar

Sathiamurthy, E., and Voris, H. K. (2006). Maps of Holocene sea level transgression and submerged lakes on the Sunda Shelf. Nat. Hist. J. Chulalongkorn Univ. Suppl. 2, 1–43.

Google Scholar

Sheue, C. R., Liu, H. Y., Tsai, C. C., Rashid, S. M. A., Yong, J. W. H., and Yang, Y. P. (2009). On the morphology and molecular basis of segregation of Ceriops zippeliana and C. decandra (Rhizophoraceae) from Asia. Blumea 54, 220–227. doi: 10.3767/000651909X476193

CrossRef Full Text | Google Scholar

Shi, S., Huang, Y., Tan, F., He, X., and Boufford, D. E. (2000). Phylogenetic analysis of the Sonneratiaceae and its relationship to Lythraceae based on ITS sequences of nrDNA. J. Plant Res. 113, 253–258. doi: 10.1007/PL00013926

CrossRef Full Text | Google Scholar

Su, G., Huang, Y., Tan, F., Ni, X., Tang, T., and Shi, S. (2007). Conservation genetics of Lumnitzera littorea (Combretaceae), an endangered mangrove, from the Indo-West Pacific. Mar. Biol. 150, 321–328. doi: 10.1007/s00227-006-0357-6

CrossRef Full Text | Google Scholar

Tamura, K., Stecher, G., Peterson, D., Filipski, A., and Kumar, S. (2013). MEGA6: molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729. doi: 10.1093/molbev/mst197

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, J. D., Gibson, T. J., Plewniak, F., Jeanmougin, F., and Higgins, D. G. (1997). The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25, 4876–4882. doi: 10.1093/nar/25.24.4876

PubMed Abstract | CrossRef Full Text | Google Scholar

Timm, J., and Kochzius, M. (2008). Geological history and oceanography of the Indo-Malay Archipelago shape the genetic population structure in the false clown anemonefish (Amphiprion ocellaris). Mol. Ecol. 17, 3999–4014. doi: 10.1111/j.1365-294X.2008.03881.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomlinson, P. B. (1986). The Botany of Mangroves. Cambridge tropical biology series. Cambridge: Cambridge University Press.

Google Scholar

Triest, L. (2008). Molecular ecology and biogeography of mangrove trees towards conceptual insights on gene flow and barriers: a review. Aquat. Bot. 89, 138–154. doi: 10.1016/j.aquabot.2007.12.013

CrossRef Full Text | Google Scholar

Urashi, C., Teshima, K. M., Minobe, S., Koizumi, O., and Inomata, N. (2013). Inferences of evolutionary history of a widely distributed mangrove species, Bruguiera gymnorrhiza, in the Indo-West Pacific region. Ecol. Evol. 3, 2251–2261. doi: 10.1002/ece3.624

PubMed Abstract | CrossRef Full Text | Google Scholar

Voris, H. K. (2000). Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations. J. Biogeogr. 27, 1153–1167. doi: 10.1046/j.1365-2699.2000.00489.x

CrossRef Full Text | Google Scholar

Wang, R., and Chen, Z. (2002). Systematics and biogeography study on the family Sonneratiaceae. Guihaia 22, 214–219.

Google Scholar

Watterson, G. A. (1977). Heterosis or neutrality? Genetics 85, 789–814.

Google Scholar

Wee, A. K. S., Takayama, K., Chua, J. L., Asakawa, T., Meenakshisundaram, S. H., Adjie, B., et al. (2015). Genetic differentiation and phylogeography of partially sympatric species complex Rhizophora mucronata Lam. and R. stylosa Griff. using SSR markers. BMC Evol. Biol. 15:57. doi: 10.1186/s12862-015-0331-3

PubMed Abstract | CrossRef Full Text | Google Scholar

White, C. T. (1926). A variety of Ceriops tagal C. B Rob. (= C. candolleana W. & A.). J. Bot. 64, 220–221.

Google Scholar

Williams, S. T., and Benzie, J. A. H. (1998). Evidence of a biogeographic break between populations of a high dispersal starfish: congruent regions within the Indo-West Pacific defined by color morphs, mtDNA, and allozyme data. Evolution 52, 87–99. doi: 10.2307/2410923

CrossRef Full Text | Google Scholar

Woerner, A. E., Cox, M. P., and Hammer, M. F. (2007). Recombination-filtered genomic datasets by information maximization. Bioinformatics 23, 1851–1853. doi: 10.1093/bioinformatics/btm253

PubMed Abstract | CrossRef Full Text | Google Scholar

Wyrtki, K. (1961). Scientific Results of Marine Investigations of the South China Sea and the Gulf of Thailand 1959-1961. NAGA Report 2. San Diego, CA: The University of California, 164–169.

Google Scholar

Yahya, A. F., Hyun, J. O., Lee, J. H., Kim, Y. Y., Lee, K. M., Hong, K. N., et al. (2014). Genetic variation and population genetic structure of Rhizophora apiculata (Rhizophoraceae) in the greater Sunda Islands, Indonesia using microsatellite markers. J. Plant Res. 127, 287–297. doi: 10.1007/s10265-013-0613-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Yang, S., Fang, L., Li, J., Zhong, C., Zhou, R., et al. (2015). Phylogenetic position of Sonneratia griffithii based on sequences of the nuclear ribosomal internal transcribed spacer and 13 nuclear genes. J. Syst. Evol. 53, 47–52. doi: 10.1111/jse.12102

CrossRef Full Text | Google Scholar

Zhou, R., Shi, S., and Wu, C.-I. (2005). Molecular criteria for determining new hybrid species—an application to the Sonneratia hybrids. Mol. Phylogenet. Evol. 35, 595–601. doi: 10.1016/j.ympev.2005.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R., Zeng, K., Wu, W., Chen, X., Yang, Z., Shi, S., et al. (2007). Population genetics of speciation in nonmodel organisms: I. Ancestral polymorphism in mangroves. Mol. Biol. Evol. 24, 2746–2754. doi: 10.1093/molbev/msm209

PubMed Abstract | CrossRef Full Text | Google Scholar