Genetic Structure and Population Demographic History of a Widespread Mangrove Plant Xylocarpus granatum J. Koenig across the Indo-West Pacific Region

Xylocarpus granatum J. Koenig is one of the most widespread core component species of mangrove forests in the Indo-West Pacific (IWP) region, and as such is suitable for examining how genetic structure is generated across spatiotemporal scales. We evaluated the genetic structure of this species using maternally inherited chloroplast (cp) and bi-parentally inherited nuclear DNA markers, with samples collected across the species range. Both cp and nuclear DNA showed generally similar patterns, revealing three genetic groups in the Indian Ocean, South China Sea (with Palau), and Oceania, respectively. The genetic diversity of the Oceania group was significantly lower, and the level of population differentiation within the Oceania group was significantly higher, than in the South China Sea group. These results revealed that in addition to the Malay Peninsula—a common land barrier for mangroves—there is a genetic barrier in an oceanic region of the West Pacific that prevents gene flow among populations. Moreover, demographic inference suggested that these patterns were generated in relation to sea level changes during the last glacial period and the emergence of Sahul Shelf which lied northwest of Australia. We propose that the three genetic groups should be considered independent conservation units, and that the Oceania group has a higher conservation priority.


Introduction
Mangrove is a unique forest ecosystem that connects different habitats: land, freshwater, and ocean. Mangrove forests form the center of mangrove ecosystems, and are distributed in tropical and subtropical coastal and riverine regions of the world [1]. The ecological services and economic contribution of mangrove forests are important to many tropical and subtropical countries and provide tremendous benefits to human beings [1,2]. However, because mangroves forests are distributed close to areas of high anthropogenic disturbance and land use change, they are rapidly declining and under destruction, which results in the severe fragmentation of forests [3,4]. Thus, it is well recognized that mangrove forests have a high priority in conservation biology and forest management in many countries, not only due to human activities, but also because of recent global warming [5].
Despite the wide distribution range of mangrove forests, the number of component species is relatively small compared with other tropical forest habitats [6]. Only 30-40 species are recognized as core mangrove component species [1,2]. Even though this may be a small number, most species are widespread, enabling mangrove forests to be distributed globally across the tropics. Core groups of mangroves, such as Avicennia spp., Rhizhophora spp., Bruguiera spp., and Xylocarpus spp., are found across different continents and oceanic regions [1,3]. Each group is an important component of mangrove ecosystems, and contributes to forest formation and succession processes. One of the mechanisms by which the core mangrove species can maintain their wide distribution is due to their effective method of seed dispersal, that is, sea-dispersal. By dispersing buoyant propagules via ocean currents, mangrove species have established wide distribution ranges spanning over continents. However, the specific processes involved in the formation and maintenance of the wide distribution of core mangrove species are not well understood.
Recently, molecular studies have suggested the presence of genetically distinct units within widespread mangrove species that are formed by barriers to gene flow. A number of studies have been performed in the Indo-West Pacific (IWP) region where a higher level of mangrove species diversity exists, e.g., Ceriops tagal (Perr.) C.B.Rob. [7], Bruguiera gymnorhiza (L.) Lam. [8,9], and Rhizophora spp. [10][11][12]. These studies showed strong genetic structure in the IWP, especially across the Malay Peninsula, which is a significant biogeographic barrier that has shaped genetic structure between the Indian and Pacific oceans. An additional "cryptic barrier" [12] in the West Pacific is reported in a study of Rhizophora stylosa Griff. Distinct genetic structure between the South China Sea and Oceania has also been found in Ceriops tagal [13] and Bruguiera gymnorhiza [14]; this concordance in phylogeography may indicate shared biogeographic histories. To understand the genetic structure of mangroves across the entire IWP region, it is necessary to focus on species that are widespread across the region. Moreover, although coalescent-based population demographic inferences have started to provide a deeper understanding of the genetic diversity and structure of species and their history (e.g., [15][16][17]), this kind of approach has not yet been well conducted in mangrove species. Thus, further inferences of past demographic history of mangrove species can shed new light on not only the population genetics of mangrove species, but also mangrove forest management and conservation, as demonstrated in forest tree species [16,[18][19][20].
Xylocarpus granatum J. Koenig is one of the most widespread core mangrove species in the IWP, and can provide an ideal system that allows us to test the presence of genetic structures across the region. The genus Xylocarpus is the only mangrove genus in the large mahogany family, Meliaceae. The genus occurs solely in the IWP region, from east Africa westward to Oceania, and is comprised of three accepted species, X. granatum, X. moluccensis M.Roem., and X. rumphii (Kostel.) Mabb. [21].
Xylocarpus species have small, white flowers that are most likely insect-pollinated. Unlike the iconic mangrove species in the Rhizophoraceae that disperse viviparous propagules, Xylocarpus produces large, cannonball-like, multi-seeded fruits, and their seeds are buoyant for dispersal via the ocean. Within the genus, X. granatum is the most widespread species, but the genetic diversity and structure of this species has not been examined. Therefore, this species can be a suitable model to evaluate in order to provide a deeper understanding of the genetic structure and past demographic history of the mangrove species across the IWP region.
The aims of this study are: (1) to examine the genetic structure of X. granatum over a wide distribution range across the IWP using maternally inherited cpDNA and bi-parentally inherited nuclear DNA markers; (2) to infer the past demographic history of the species to further understand any genetic structure and barriers detected; and (3) to propose conservation units to maintain local genetic diversity.

Sample Collection and DNA Extraction
Leaf samples of 111 individuals of X. granatum from 18 populations were collected, covering 14 countries across the entire distribution range of X. granatum in the Indo-West Pacific (Table 1, Figure 1). Total DNA was extracted from silica-dried leaves using a modified cetyltrimethylammonium bromide (CTAB) method [22]. The extracted DNA was eluted in Tris-EDTA Buffer (TE) and stored at −20 • C until sequencing and genotyping. Xylocarpus moluccensis and X. rumphii were also examined in this study to evaluate whether cpDNA variation was shared among closely related species. For X. moluccensis, which is also widespread in SE Asia, 17 individuals from five populations were examined.   Table 1. The inset (bottom) denotes the haplotype network of the four haplotypes extracted from Figure 2.

Chloroplast DNA Sequencing and Data Analysis
Three cpDNA intergenic regions, trnD-trnT [23], trnL-trnF [24], and accD-psaI [25], were sequenced for all of the samples using published universal primer pairs. Each polymerase chain reaction (PCR) was performed in a total volume of 10 μL with 1.0 μL of 10× Ex Taq buffer (Takara Bio Inc., Shiga Prefecture, Japan), 0.8 μL of each 10 mM deoxynucleotide (dNTP) mixture (Takara Bio Inc., Shiga Prefecture, Japan), 0.4 μL of each 10 pM primer pair, 0.05 μL of TaKaRa Ex Taq (Takara Bio Inc., Shiga Prefecture, Japan), and 6.35 μL of sterilized water mixed with 1.0 μL of template DNA. The PCR protocol used was initial denaturation at 95 °C for 60 s, followed by 35 cycles of 95 °C for 45 s, 55 °C for 45 s, and 72 °C for 60 s, and a final elongation at 72 °C for 10 min. The PCR products were purified using ExoSAP-IT (USB Corporation, Cleveland, OH, USA) and cycle sequencing was performed on the purified product using the BigDye cycle sequencing kit ver. 3.1 (Applied Biosystems, Foster City, CA, USA). Subsequently, the nucleotide sequence was determined with an ABI3500 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). Raw sequence data was imported into Auto Assembler ver. 2.1.1 (Applied Biosystems, Foster City, CA, USA), and peak calling was visually corrected. The sequence data of the three regions were concatenated and aligned using the ClustalW algorithm [26] incorporated within MEGA5 [27]. The haplotype of each individual was identified using the final concatenated sequence by considering only base substitution and indels. A statistical parsimony network was constructed, using the final concatenated sequences of 1702 bp by TCS 1.21 [28] to visualize the relationships among the cpDNA haplotypes. Genetic diversity (h) was calculated for each population. Genetic differentiation among populations was evaluated by calculating FST [29] and its standardized value, F'ST, which always ranges from 0 to 1 [30], using GenAlEx 6.5 (hereafter, GenAlEx [31]).

Chloroplast DNA Sequencing and Data Analysis
Three cpDNA intergenic regions, trnD-trnT [23], trnL-trnF [24], and accD-psaI [25], were sequenced for all of the samples using published universal primer pairs. Each polymerase chain reaction (PCR) was performed in a total volume of 10 µL with 1.0 µL of 10× Ex Taq buffer (Takara Bio Inc., Shiga Prefecture, Japan), 0.8 µL of each 10 mM deoxynucleotide (dNTP) mixture (Takara Bio Inc., Shiga Prefecture, Japan), 0.4 µL of each 10 pM primer pair, 0.05 µL of TaKaRa Ex Taq (Takara Bio Inc., Shiga Prefecture, Japan), and 6.35 µL of sterilized water mixed with 1.0 µL of template DNA. The PCR protocol used was initial denaturation at 95 • C for 60 s, followed by 35 cycles of 95 • C for 45 s, 55 • C for 45 s, and 72 • C for 60 s, and a final elongation at 72 • C for 10 min. The PCR products were purified using ExoSAP-IT (USB Corporation, Cleveland, OH, USA) and cycle sequencing was performed on the purified product using the BigDye cycle sequencing kit ver. 3.1 (Applied Biosystems, Foster City, CA, USA). Subsequently, the nucleotide sequence was determined with an ABI3500 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). Raw sequence data was imported into Auto Assembler ver. 2.1.1 (Applied Biosystems, Foster City, CA, USA), and peak calling was visually corrected. The sequence data of the three regions were concatenated and aligned using the ClustalW algorithm [26] incorporated within MEGA5 [27]. The haplotype of each individual was identified using the final concatenated sequence by considering only base substitution and indels. A statistical parsimony network was constructed, using the final concatenated sequences of 1702 bp by TCS 1.21 [28] to visualize the relationships among the cpDNA haplotypes. Genetic diversity (h) was calculated for each population. Genetic differentiation among populations was evaluated by calculating F ST [29] and its standardized value, F' ST , which always ranges from 0 to 1 [30], using GenAlEx 6.5 (hereafter, GenAlEx [31]).

Nuclear SSR Genotyping and Data Analysis of X. Granatum
We obtained Short Sequence Repeat (SSR) data for a total of 13 populations of X. granatum ( Table 1). The SSR data of nine populations for 11 loci was obtained from Tomizawa et al. [32] ( Table 1). The remaining four populations-5-Phuket, 14-Airai, 17-Baie de Tare, and 18-Daintree-were genotyped with the eleven SSR loci based on the PCR and fragment analysis protocol in Tomizawa et al. [32]. Allele calling was done using GeneMapper ver. 4.1 (Applied Biosystems, Foster City, CA, USA). The apportionment of genetic variation among group (F RT ), among populations within a group (F SR ), and among individuals within populations (F ST ) were evaluated based on a Bayesian model-based clustrering analysis results of K = 2 (Indian Ocean-South China Sea, and Oceania groups) and K = 3 (Indian Ocean, South China Sea, and Oceania groups) (see results for details) using an analysis of molecular variance (AMOVA, [33]) implemented in GenAlEx. The AMOVA was also conducted for each pair of the three groups at K = 3. The standardized values of F' RT and F ST were also calculated using GenAlEx. The genetic relationships among populations were evaluated by generating a neighbor-joining (NJ) tree based on the Nei's genetic distances (DA) [34] using Populations 1.2.30 beta software [35]. The statistical confidence in the topology of the tree was evaluated by 1000 bootstraps derived using the same software ( Figure S1). The NJ tree was reconstructed on a topographic map using Mapmaker and GenGIS2 software [36]. Individual-level genetic relationships were analyzed with principal coordinates analysis (PCoA) using GenAlEx. Individual-based genetic structure was evaluated using a Bayesian model-based clustering algorithm implemented in the software STRUCTURE v. 2.3.4 [37,38]. STRUCTURE analysis was performed using the admixture model [38], with 20 runs for each number of subpopulations (K), from K = 1 to 13. Each run consisted of 30,000 replicates via the Markov chain Monte Carlo (MCMC) method after a burn-in of 20,000 replicates. STRUCTURE HARVESTER [39] was employed to evaluate the probability of the data (LnP (D)) for each K, and to calculate ∆K according to the method described by Evanno et al. [40]. Multimodality among runs and major clustering patterns at each K were evaluated using the CLUMPAK server [41].

Inference of Demographic History
The software DIYABC v2.0 [42,43] was used to analyze the nuclear SSR data of X. granatum to infer demographic history based on the approximate Bayesian computation (ABC) approach. To keep the scenarios in the ABC analysis simple, three groups were defined based on the result of the STRUCTURE analysis: Pop I (Pop IDs 2, 4, 5 and 6), Pop II (Pop IDs 7,8,11,13 and 14), and Pop III (Pop IDs 15, 16, 17 and 18). The main aims were: 1) to infer population size change through time for each group, and 2) to infer the order of splitting into the three groups with their time scale, to understand how the modern genetic structure of X. granatum, as shown in the NJ tree and the STRUCTURE analysis, was generated (Figure 3a,b). Thus, we conducted two analyses: 'ABC1-Inference of population size change in each regional group', and 'ABC2-Inference of population demographic history among regional groups'. In all of the scenarios, t# represented the time scale, measured by generation time, and N# represented the effective population size of the corresponding populations (PopI, II, and III) during the relevant time period (e.g., 0-t1, t1-t2, and t2-t3).
ABC1-Inference of population size change in each regional group. Four simple scenarios were tested for each group (Figure 4a). The scenarios were as follows: ABC2-Inference of population demographic history among regional groups. Six simple population demographic scenarios were tested with considering all the possible hierarchical population splits and populations to be traced back to the common ancestral population ( Figure 4b). As population shrinkage was detected in Pop III in ABC1 (see results), population shrinkage for Pop III at t1 was added to all the scenarios, giving Na and N3 for the effective population size before and after the population size reduction, respectively. The change in the number of repeats followed a generalized stepwise mutation model (GSM; [44]) and single nucleotide indels (SNI) were also allowed. The mutation rate of the former was assumed to be higher than the mutation rate of the latter. The default values of the priors were used for all of the parameters (Table S1) except t2, according to our preliminary test runs (data not shown). The mean values of the expected heterozygosity (H E ), the number of alleles, the size of variance among the alleles, and Garza and Williamson's M [45] were used as summary statistics for each of the three populations in ABC1 and 2. As for summary statistics for each of the population pairs in ABC2, we used H E , A, Garza and Williamson's M, F ST , genotype likelihood, shared allele distance, and (δµ) 2 distance [46]. One million simulations were performed for each scenario, and the most likely scenario was evaluated by comparing posterior probabilities with the logistic regression method. The goodness-of-fit of the scenarios was also assessed by principal component analysis (PCA) using the option "model checking" in DIYABC. Six simple population demographic scenarios were tested with considering all the possible hierarchical population splits and populations to be traced back to the common ancestral population (Figure 4b). As population shrinkage was detected in Pop III in ABC1 (see results), population shrinkage for Pop III at t1 was added to all the scenarios, giving Na and N3 for the effective population size before and after the population size reduction, respectively. The change in the number of repeats followed a generalized stepwise mutation model (GSM; [44]) and single nucleotide indels (SNI) were also allowed. The mutation rate of the former was assumed to be higher than the mutation rate of the latter. The default values of the priors were used for all of the parameters (Table S1) except t2, according to our preliminary test runs (data not shown). The mean values of the expected heterozygosity (HE), the number of alleles, the size of variance among the alleles, and Garza and Williamson's M [45] were used as summary statistics for each of the three populations in ABC1 and 2. As for summary statistics for each of the population pairs in ABC2, we used HE, A, Garza and Williamson's M, FST, genotype likelihood, shared allele distance, and (δμ) 2 distance [46]. One million simulations were performed for each scenario, and the most likely scenario was evaluated by comparing posterior probabilities with the logistic regression method. The goodness-of-fit of the scenarios was also assessed by principal component analysis (PCA) using the option "model checking" in DIYABC.

Phylogenetic Relationships of the cpDNA Haplotypes in Xylocarpus
A total of 1702 bp of cpDNA sequence was used for analysis, including 687 bp of trnD-trnT Intergenic Spacer (IGS), 435 bp of trnL-trnF IGS, and 580 bp of accD-psaI. In X. granatum, four haplotypes (XG1-XG4) were characterized based on four nucleotide substitutions; while in X. moluccensis, five haplotypes (XM1-XM5) were characterized based on six nucleotide substitutions and one indel ( Table 2). The two species did not share any haplotypes, but the haplotype of X. rumphii used in this study was the same as a haplotype of X. moluccensis (XM4).
The statistical parsimony network constructed by TCS showed that X. granatum and X. moluccensis are clearly differentiated. In X. granatum, haplotype XG2 and XG4 were derived from the major haplotype XG3, and haplotype XG1 was derived from XG2 ( Figure 2). There were only one or two steps between neighboring haplotypes. In X. moluccensis and X. rumphii, nearby haplotypes were differentiated by one to three steps.

Phylogenetic Relationships of the cpDNA Haplotypes in Xylocarpus
A total of 1702 bp of cpDNA sequence was used for analysis, including 687 bp of trnD-trnT Intergenic Spacer (IGS), 435 bp of trnL-trnF IGS, and 580 bp of accD-psaI. In X. granatum, four haplotypes (XG1-XG4) were characterized based on four nucleotide substitutions; while in X. moluccensis, five haplotypes (XM1-XM5) were characterized based on six nucleotide substitutions and one indel ( Table 2). The two species did not share any haplotypes, but the haplotype of X. rumphii used in this study was the same as a haplotype of X. moluccensis (XM4).
The statistical parsimony network constructed by TCS showed that X. granatum and X. moluccensis are clearly differentiated. In X. granatum, haplotype XG2 and XG4 were derived from the major haplotype XG3, and haplotype XG1 was derived from XG2 ( Figure 2). There were only one or two steps between neighboring haplotypes. In X. moluccensis and X. rumphii, nearby haplotypes were differentiated by one to three steps.

Geographical Distribution of cpDNA Haplotypes
The geographical distribution of the haplotypes for the 18 populations of X. granatum revealed that haplotype XG3 was widely distributed across the Indo-West Pacific (Table 1, Figure 1). The distribution of haplotypes demonstrated no clear genetic break between the Indian and Pacific Oceans (due to the presence of XG3 in both). Instead, we observed a localized distribution in the other three haplotypes, with XG1 and XG2 largely confined to the Andaman Sea (4-Ayeyarwady, 5-Phuket, 6-Klang, and 12-Dong Rui) and in Oceania (14-Airai, 15-Kosrae, 16-Malatie, 17-Baie de Tare, and 18-Daintree), respectively, while XG4 was found in only one individual from Mozambique (1-Maputo). Haplotypes XG1 and XG2 were fixed in one population (4-Ayeyarwady), and four populations (15-Kosrae, 16-Malatie, 17-Baie de Tare, and 18-Daintree), respectively. The values of F ST and F' ST calculated among populations by cpDNA data were 0.853 and 0.933, respectively, suggesting the level of population differentiation was high.

Genetic Structure Revealed by SSR Markers
The total number of alleles detected was five to 20 per locus. Although the highest ∆K was detected at K = 2, the Ln P(D) values increased together with the number of K up to around K = 9 in the STRUCTURE analysis ( Figure S2), and genetic clustering patterns were clear from K = 2 to 9. At K = 2, the geographic pattern of the two clusters was clear; cluster 1 (C1 in blue color) was distributed in the Indian Ocean and South China Sea, while cluster 2 (C2 in orange color) was distributed in Oceania, consisting of 15-Kosrae (Micronesia), 16-Malatie (Vanuatu), 17-Baie de Tare (New Caledonia), and 18-Daintree (Australia). One exception was the 14-Airai (Palau) population, which clustered into C1 despite its location in Oceania (Figure 3a). At K = 3, C1 at K = 2 was further divided into two clusters, which corresponded to the Indian Ocean (C3 in dark blue color), consisting of 2-Quelimane (Mozambique), 4-Ayeyarwady (Myanmar), 5-Phuket (Thailand), and 6-Klang (Malaysia); and the South China Sea, consisting of 7-Kemaman (Malaysia), 8-Sabah (Malaysia), 11-Bali (Indonesia), 13-Panay (Philippines), and 14-Airai (Palau). A similar pattern to the results of K = 3 was also detected by the NJ tree for 13 populations (Figure 3b). Although more local genetic clusters were detected up to K = 9, the NJ tree for nine clusters revealed the same three genetic groups detected at K = 3 and the NJ tree for 13 populations. An individual level PCoA analysis demonstrated three well-segregated groups in the plot of the second two axes (1 vs. 3 explained 26.56% of the total variation) (Figure 3c) and the first two axes (1 vs. 2 explained 27.92%). Thus, finally, the 13 populations were clustered into three groups: Pop I, II, and III, prior to conducting ABC analyses (as stated in the Materials and Methods). In concordance with the cpDNA data, a genetic break was detected at the boundary between the South China Sea and Oceania. The cpDNA haplotypes from both regions were detected in population 14-Airai, which is located between the two regions, and this population also revealed a small amount of admixture between the Indian Ocean-South China Sea cluster and the Oceania cluster in the nuclear DNA data.
The overall value of F ST was 0.404, suggesting that the level of population differentiation among the 13 populations was high, especially when considering the standardized F' ST value of 0.802 (Table 3). The value of F ST among the four populations in Oceania was significantly higher than among the five populations in the South China Sea, while significant differences were neither detected among the three groups nor between other pairs of groups. There was relatively high genetic differentiation among the populations within groups (F' SR values were 0.678 and 0.610 for K = 2 and 3, respectively), although larger differentiation was detected among groups (F' RT values were 0.748 and 0.648 for K = 2 and 3, respectively) ( Table 3). For all of the comparisons, F ST values are always two to three times higher, which can be due to higher variation among individuals than among populations or groups. Moreover, the AMOVA for each group pair suggested that genetic differentiation between the South China Sea and Oceania (F' RT = 0.815) was much higher than that between the Indian Ocean and the South China Sea (F' RT = 0.389), and the Indian Ocean and Oceania (F' RT = 0.416).

Inferences of Demographic History
When considering the changes in the temporal population size in each regional group in ABC1, scenario 1 (constant population size model) showed the highest posterior probability in Pop1 (Indian Ocean) and 2 (South East Asia), and summary statistics and the PCA suggested that this scenario best described the data (Table 4, Figures S3b and S4b). However, the 95% confidence interval (CI) of scenario 1 overlapped with other scenarios in these two groups (Table 4). Regarding Pop3 (Oceania), scenario 2 (population shrinkage model) showed the highest posterior probability (0.3917) and its 95% CI (0.3858-0.3977) did not overlap with other scenarios. In this scenario, the effective population size after the shrinkage was estimated to be 2310 (95% hyper probability density (HPD), 806-6960), and the population shrinkage was inferred to have occurred 12800 (95% HPD, 888-29,200) generations ago (Table S3, Figure S5a). This time scale can be translated to 128,000 years Before Present (BP) (95% HPD, 8880-292,000) under the assumption of a generation time of 10 years for X. granatum. The median values of the mean mutation rate of SSR, mean P (the parameter of the geometric distribution to generate multiple stepwise mutations) and SNI were estimated to be 4.53 × 10 −4 (95% CI: 1.30 × 10 −4 -9.65 × 10 −4 ), 0.254 (95% CI: 0.128-0.300) and 2.54 × 10 −7 (95% CI: 1.18 × 10 −8 -7.72 × 10 −6 ), respectively. The summary statistics and the PCA suggested the good fit of the scenario (Table S4, Figure S5b). All four summary statistics in this scenario were not significantly differentiated from the observed value, and the PCA showed that the observed data point (yellow) was in the middle of a small cluster of datasets from the posterior predictive distribution (red, Figure S5b), suggesting scenario 2 fits the observed data well in Pop3. Table 4. Posterior probability of each scenario in ABC1 and its 95% hyper probability density (HPD) based on the logistic estimate by DIYABC v2.0 [42,43]. The scenarios tested in ABC1 are shown in Figure 4. In the inference of population demographic history among regional groups (ABC2), the highest posterior probability was detected for scenario 1 (0.3794, 95% CI, 0.3712-0.3877) without overlapping with the 95% CI of other scenarios (Table 5). In this scenario, the median values of N1, N2, N3, and Na were estimated as 4090 (95% HPD, 1680-7610), 7880 (95% HPD, 4140-9850), 3090 (95% HPD, 822-8510) and 7840 (95% HPD, 3740-9910), respectively. However, Na was not well estimated (Table S5, Figure  S6a). The median values of the time scales for the population shrinkage of Pop3 (t1), the splitting of Pop1 and Pop2 (t2), and the splitting of Pop1 and Pop3 (t3) were 4810 (95% HPD, 191-25,400), 4990 (95% HPD, 1460-18,500), and 15500 (95% HPD, 4590-44,300) generations ago, respectively, corresponding to 48,100 (95% HPD, 1910-254,000), 49,900 (95% HPD, 14,600-185,000), and 155,000 (95% HPD, 45,900-443,000) years BP, respectively, although t1 was not well estimated. The median values of the mean mutation rate of SSR, mean P (the parameter of the geometric distribution to generate multiple stepwise mutations), and SNI were estimated to be 3.47 × 10 −4 (95% CI: 1.58 × 10 −4 -7.82 × 10 −4 ), 0.248 (95% CI: 0.134-0.300), and 1.13 × 10 −7 (95% CI: 1.11 × 10 −8 -2.89 × 10 −6 ), respectively. All 36 summary statistics for this scenario were not significantly differentiated from the observed value, and the PCA suggested a good fit of the scenario (Table S6, Figure S6b). Moreover, the PCA showed that the observed data point (yellow) was in the middle of a small cluster of datasets from the posterior predictive distribution (light green, Figure S6b), suggesting that scenario 1 fits the observed data well. Table 5. Posterior probability of each scenario in ABC2 and its 95% hyper probability density (HPD) based on the logistic estimate by DIYABC v2.0 [42,43]. The scenarios tested in ABC2 are shown in Figure 4. This study detected clear genetic structure with three genetic clusters across the distribution of Xylocarpus granatum, namely the Indian Ocean, the South China Sea, and Oceania clusters, in both cp and nuclear DNA datasets. In particular, a genetic break with high genetic differentiation (F' RT = 0.815) was detected between the South China Sea and Oceania in the West Pacific, suggesting limited gene flow by sea-dispersed propagules between these two regions. This genetic break could be recognized as the "cryptic barrier" [12] that might have historically prevented gene flow between South East Asia and Oceania, as shown in Rhizophora stylosa [12], Bruguiera gymnorhiza [14], Sonneratia alba Griff. [47], and Vigna marina (Burm.) Merr. [48].
Although it is difficult to compare the exact location of genetic breaks among these mangrove species due to studies having different sampling locations, the presence of genetic breaks in a similar oceanic region is surprising, because all of the core mangrove species are sea-dispersed plants that can utilize oceans as corridors for migration. Sea dispersal is considered to be one of the most effective modes of seed dispersal for land plants, with dispersal ranges reaching over 100 km [49]. Indeed, long-distance dispersal by sea dispersal is well supported in population genetic studies in extremely widespread plants, especially for the pantropical plants with sea drifted seeds that have vast distribution ranges across the tropics, such as Hibiscus tiliaceus L. [50,51], Ipomoea pes-caprae (L.) R. Br. subsp. brasiliensis (L.) van Ooststr. [52] and the genus Rhizophora [10]. However, topographic barriers can act as genetic barriers even in these species. For example, Takayama et al. [10] detected the Central American Isthmus as a genetic barrier dividing the Pacific and Atlantic groups of Rhizophora mangle L. and Rhizophora racemose G. Mey. in both cp and nuclear DNA, while extended gene flow was suggested within those groups.

Inferences of the Demographic History of X. granatum
Our results suggest a clear genetic break between the South China Sea and Oceania, despite no obvious topographic barrier in this oceanic area. The ABC inference of past demographic history is helpful to try and understand this pattern. The inferred divergence time of 155,000 (95% HPD, 45,900-443,000) years BP between the South China Sea and Oceania groups suggests that the shallow sea between the Australian continent and New Guinea (the Torres Strait) might have acted as a significant land geographic and genetic barrier during and since the last glacial maximum (LGM, ca. 20,000 years BP), when the entire Sahul Shelf was exposed [53]. The past influence of the Sahul Shelf, compounded by the present large geographic distance among oceanic islands, may have resulted in the observed genetic differentiation of multiple core mangrove species over the West Pacific region. Furthermore, the westward Equatorial Current that bifurcates north and south as it approaches the Indo-Malay Archipelago [54] could be a strong barrier preventing the gene dispersal of haplotype XG2 beyond Palau. Thus, although we need to be cautious about uncertainties in the ABC inferences, such as the generation time of species, overlapping of generations, and wide 95% CI of inferred parameters in the assumed model [16], the present results suggest that genetic differentiation caused by a genetic barrier between the South China Sea and Oceania was generated by episodes of eustatic change in sea levels in relation to ice ages over the past several hundred thousand years. Moreover, population shrinkage was detected in the Oceania group (estimated 128,000 years BP (95% HPD, 8880-292,000)), suggesting that effective population size was reduced following genetic divergence from the South China Sea. In addition, as the mean value of F ST among populations within Oceania (0.4090) was significantly higher than within the South China Sea (0.1800), higher genetic isolation within the Oceania group might have also contributed to a decrease in the regional effective population size in Oceania. Regarding the assumed model in the ABC, DIYABC does not consider gene flow after divergence, and it may bias the inferred temporal parameters [16,20]. However, as genetic differentiation between the South China Sea and Oceania is high enough with less gene flow (F' RT = 0.815), this bias may be limited. In addition, since no consideration of gene flow after divergence may underestimate divergence times [16], the main discussion here would not be changed with the population split still occurring before the Last Glacial Maximum (LGM). Finally, although the ABC approach in this study revealed the demographic history of the species, the results obtained in this study might be based on the limited sample size in this study. Thus, further analysis would be required with more loci and more sample sizes in future studies in order to obtain a deeper understanding of the demography of the species.
The population in Palau (14-Airai) revealed a unique genetic composition, and the cpDNA showed a partial mixture of two haplotypes (XG2 and 3), which were dominant in the South China Sea and Oceania, respectively. In addition, nuclear DNA revealed that this population was clustered with the South China Sea group, although it is located in Oceania. As Palau is geographically located close to the boundary between the South China Sea and Oceania, a mixture of cpDNA haplotypes might be reasonable when seed flow between them is considered. Indeed, a mixture of cpDNA halotypes from different lineages has been commonly detected around the boundary areas of lineages in many plant populations (e.g., [48,55,56]). However, nuclear DNA showed only a small amount of admixture between the South China Sea and Oceania clusters, and was basically grouped with the South China Sea cluster. Thus, although the limited sample size examined in this study should be considered, this pattern could be due to differential gene flow between pollen and seeds between the South China Sea and Oceania groups, or incomplete sorting of cpDNA haplotypes. These scenarios need to be further tested with more extended sample collections that cover the species distribution.

Malay Peninsula as a Common Barrier to Widespread Mangrove Plants
The results of the NJ tree, PCoA, and the STRUCTURE analysis of nuclear SSR data showed clear genetic structure between the Indian Ocean and South China Sea, suggesting that the Malay Peninsula also acts as a barrier for gene flow in X. granatum. The Malay Peninsula has been reported as a clear land barrier for several mangrove species, such as B. gymnorhiza [8,9,14] and Ceriops tagal [7]. This land barrier emerged in this region during the LGM, and has been routinely evoked as an explanation for the genetic structure of mangrove plants across the Malay Peninsula [57]. Indeed, the inferred divergence time of 49,900 (95% HPD, 14,600-185,000) years BP in the ABC generally supports this hypothesis. Although cpDNA generally showed similar patterns to the nuclear DNA, cpDNA did not suggest clear genetic structure across the Malay Peninsula, especially when considering the presence of the widespread haplotype, XG3, from East Africa to Palau (Figure 1). The presence of the widespread haplotype over this wide range implies either frequent long-distance migration, or less variation (resolution) of the molecular marker.

Relationships with Other Xylocarpus Species
Similar to other mangrove groups, Xylocarpus granatum has its closely related species, X. moluccensis, in a core group of mangroves that has overlapping distribution with X. granatum [21]. Another species, X. rumphii, is also recognized as a mangrove, but is not very common, and the distribution is scattered. In these results, X. granatum and X. moluccensis formed two clear clades, and X. rumphii shared the same haplotype with X. moluccensis ( Figure 2, Table 2). Perhaps it is necessary to further study the taxonomic status of X. rumphii with extensive sampling from its distribution range. Only cpDNA data was collected for X. moluccensis, and the haplotype distribution map shows genetic structure across the Malay Peninsula ( Figure S7), which also suggests that the Peninsula plays a role as a barrier for the seed dispersal of mangroves, as discussed above. It is interesting that no evidence of hybridization between X. granatum and X. moluccensis was detected, from field observations or with molecular data, despite collecting samples from the same area and population in some cases (Philippines and Malaysia). In general, core mangrove species tend to have a sister mangrove species, e.g., Bruguiera gymnorhiza-B. sexangula, Rhizophora stylosa-R. mucronata, and Rhizophora mangle-R. racemosa [58]. The distribution ranges of the sister species often overlap, and when they grow in the same location, hybrids are sometimes formed. Hybrid formation is quite common, especially for the mangrove species of Rhizophoraceae. Hybrid formation can be an important phenomena for plants to increase variation and promote speciation, but may not occur in Xylocarpus.

Conclusions
In this study, a clear genetic structure was detected within the distribution range of one of the most widespread and core component species of mangrove forests in the IWP. The genetic structure was shaped not only by the Malay Peninsula, a common land barrier for mangroves, but also by a cryptic barrier across the West Pacific. Moreover, the demographic inferences revealed that genetic divergence across the West Pacific was older than that across the Malay Peninsula, implying that a common factor such as the Sahul Shelf acted to shape the genetic structure of sea-dispersed mangrove plants. Other factors that generated genetic structure also need to be considered and evaluated in future studies, including adaptation to the local environment, ancient ocean currents, the location of refugia during the ice ages, and expansion routes during oscillation cycles. However, although further study would be needed, this study has provided tentative information on conservation genetics from a neutral genetic variation point of view, and the three groups detected in this study should be treated as independent conservation units in order to maintain local genetic diversity. In addition, conservation priority should be given to the Oceania region, as genetic diversity is lower and genetic isolation among populations is more prominent than in other areas.
Supplementary Materials: The following are available online at www.mdpi.com/1999-4907/8/12/480/s1, Figure S1: Neighbor-joining (NJ) tree for 18 populations of Xylocarpus granatum based on nuclear SSR dataset. Genetic distances were calculated by Populations 1.2.30. Statistical confidence of the topology was evaluated by 1000 bootstrap replicates. Figure S2: Plot of LnP(D) and ∆K calculated for K ranging from 1 to 13. See text for detail. Figure S3: (a) The prior and posterior distributions for each parameter and (b) Principal Component Analysis (PCA) of "model checking" for the scenario 1 in ABC1 obtained for Pop1 using DIYABC. Priors and posterior probability densities are shown in the Y-axis., Figure S4: (a) The prior and posterior distributions for each parameter and (b) Principal Component Analysis (PCA) of "model checking" for the scenario 1 in ABC1 obtained for Pop2 using DIYABC. Priors and posterior probability densities are shown in the Y-axis., Figure S5: (a) The prior and posterior distributions for each parameter and (b) Principal Component Analysis (PCA) of "model checking" for the scenario 2 in ABC1 obtained for Pop3 using DIYABC. Priors and posterior probability densities are shown in the Y-axis. Figure S6: (a) The prior and posterior distributions for each parameter and (b) Principal Component Analysis (PCA) of "model checking" for the scenario 1 in ABC2 using DIYABC. Priors and posterior probability densities are shown in the Y-axis. Figure S7: Geographical distributions of the five cpDNA haplotypes of Xylocarpus moluccensis. The size of the pie charts indicates sample size of the population. Population ID are defined in Table 1. The inset (bottom) denotes the haplotype network of XM1-5 extracted from Figure 2. Table S1: Prior distributions of the parameters used in ABC1, Table S2: Prior distributions of the parameters used in ABC2, Table S3: Demographic parameters of the most likely scenario in each region in ABC1 obtained by DIYABC, Table S4: Comparison of summary statistics for the observed data set and posterior simulated data sets in the most likely scenario in each region in ABC1, Table S5: Demographic parameters of the scenario 1 in ABC2 obtained by DIYABC, Table S6: Comparison of summary statistics for the observed data set and posterior simulated data sets in the scenario 2 in ABC2.