To cover geographic distributions of the three Rhizophora species as wide as possible, we collected and genotyped a total 935 individuals from 54 populations across the whole IWP region (Table 1 and Figure 1). The present samples represent an expansion of those included in Lo et al. (2014), with much larger sample sizes at the population level for all three Rhizophora species. Fresh leaves were sampled from individual trees located at >10 m intervals to avoid non-random inclusion of closely related individuals within each wild population, dried with silica gel in sealed plastic bags, and preserved in an electronic auto-dry cabinet (Weifo, Taiwan).

The dried leaf pieces of each individual sample were disrupted in a 2 ml microcentrifuge tube using QIAGEN TissueLyser II (QIAGEN, Hilden, Germany), and total genomic DNA was extracted using the standard CTAB method (Doyle and Doyle, 1987). Thirteen microsatellite loci, including RM11, RM38, RM41, and RM46 (Rosero-Galindo et al., 2002), Rhst01, Rhst02, Rhst11, Rhst13, Rhst15, and Rhst19 (Islam et al., 2004, 2013), RS19, RS59, and RS67 (Takayama et al., 2009), were genotyped for each sampled individual of all populations and species. Multiplex PCRs were carried out in a final volume of 5 μL with florescent-labeled primers using QIAGEN Multiplex PCR Kit (Qiagen, Hilden, Germany) following manufacturer's recommendation for microsatellite amplification. Specifically, the reaction mix contained 1 × QIAGEN Multiplex PCR Master Mix, 0.2 μM of each primer, and less than 0.1 μg template DNA. The PCR cycling protocol was as follows: an initial activation step for 15 min at 95°C, followed by 30 cycles of 30 s denaturation at 94°C, 90 s annealing at 57°C, and 60 s extension at 72°C, and a final extension step of 30 min at 60°C. Three multiplex mixes were applied in genotyping according to the compatibility of primers and the length of PCR product facilitated by Multiplex Manager 1.2 (Holleley and Geerts, 2009). RM11 (NED), RM38 (VIC), RM41 (6-FAM), and Rhst15 (PET) formed the first combination, RM46 (VIC), Rhst01 (6-FAM), Rhst02 (NED), and Rhst11 (PET) formed the second combination, and Rhst13 (VIC), Rhst19 (6-FAM), RS19 (6-FAM), RS59 (PET), and RS67 (NED) formed the last combination. The forward or reverse primer of each primer pair was labeled with a different dye (6FAM, Integrated DNA Technologies; NED, PET, and VIC, Applied Biosystems) as indicated in the parentheses following the locus name. PCR products were analyzed with an ABI 3130xl capillary sequencer following the manufacturer's instructions, with the GeneScan 500 LIZ size standard (Applied Biosystems) as internal marker. Alleles were scored using GENEMAPPER software version 4.0 (Applied Biosystems) and were manually checked for accuracy.

Null alleles and potential stuttering were checked with MICRO-CHECKER version 2.2.3 (Van Oosterhout et al., 2004). The frequencies of null alleles were estimated with FREENA (Chapuis and Estoup, 2007) with 10,000 replications, using the implemented expectation maximization algorithm of Dempster et al. (1977). Linkage disequilibrium between each pair of loci across populations within each species was tested with FSTAT v2.9.3 (Goudet, 2001). No significant linkage disequilibrium was detected between loci in these populations.

Deviation from Hardy-Weinberg equilibrium (HWE) was tested for each of the locus—population combinations using GENEPOP (http://wbiomed.curtin.edu.au/genepop). The exact p-values of HWE tests were estimated with a Markov chain algorithm described by Guo and Thompson (1992), and the parameters of Markov chain were set as 10,000 dememorization steps and 1000 batches, with 10,000 iterations per batch. A global test of heterozygote deficiency was also carried out for each population using the same set of Markov chain parameters. The p-values of significance tests were subsequently adjusted for multiple tests with sequential Bonferroni correction (Holm, 1979; Rice, 1989).

The observed number of alleles per locus (A), rarefied number of alleles per locus (A_R), and inbreeding coefficient (F_IS) within populations were estimated using FSTAT 2.9.3 (Goudet, 2001). To reduce the bias arising from unequal population sample sizes, we used a standardized sample size of 10 to calculate A_R following the rarefication process of Petit et al. (1998). We further calculated the effective number of alleles per locus (A_E) with GenAlEx version 6.5 (Peakall and Smouse, 2012), which corrects for uneven allele frequencies at each locus. The observed heterozygosity (H_O) and expected heterozygosity (H_E) within populations, and the total gene diversity (H_T) at the species level, inbreeding coefficient within individuals (F_IS) and genetic differentiation among populations (F_ST) were also computed using GenAlEx version 6.5 (Peakall and Smouse, 2012). The multilocus outcrossing rate (t) was estimated based on the mean F_IS over loci within each population using the equation t = (1 − F_IS)/(1 + F_IS) (Wright, 1969). The estimate of t can exceed 1 or 100% if F_IS is negative, indicating heterozygote excess compared to random mating.

To analyze population genetic structure, we first carried out an allele-size permutation test to see whether mutation was non-negligible compared with genetic drift and migration (Hardy et al., 2003) before applying any F_ST estimation. We compared observed R_STvalues (before randomization) with pR_ST values obtained from 10,000 allele size permutations using SPAGeDi 1.4 (Hardy and Vekemans, 2002), and found no significant difference between them, confirming that mutation could not be the main cause of genetic differentiation among these populations. We then estimated F_ST (Weir, 1996) between all population pairs within each species and determined the significance of the estimate using FSTAT 2.9.3 (Goudet, 2001).

To further elucidate genetic relationships among populations at both intra- and interspecific levels, we performed the following three analyses. First, we conducted a discriminant analysis of principal components (DAPC, Jombart et al., 2010) for each species and for all three species combined using the R package adegenet (Jombart, 2008). We then examined genetic structure of each species and the relationships between species using the Bayesian model-based clustering method implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000). Five independent assessments were carried out for each of the three species separately, for populations in all three species combined, and for populations of R. mucronata and R. stylosa combined. To determine the optimal K-value that best fits the data for each analysis, 20 iterations were performed for each K-value, from K = 1 to K = 10. Each run consisted of 2 × 10⁶ Markov chain Monte Carlo (MCMC) iterations with a burn-in period of 5 × 10⁵, employing the admixture model with correlated allele frequencies between clusters (Falush et al., 2003). Selection of the best K-value was done based on the ΔK parameter proposed by Evanno et al. (2005) and implemented in the Structure Harvester program version 0.6.94 (Earl and Vonholdt, 2012). For comparing the implied best Ks and selected Ks, 20 iterations were carried out with 2 × 10⁶ replicates of the MCMC after a burn-in period of 2 × 10⁶. We used the program CLUMPAK (Kopelman et al., 2015) to produce graphical displays. The third analysis was to characterize genetic relationships among populations using the neighbor-joining (NJ) method based on Nei's genetic distance (Nei et al., 1983) between populations with POPTREE2 (Takezaki et al., 2010), using a bootstrapping value of 10,000.

To estimate recent gene flow (m, proportion of individuals with migrant origin) between populations and between species, we used BayesAss v3.0.3 (Wilson and Rannala, 2003). The populations of this study were sampled to cover a vast geographic area of the three species' distributional ranges, but also included those from the same coastline and located adjacent to each other. In some locations, the same neighboring populations were sampled for more than one species. These population pairs provided rare opportunities for testing whether the direction of local oceanic current affected the direction of gene flow. Thus, we considered BayesAss analysis was well suited for detecting gene flow between closely located pairs of populations for comparison with the level of gene flow between those populations situated across vast stretches of oceans or on opposite sides of large land masses, in order to provide the F_ST estimates that reflect species-wide genetic differentiation among populations. We estimated recent gene flow between populations within each species with 3 × 10⁸ iterations, a burn-in value of 3 × 10⁷, and a sampling interval of 15,000. The mixing parameters for MCMC chains were adjusted in order to reach a swap acceptance rate between 20 and 60% as suggested by Wilson and Rannala (2003). The mixing parameters for estimating migration rate, inbreeding coefficient, and allele frequencies were set to 0.4, 1.0, and 0.8 respectively. Ten replicate runs of the MCMC with different seeds were carried out to check for convergence for each analysis. Visualization of chain mixing and convergence were accomplished with TRACER v1.6 (Rambaut et al., 2014). The results presented here were from the runs with the largest effective sample sizes.

All the estimates of genetic diversity were relatively low at the population level in all three IWP Rhizophora species (Table 1), compared to the population mean values measured with microsatellite markers for other plant species (see Section Discussion below). The mean values of allelic richness were similar among the three species, such as the number of alleles per locus (A), averaged 0.315, 0.315, and 0.305 in R. apiculata, R. mucronata, and R. stylosa respectively, although the range of the estimates among populations was much broader, from 1.85 (AUS.TRI) to 4.54 (THL.MNP) in R. apiculata, from 1.92 (SRL.RKW) to 4.31 (PHP.PNL) in R. mucronata, and from 1.77 (CHN.DNZ) to 6.15 (AUS.CTR) in R. stylosa. At the species level, similar amount of total gene diversity (H_T) was found in R. apiculata (0.621), R. mucronata (0.648), and R. stylosa (0.620) (Table 2), and nearly half of this diversity resided between populations, as indicated by the highly significant F_STestimates in all the three species.

Compared to populations located in or near the centers of these species' geographic ranges, the levels of genetic diversity were often lower in peripheral populations. For example, in R. apiculata, lower genetic diversity was found in populations located in its geographically margins, such as those in East Australia (AUS.DTR and AUS.TRI), China (CHN.QLH and CHN.TLH), and Pacific islands (USA.GM, MCN.CHK, and MCN.KSR; with a rare exception of MCN.YAP, H_E = 0.466), whereas relatively high estimates of allelic richness and heterozygosity were found in South East Asia, such as populations from Indonesia (IDN.NSL and IDN.TRK), Malaysia (MLS.BLG), Philippines (PHP.PNL), and Thailand (THL.RNG and THL.PNB). In R. stylosa, high levels of heterozygosity were found in two populations from North Australia (AUS.CTR, H_E = 0.650; and AUS.EMB, H_E = 0.600), whereas much lower levels were found in populations from China, Fiji, Guam, and Kiribati. However, there were exceptions to the general pattern, such as in R. mucronata the highest level of heterozygositywas found in East Australia (AUS.TRI, H_E = 0.583), while the lowest estimate was found in a population from Thailand (THL.PNB, H_E = 0.217).

Among a total of 702 locus-population combinations (calculated based on the 13 microsatellite loci × 54 populations studied), 164 of the combinations were monomorphic as some of the loci were not variable in all 54 populations. Significant deviations from HWE were detected in 30 of the remaining 538 locus-population combinations, and the global test of heterozygosity deficiency revealed that 13 out of 54 populations significantly deviated from HWE, corresponding to the significantly positive F_IS found in four populations (CHN.QLH, MCN.YAP, SRL.PBL, and THL.MNP) of R. apiculata, four populations (AUS.TRI, IDN.TRK, MCN.YAP, and PHP.PNL) of R. mucronata, and five populations (AUS.CDR, AUS.TBC, FJ.SVP, MCN.CHK, and PHP.PNL) of R. stylosa (Table 1). However, a correlation analysis of the estimated null allele frequencies (Supplementary Table 1) and F_IS values at each of the 13 microsatellite loci in these populations indicated that the presence of null alleles at some of the loci could result in an underestimate of H_O and thus an overestimate of F_IS. We found that when null allele was not present at a locus, the estimate of F_IS was near zero or even negative, which indicates outcrossing. Because a species' mating system affects all loci equally, whereas null alleles affect only those loci where they are present, zero or negative F_IS at the loci without null alleles provide strong evidence for an outcrossing mating system in the same population. Thus, the multilocus outcrossing rate t, computed based on the mean F_IS over all loci (including those with null alleles), must represent an underestimate of outcrossing or an overestimate of inbreeding within each population. The estimates of t ranged from 0.441 to 1.567 among populations of R. apiculata, from 0.338 to 1.225 among populations of R. mucronata, and from 0.504 to 1.894 among populations of R. stylosa. The average outcrossing rate over populations was 0.821, 0.776, and 1.051 in R. apiculata, R. mucronata, and R. stylosa, respectively, indicating that they are predominantly outcrossing species, consistent with a primarily wind-pollinated mating system.

Most of the intraspecific population pairs were significantly differentiated from each other (P < 0.05 after Bonferroni correction; Supplementary Table 2). However, the pairwise F_ST values were highly heterogeneous, ranging from 0.01 (between AUS.DTR and AUS.TRI) to 0.75 (between AUS.DTR and CHN.QLH) in R. apiculata; from 0.07 (MLS.KPS and THL.PNB) to 0.72 (AUS.DTR and SRL.RKW) in R. mucronata; and from 0.03 (AUS.TDR and AUS.MTB) to 0.76 (USA.GM and FJ.SVP) in R. stylosa. As expected, the high F_ST values were found between geographically distant populations, whereas the low F_ST values were found between geographically adjacent populations. At the species level, the mean F_ST values over all 13 loci were similar, 0.461 in R. apiculata, 0.452 in R. mucronata, and 0.484 in R. stylosa (Table 2). The high levels of population genetic differentiation were accordant with the presence of private alleles in many of the populations. Private alleles were detected in 9 of 20 populations of R. apiculata, 11 of 12 populations of R. mucronata, and 14 of 22 populations of R. stylosa. The number and frequencies of private alleles in these populations are given in Supplementary Table 3.

The DAPC results showed that individuals of R. apiculata could be clearly separated from those of R. mucronata and R. stylosa without ambiguity (Figure 2). However, the clusters of R. mucronata and R. stylosa showed some degree of overlapping (Figure 2A). Three clusters were identified in R. apiculata, namely Australia, Pacific islands, and SE Asia (Figure 2B). This pattern was also detected in R. stylosa, but the separation between the three geographical groups was not as well-defined as in R. apiculata (Figure 2C). In R. mucronata, individuals from Australia and Kenya were grouped together, individuals from SE Asia and Pacific islands were grouped together, and individuals from Sri Lanka were separated from the other two clusters (Figure 2D).

The inferred best K-values from STRUCTURE cluster analysis (Figure 3; Supplementary Figure 1) were different between the species, with K = 7 clusters for R. apiculata (Figure 3A), 2 clusters for R. mucronata (Figure 3B), and 2 for R. stylosa (Figure 3C). When adopting the same K = 2 clusters for interspecific comparisons, populations of R. apiculata from Australia and Pacific Islands formed one cluster, and all populations from SE Asia formed the other (Supplementary Figure 2A). The minor mode of K = 7 in R. apiculata (supported by 7 out of 20 iterations; Supplementary Figure 2B) was similar to the major mode (12 out of 20; Figure 3A), but with all Australian populations placed in one group. In R. mucronata, the major cluster mode (10 out of 20; Figure 3B) separated most of the populations from Australia, Indonesia, and Micronesia from others, but with considerable admixture in several of the populations. However, a minor cluster mode (7 out of 20; Supplementary Figure 2C) separated the two Australian populations from all others, whereas the second minor mode (3 out of 20; Supplementary Figure 2D) placed populations from Australia, Sri Lanka and Kenya into one group, and all other populations formed the second cluster. In contrast, populations of R. stylosa from Australia and Fiji were well separated from populations of SE Asia and Pacific Islands, with little admixture.

The NJ analysis based on allele frequencies in populations of all three species showed three major groups (Figure 4), with each group consisting of largely conspecific populations. Within the R. apiculata group, populations can be further divided into two subgroups, one consisting of populations from South East Asia, and the other including populations form Australia and Pacific islands. Populations of R. mucronata were grouped together with two exceptions, AUS.TRI and AUS.DTR from Australia, which were grouped with populations of R. stylosa from the sympatric sites. The R. stylosa group was further divided into two subgroups, with populations from Pacific islands and South East Asia forming one group and all Australian populations and Fiji forming the other. Among the Australian populations, the two populations from North Australia (AUS.EMB and AUS.CTR) were clearly divergent from those of East Coast.

The estimates of gene flow/migration rate between intraspecific populations, ranged from 0.0068 to 0.1721 in R. apiculata; from 0.0088 to 0.2052 in R. mucronata, and from 0.0062 to 0.1513 in R. stylosa (Supplementary Table 4). With rare exceptions, the relatively higher rate of gene flow was detected for geographically adjacent population pairs, e.g., from AUS.TRI to AUS.DTR (0.1389) and from THL.RNG to THL.PNB (0.1721) in R. apiculata; from AUS.TRI to AUS.DTR (0.0866) and from SRL.RKW to SRL.PBL (0.1425) in R. mucronata; from AUS.EMB to AUS.CTR (0.1419) and from AUS.CDR to AUS.TDR (0.1513) in R. stylosa. Furthermore, strong directional differences in gene flow were detected between some of these adjacent populations, e.g., gene flow from AUS.EMB to AUS.CTR was much higher than that from the opposite direction in both R. apiculata (0.1285 vs. 0.0077) and R. stylosa (0.1419 vs. 0.0072); gene flow from AUS.TRI to AUS.DTR was up to 18 times higher than that from AUS.DTR to AUS.TRI and the same pattern of directionality was observed in all three species (0.1389 vs. 0.0072 in R. apiculata; 0.0888 vs. 0.0089 in R. mucronata; and 0.1243 vs. 0.0070 in R. stylosa). Strong directional differences in gene flow were detected in seven population pairs in R. apiculata, four in R. mucronata, and seven in R. stylosa (highlighted in the Supplementary Table 4), suggesting that directionality of wind and ocean currents must have played a major role in the dispersal of water-floating propagules or pollen between adjacent populations of Rhizophora.

Comparative analyses of allelic diversity, observed and expected heterozygosities revealed similarly low genetic variation at the population and species levels in the three IWP Rhizophora species, which is consistent with their common life history traits, similar environmental conditions and evolutionary history. However, in contrast to the expectation that long-lived woody species maintain more variation within species and within populations but have less variation among populations than species with other life forms (Hamrick et al., 1992), we detected an opposite pattern of population genetic diversity in the three Rhizophora species: mostly low variation within populations but high differentiation among populations. Compared to the mean values of genetic variation at the population level estimated with genomic microsatellite markers in species belonging to the same taxonomic group rosids of eudicots (A = 6.21, H_O = 0.480, and H_E = 0.579; Merritt et al., 2015), all the corresponding measures of genetic variation were much lower in the populations of IWP Rhizophora (Table 1), despite being primarily wind-pollinated, long-lived trees with ocean current-dispersed seedlings and wide geographic distributions. We considered possible sampling effects, but did not find any significant relationships between population sample size and genetic diversity parameters (correlation test p = 0.276 and 0.145 for A and H_E, respectively). Repeated founder effects or bottleneck events associated with geo-climatic history in the IWP region and genetic drift in small and geographically isolated populations are most likely the major cause of the low genetic diversity within populations and within species of Rhizophora, especially in marginal populations of the species' distribution ranges. With rare exceptions, the lowest estimates of genetic diversity were found in the peripheral populations of each species. Low genetic diversity in marginal populations is also shown in other studies of R. stylosa. For example, Wilson (2009) sampled 20 populations of R. stylosa from eastern Australia (representing the species' southern geographic range limit), and reported much lower estimates of heterozygosity (H_O = 0.19 and H_E = 0.26). Islam et al. (2013) reported a very low level of H_O (0.130) but higher H_E (0.338) in 16 populations of R. stylosa from Sakishima islands of the Japanese archipelago (representing the species' northern distributional edge). These marginal populations of Rhizophora may represent descendants of a few ancestral propagules dispersed from the more centrally located populations after the last glacial maximum (LGM), resulting in further reduction of genetic variation due to the founder effect. Low genetic diversity in peripheral populations was also detected in other mangrove tree species, such as Bruguiera gymnorhiza (Rhizophoraceae) from its northern geographic limit in Japan (Islam et al., 2014), and Rhizophora mangle and Avicennia germinans (Avicenniaceae) in the AEP region (Pil et al., 2011; Sandoval-Castro and Muñiz-Salazar, 2012; Sandoval-Castro et al., 2014).

Geographic distributions of mangrove trees are strongly influenced by air temperature and availability of intertidal habitats, with the latter being tightly linked to sea level fluctuations (Duke et al., 1998). During the Pleistocene epoch, the sea level fluctuated to a great extent, up to 125 m lower than the present level (Voris, 2000). Along with every cycle of warming and cooling, mangrove populations, especially those located in their distributional limits, would have experienced repeated local extinction and recolonization by a few founders, leading to a very small long-term effective population size (N_E), which is often a small fraction of the census population size (Woodruff, 2001; Frankham et al., 2002). The generally small N_E alone could account for the observed pattern of low variation within populations and high divergence between populations due to enhanced genetic drift in mangroves. A recent large-scale comparative genomic study has shown that mangrove populations generally have much smaller N_E and lower level of heterozygosity than non-mangrove plants based on genome-wide single nucleotide polymorphism (SNP) data (S. Shi and C-I Wu, personal communications). It can be concluded that genetic drift might have played a dominant role in the evolutionary history of mangrove trees, as evidenced at both SNP and microsatellite loci across the genome.

Compared to the present estimates, much lower levels of allelic diversity and heterozygosity in populations of R. mucronata and R. stylosa (e.g., H_O averaged 0.108 and 0.097 respectively) were reported in Wee et al. (2015). They interpreted the low heterozygosity as a result of inbreeding in the two species. However, other factors need to be taken into consideration before attributing the lack of heterozygosity to inbreeding. Furthermore, the characteristics of microsatellite markers can affect the estimates of genetic variation (Merritt et al., 2015). For comparing the level of genetic variation between species or between studies, it is necessary to use the same set of marker loci. Nevertheless, despite the difference in marker loci employed between studies, the same pattern of genetic diversity: low within- but high between-population variation, was found in R. mucronata and R. stylosa in this study and Wee et al. (2015) as well as in R. apiculata, indicating that the same evolutionary processes operated on genetic diversity in these species.

Rhizophora species are self-compatible, weakly protandrous, and predominantly wind-pollinated, though pollination by small insects may also exist (Tomlinson, 1986; Kondo et al., 1987; Nadia and Machado, 2014). According to their population outcrossing rates (t), self-compatible plant species can be further characterized as predominantly outcrossing (t > 80%), mixed mating (t = 20–80%), or predominantly selfing (t = 0−20%). Although t is commonly computed using progeny arrays (e.g., Sun and Ganders, 1986, 1988, 1990; Sun and Ritland, 1998; Sun et al., 1998), it also can be indirectly estimated based on Wright's within-population inbreeding coefficient (Wright, 1965, 1969) using the relationship t = (1 − F_IS)/(1 + F_IS) as described by Weir (1996). Our present estimates of t based on F_IS were similarly high in R. apiculata and R. mucronata, and the estimates were still higher for populations of R. stylosa. Accordingly, the breeding systems of these species can be defined as mixed mating or predominantly outcrossing. Given that the population F_IS value might represent an overestimate of inbreeding if null alleles are present at some of the marker loci, the estimate of t could even be higher. Hence, a shift in mating system toward inbreeding is unlikely a valid explanation for the generally low levels of genetic variation within their populations.

Several previous studies reported various levels of inbreeding in populations of Rhizophora from both IWP and AEP regions (e.g., Lowenfeld and Klekowski, 1992; Dominguez et al., 1998; Cerón-Souza et al., 2012; Yahya et al., 2013), as well as in populations of other mangrove species, such as Avicennia marina (Maguire et al., 2000), Ceriops decandra (Tan et al., 2005), and Kandelia candel (Chen et al., 1996; Sun et al., 1998). For example, field observations of bee activities suggested that there might be substantial geitonogamous selfing in populations of K. candel and thus quantitative analysis of the mating system parameters was performed using progeny arrays assayed for allozyme markers (Sun et al., 1998). Although the breeding system of K. candel can be regarded as mixed-mating based on the estimated outcrossing rate in two of its populations (0.697 ± 0.091 and 0.797 ± 0.062 respectively), the levels of allozyme variation were extremely low at both population and species levels in comparison with other plant species with mixed-mating system. Random genetic drift and founder effects are most likely the driving force in mangrove populations, especially in the species' distributional margins, overriding the effect of breeding system on the level and pattern of genetic diversity.

With rare exceptions (e.g., Lowenfeld and Klekowski, 1992; Sun et al., 1998), the estimates of F_IS have been commonly used to infer inbreeding in many mangrove studies. This is probably due to the difficulties to assay progeny arrays with mangrove propagules. However, microsatellite-based estimates of F_IS should be interpreted with caution if null alleles are present at one or more of the sampled loci, because null alleles could result in an excess of observed homozygosity compared to HWE expectation and consequently result in an overestimate of F_IS at the locus. Null alleles are commonly present at microsatellite loci (Chapuis and Estoup, 2007), thus the estimates of F_IS at each locus should be carefully checked before inferring inbreeding at the population level.

In this study, we also detected significant F_IS in four of the 20 populations of R. apiculata, four of the 12 populations of R. mucronata, and five of the 22 populations of R. stylosa. After carefully examining the estimates of F_IS at each locus in these 13 populations, we found that the deficiencies of observed heterozygosity in all the populations with significant F_IS could be explained by the presence of null alleles at some of the marker loci in these populations, because other loci in the same individuals had nearly zero or even negative F_IS. Because a species' mating system affects all genetic loci equally, whereas null alleles affect only those loci where they are present, thus near zero or negative F_IS at some of the individual loci provides stronger evidence for random mating or outcrossing than any high values of positive F_IS at different loci in the same population as evidence for inbreeding. We think heterogeneity in F_IS across the sampled loci within the same population, rather than the mean value of F_IS over these loci, should be taken into account before inferring inbreeding at the whole population level.

We found that F_IS values were highly heterogeneous among the marker loci in the 13 populations with significant F_IS, suggesting that the significant F_IS values were likely caused by the presence of null alleles at some of the loci in these populations, rather than reflecting the actual population level of inbreeding. To prove this might be the case, we further analyzed the relationship between null allele frequencies and F_IS estimates at each of the sampled loci for all populations included in the present study, and also did the same analysis for the data given in Wee et al. (2015; original data from http://datadryad.org/resource/doi:10.5061/dryad.42711/1). We found highly significant correlations (p < 2.2e⁻¹⁶) between the two estimates in both of the studies (Figure 5). This relationship indicates caution should be exercised in inferring inbreeding from microsatellite-based F_IS estimates in all future population genetic studies. Because it is often impractical to eliminate all the loci with null alleles if they are present in only a few sampled populations but not in most others, a simple solution to this problem is to inspect heterogeneity in the F_IS estimates across all sampled loci for each population, as we have done here in data analysis and interpretation.

We detected a generally low level of gene flow between pairs of populations within each Rhizophora species, as most of these populations were separated by a vast body of sea water or land barriers. In contrast, gene flow can be relatively high between neighboring populations along the same coastline. We also detected a strong directional asymmetry in the level of gene flow between some of these adjacent populations. For example, the direction of gene flow was predominantly northward between AUS.TRI and AUS.DTR, and westward between AUS.EMB and AUS.CTR; and this same pattern of directionality was detected in all the three Rhizophora species. Therefore, to fully understand the levels of gene flow and pattern of population genetic differentiation, the simple “isolation by distance” model alone may not be applicable for mangrove species. Other factors, such as directions of wind and local current, costal topographies as well as land barriers need to be considered in addition to geographic distances between populations.

High levels of gene flow and hence low genetic differentiation between populations are often associated with long-lived trees with predominant wind pollination and water-dispersed seeds. But this is not the case in all the three IWP Rhizophora species. Despite that their large water-buoyant propagules could preserve viability after floating for long periods of time in salt water (Davis, 1940; Rabinowitz, 1978; Drexler, 2001), the presumed capacity for long-distance dispersal via ocean currents was not supported in the present study, as indicated by the generally low level of gene flow and high genetic differentiation among their populations. Wee et al. (2014, 2015) also found similarly high levels of population genetic differentiation in R. mucronata and R. stylosa. Other physical factors, such as estuarine geomorphology and the interplay between water and wind currents likely played important roles in facilitating or blocking gene flow. For instance, Van der Stocken et al. (2015) showed that whether mangrove propagules could reach the open ocean depends on the strength and direction of both outgoing water flow and winds. Other studies have found that the dispersal of Rhizophora seedlings was mostly restricted to narrow geographic scales (Cerón-Souza et al., 2012; Islam et al., 2013). High levels of genetic differentiation were also found in other water-dispersed mangrove species, including Avicennia (Maguire et al., 2000; Dodd and Rafii, 2002), Ceriops tagal (Ge and Sun, 2001; Huang et al., 2012), and B. gymnorhiza (Minobe et al., 2009). Land barriers, vast oceanic surface distance, and directionality and dynamic interactions of wind and current, whether acting alone or in synergy, apparently played important role in determining the levels of gene flow between populations in most of the mangrove species.

Although the presence of null alleles may lead to an overestimation of genetic differentiation (Chapuis and Estoup, 2007), the effects could be negligible if they were uncommon or rare, i.e., at a frequency of < 0.2 (Carlsson, 2008; Queloz et al., 2010). The MICRO-CHECKER results showed that null alleles were not present at any of the 13 microsatellite loci in 16 populations, and the estimated null allele frequencies in the remaining 38 populations were generally low, only 45 out of 702 locus-population combinations in excess of 0.2. Because none of the null alleles at any of the microsatellite loci was commonly present in these populations, all the 13 loci were used in our data analyses. Still the possible effects of null alleles on F_ST estimates should be taken into account. Thus, we compared the F_ST estimates from FSTAT (null alleles included) with those from FREENA (null alleles excluded). We found that the estimates of genetic differentiation were comparably high whether or not the null alleles were included in the data analyses (Supplementary Table 2).

It should be noted, however, the estimates of genetic differentiation between populations are also found to be closely correlated with maximum geographic distance between the analyzed populations in outcrossing species (Nybom and Bartish, 2000). This indicates the importance of sampling strategies in population genetic studies. Limited coverage of the species' geographic distribution in population sampling could potentially bias F_ST estimates downward. Thus, geographic range of the samples should closely match the species' distributional range in the analysis of population genetic structure, as we have done in the present study.

At individual species level, genetic divergence among populations may not follow the same geographic pattern. In this study, the STRUCTURE analysis showed that populations of R. apiculata were highly divergent, forming 7 small clusters among the 20 sampled populations. However, only two main clusters were present in the neighbor-joining analysis of genetic distances among these populations, showing a deep division between Australasian and Indo-Malaysian populations, coincident with the marine Wallace's Line (Wallace, 1863). This pattern of two notable clusters of populations of R. apiculata is concordant with its phylogeographic disjunction previously revealed based on morphological traits (Duke et al., 2002) and cpDNA and ITS sequences (Lo et al., 2014), suggesting the effect of ancient separation of Australia and SE Asia geological plates. Subsequent genetic differentiation between the Pacific island populations and those of Australia could be caused by the rarity of long-distance dispersal and hence lack of subsequent gene flow following the initial colonization of the islands from source populations in Australia. Within SE Asia, populations on the western side of Malay Peninsula were separable from those on the eastern side, suggesting that the Malay Peninsula acted as a land barrier to the dispersal of R. apiculata, which is consistent with previous studies on Rhizophora species (Inomata et al., 2009; Ng et al., 2015; Wee et al., 2015). This pattern of population differentiation is also evident in other sea water dispersed mangrove species, such as Ceriops species (Ge and Sun, 2001; Liao et al., 2007; Huang et al., 2008, 2012), B. gymnorhiza (Minobe et al., 2009; Urashi et al., 2013), and Lumnitzera species (Su et al., 2006, 2007; Li et al., 2016). Similarly, genetic structure of R. mangle, R. racemosa, and R. samoensis in the AEP region is influenced by the Central American Isthmus and the American continent land masses (Núñez-Farfán et al., 2002; Takayama et al., 2013).

Genetic divergence among populations of R. mucronata and R. stylosa did not exhibit the same pattern as seen in R. apiculata. All the Australian populations of both R. mucronata and R. stylosa were still well separated from those of SE Asia, but their Pacific island populations were more similar to SE Asian than to Australian populations, as reflected by DAPC, STRUCTURE, and NJ analyses. This pattern of population genetic divergence could not be explained by ancient vicariant events, which were expected to affect co-distributed taxa simultaneously and thus produce similar patterns of phylogeographic disjunction, as revealed previously in both R. apiculata and R. mucronata/stylosa in the IWP region with DNA sequence data (Lo et al., 2014). The incongruity between population genetic divergence and phylogeographic disjunction might arise from the difference in the rate of evolution at different molecular marker loci. Microsatellites evolve much faster than DNA sequences, especially faster than cpDNA sequences in plant species (Freeland et al., 2011). Rosenbaum and Deinard (1998) have pointed out some problems with attempting to resolve evolutionary relationships between taxa or among intraspecific populations through microsatellite analysis. Nevertheless, microsatellites continue to be the marker of choice in a wide variety of studies, especially for investigation of population genetic structure and paternity analysis (Merritt et al., 2015).

Our STRUCTURE and gene flow analyses of combined populations (including two or all three species together) showed that R. apiculata was indeed distinct from R. mucronata and R. stylosa, but the species boundary between the latter two taxa was blurry (Supplementary Figures 3, 4). For example, in the combined STRUCTURE analysis including both R. mucronata and R. stylosa, the Australian populations AUS.TRI and AUS.DTR of R. mucronata were grouped together with those of R. stylosa from the same sites (Supplementary Figure 3), similar to the finding based on NJ analysis (Figure 4). These findings are congruent with the extent of natural hybridization as well as phylogenetic relationships among the three species. Natural hybridization between R. apiculata and R. stylosa or R. mucronata resulted in sterile F1 hybrids, named as R. x lamarckii or R. annamalayana (e.g., Kathiresan, 1995; Chan, 1996; Ng et al., 2013). These hybrid sterilities indicate a strong postzygotic isolation mechanism maintaining morphological and genetic distinctions between sympatric R. apiculata and R. mucronata/R. stylosa. In contrast, R. mucronata and R. stylosa are morphologically similar and difficult to distinguish even by experts in the field (Duke et al., 2002). When all 54 populations of the three species were combined in the same analysis, the level of gene flow between R. mucronata and R. stylosa was many times higher than between R. mucronata/R. stylosa and R. apiculata (Supplementary Figure 4). However, no putative hybrids between R. mucronata and R. stylosa have been described based on morphological traits, as the differences between the parental species are subtle and the identification of morphologically intermediate individuals is nearly impossible. Ng and Szmidt (2015) provided genetic evidence for substantial introgression between R. mucronata and R. stylosa at all studied locations where they co-occur, indicating the absence of reproductive barrier. Without pre- or postzygotic barriers, interbreeding between R. mucronata and R. stylosa would be equivalent to gene flow within the same species. Any existing morphological or genetic divergence between the two parental taxa, in the absence of strong local selection for their maintenance, would tend to disappear at the sites of sympatry. This could explain the difficulties in separating R. mucronata from R. stylosa in the two Australian populations AUS.TRI and AUS.DTR where they coexist.

Despite the difficulties in distinguishing R. mucronata and R. stylosa based on morphological keys (such as leaf size and length of flower styles; Duke et al., 2002), there seems to be differential habitat preference: R. stylosa frequently occurs in hard sandy soil and more exposed offshore sites, whereas R. mucronata occurs mostly in areas subject to regular freshwater flows, at least in the eastern part of its range (Duke, 2006). These observations suggest that environmental conditions might have exerted heavy influence on the morphological characters used to delimit R. stylosa and R. mucronata. Also, the morphological or genetic divergence between the two taxa likely developed independently in multiple localities, resulting in polyphyletic origins of both R. mucronata and R. stylosa (Lo et al., 2014). Thus, it seems appropriate to consider R. mucronata and R. stylosa as different ecotypes of the same species, or at most as incipient ecological species, assuming that the speciation process is not reversed due to ongoing hybridization/gene flow. Genome wide SNP markers can be explored to shed further light on the pattern of genetic divergence, hybridization, and introgression between populations of R. mucronata and R. stylosa.

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.