Population genetic patterns of a mangrove‐associated frog reveal its colonization history and habitat connectivity

To understand the impact of historical and contemporary habitat distributions and connectivity on the spatial patterns of a species' genetic variation and divergence, we examined the phylogenetic relationship, population genetic structure and demographic history of a mangrove‐specialist, the crab‐eating frog (Fejervarya cancrivora) from 10 geographic populations in China and northern Vietnam.


| INTRODUC TI ON
Contemporary spatial patterns of a species' genetic diversity are often shaped by a combination of historical and current habitat availability and connectivity. Historical processes such as climate changes, glaciations and palaeogeographic formations may be major drivers of the extant patterns of population genetic diversity (Julio Pineros & Gutierrez-Rodriguez, 2017;Nali et al., 2020;Pinceel et al., 2013;Strugnell et al., 2012). Contemporary factors, including environmental conditions and anthropogenic disturbances, can also strongly impact the geographic distribution of population genetic variation (Cole et al., 2016;Coleman et al., 2018;Wang et al., 2017).
Determining how past and present factors contributed to observed genetic structures enables a deeper understanding of species' evolutionary histories and ecological adaptations. It also contributes to the formulation of appropriate conservation and management strategies (McCartney-Melstad & Shaffer, 2015).
Mangrove forests occur in intertidal zones of tropical and subtropical coastlines and provide essential ecosystem services as well as critical habitats for land and marine fauna (Gopal & Chauhan, 2006;Kristensen et al., 2008;Nagelkerken et al., 2008). However, because they are particularly vulnerable to environmental disturbances, sealevel fluctuations and climate change, mangrove forests are rapidly disappearing and now rank as one of the world's most threatened ecosystems (Alongi, 2002;Duke et al., 2007). Additionally, many mangrove-associated species are evolutionarily unique and of important conservation value (Daru et al., 2013). Conversely, understanding the evolutionary history and population genetic patterns of mangrove-dependent species will help reveal the past distributions and recent conditions of mangrove forests and inform their conservation management.
The crab-eating frog, Fejervarya cancrivora (Gravenhorst, 1829), known for its exceptional adaptation to high salinity habitats throughout all its life stages (Dunson, 1977;Gordon et al., 1961;Uchiyama et al., 1990), is commonly found in tidal mangrove swamps and beaches (Gordon et al., 1961) and is widely distributed in Southeast Asia, including the coastal regions of southern China, Vietnam, Thailand, Bangladesh, Malaysia, Indonesia and the Philippines (Kurniawan et al., 2010). Highly sensitive to habitat disturbance, F. cancrivora has been suggested as an indicator species for mangrove habitat conditions (Hong et al., 2011). This frog's unique adaptations make it an ideal species to investigate the dynamics of mangrove forests. However, the taxonomic relationships and evolutionary history of F. cancrivora's various geographic populations remain a subject of debate. For instance, previous studies have suggested the inclusion of several cryptic species within this lineage. Based on allozyme and mitochondrial DNA (mtDNA) evidence, Kurniawan et al. (2010) identified three genetically divergent groups within the currently recognized F. cancrivora species: the mangrovetype from Bangladesh, Thailand and the Philippines; the large-type from Malaysia and Indonesia; and the Pelabuhan ratu/Sulawesitype (P/S-type) from Indonesia. The mangrove-type shows strict confinement to the mangrove habitat and strong tolerance for salt water, whereas the large-type and the P/S-type do not (Kurniawan et al., 2010). Morphological analyses and crossing experiments confirmed divergences among these groups and supported the proposal that each should be considered a distinct species (Islam et al., 2008;Kurniawan et al., ,2010Kurniawan et al., , , 2011. To date, most phylogeographic studies of the species have focused on Southeast Asian countries, whereas the phylogenetic position and evolutionary history of F. cancrivora populations in China remain obscure. A single specimen from Hainan Island in China was grouped with the mangrove-type based on a phylogenetic tree derived from the mitochondrial 16S rRNA gene, and the same study nested a specimen from Taiwan within the large-type clade (Kurniawan et al., 2010). This discrepancy has raised questions both about the phylogenetic position of China's populations within the species complex and about this species' historical colonization routes in China. Moreover, no systematic population surveys of the geographic distributions and population statuses of this species in China have been conducted. As mangrove forests are disappearing from southern China's coastal regions at alarming rate (Chen et al., 2009;Wang & Wang, 2007), F. cancrivora is likely to have experienced population declines and segregations. Yet, there is a lack of knowledge about the impact of mangrove losses on population genetic patterns of this and other mangrove-dependent species.
In this study, we carried out a nationwide population survey and sampling of F. cancrivora in order to resolve these populations' phylogenetic position within its range and to investigate its population genetic patterns and possible historical and contemporary drivers. To accomplish this, we used a combination of three mtDNA regions and 13 nuclear microsatellite loci of frogs collected from the coasts of southern China and northern Vietnam to reconstruct phylogenetic relationships, assess population genetic structure and infer population demographics and divergence history. In addition to revealing the historical coastal configurations and present mangrove habitat populations. This study provides a framework for understanding the population history of mangrove-associated biota and thus aids conservation management of this critical ecosystem.

K E Y W O R D S
approximate Bayesian computation, gene flow, habitat fragmentation, mangrove biodiversity, phylogeography, population divergence, Southeast Asia coasts conditions in China, our results aid conservation planning that can preserve and restore mangroves and their associated biodiversity.

| Field survey and sampling
We searched for F. cancrivora along the coastlines of Guangxi, Guangdong and Hainan Provinces between August and September 2013 to 2016, focusing on areas with moderate to large patches of mangrove forests. Those three provinces include most of the tropical coastal regions in China and 94% of all existing mangrove forests in the country (Chen et al., 2009). Field surveys were carried out in the evening (8:00-11:00 p.m.) and researchers with flashlights followed line transects near mangroves and shorelines for about three hours, collecting putative F. cancrivora samples. We found the species at nine of 14 localities in the three provinces ( Figure 1a) and obtained tissue samples from 79 individuals (Table 1). During our searches, we observed F. cancrivora on mudflats adjacent to or inside relatively large and continuous mangrove forests, but never in inland, freshwater environments or far away from mangroves. Additionally, four samples of F. cancrivora were collected from a coastal site in Nam Dinh Province of northern Vietnam (Table 1). We were unable to estimate population size at each site due to limited transect surveys and we failed to detect the species at five of the surveyed sites. The mangrove habitats at three of those sites (sites 4, 12 and 13; Table 1; Figure 1a) showed small or highly fragmented forest patches, intensive human disturbance and visible pollutions (e.g. oil slick and accumulated shoreline litter). Mangrove forests at the other two sites (Sites 5 and 14) were in good condition (i.e. large forest areas, continuous distribution and limited anthropogenic modification) and likely contained remnant F. cancrivora populations, but stormy weather during the survey period may have thwarted our ability to detect them. We also sampled a number of rice frogs (Fejervarya limnocharis) from freshwater ponds and wetlands near the sampling sites.
Identities of all samples were later verified using molecular species identification (see below). Tissue samples were preserved individually in 95% ethanol and stored at −20°C for later DNA extraction.

| Molecular analysis
Total DNA was extracted from tissue samples using the EasyPure Genomic DNA Kit (TransGen Biotech) then DNA quality and concentration were checked by 1% agarose gel electrophoresis and with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific Inc.).
For each sample, we first determined species identity based on a 490-bp segment of the mtDNA cytochrome b gene (Cytb) sequence (see below). The two morphologically similar species, F. cancrivora and F. limnocharis, diverged by 19.0% (93 bp) in this sequence.
To find phylogenetic relationships and genetic diversity of F. cancrivora, we analysed three mtDNA fragments (Cytb, using primers as described in Table S1, Appendix S1. The designs of F I G U R E 1 (a) Map of field survey locations and geographic distributions of Fejervarya cancrivora mtDNA haplotypes. Filled symbols indicate survey localities where samples were collected and open circles indicate sites where no F. cancrivora were detected. Different symbols show genetic clustering into three clusters, as revealed in population structure analyses based on microsatellite data. Colours indicate different haplotypes (ChMT-1-9) of concatenated mtDNA sequences (total length = 1,818 bp) and pie charts show the haplotype proportions found in each sampling locality in this study. (b) Median-joining haplotype network of ChMT-1-9 in (a). Each unique haplotype is represented by a circle. Circle size indicates the number of individuals (varying between 1 and 46) and each transverse dash on a line indicates a substitution between haplotypes our primers for Cytb and the D-loop region were based on F. cancrivora mitochondrial genome sequences (GenBank accession no. EU652694; Ren et al., 2009). We conducted PCR in total volumes of 25 μl, containing approximately 2 ng tissue DNA, 10 μM of each primer, and 12.5 μl 2×EasyTaq PCR SuperMix (TransGen Biotech).
The PCR conditions included an initial denaturing at 95°C for 5 min; followed by 35 cycles at 95°C for 30 s, the annealing temperature for each primer pair (see Table S1, Appendix S1) for 30 s and 72°C for 1 min; and a final extension step at 72°C for 10 min. We purified the PCR products using the EasyPure Genomic DNA Kit For all samples, we amplified a total of 13 microsatellite loci (A33, A127, B5, B19, B34, B43, B46, B47, B72, B147, B152, C46 and C63) using fluorescent-labelled forward primers and PCR conditions following Zheng et al. (2016). These loci were selected from a panel of variable microsatellite markers developed specifically for F. cancrivora in China. We analysed the PCR products on a 3730 XL DNA Analyzer with a GeneScan 500 LIZ size standard (Applied Biosystems) and genotyped them using GENEMAPPER version 3.7 (Applied Biosystems).

| Genetic diversity and phylogenetic relationships inferred from mtDNA
We used CLUSTAL W (Larkin et al., 2007) to align mtDNA sequences.

MK396081-MK396095
We estimated the phylogenetic relationships between different F. cancrivora clades using the maximum likelihood (ML) and Bayesian inference (BI) methods based on the Cytb and 16S haplotypes from both China and other Southwest Asian countries (see Table S2, Appendix S1). Corresponding F. cancrivora D-loop sequences from other countries were unavailable for this analysis. The corresponding sequences from Fejervarya iskandari (GenBank acc. nos. AB296085 and AB277303) and F. limnocharis (GenBank acc. nos. MK411008 and AB070737) were used as outgroups. Using  In addition, we estimated the divergence times between different F. cancrivora clades with BEAST version 2.6.3 (Bouckaert et al., 2019) using the Cytb + 16S haplotypes. Due to the lack of fossil-based calibration points, we employed a lognormal relaxed clock (uncorrelated), an average mtDNA substitution rate in anurans of 1% per million years (Yan et al., 2013;Zhang et al., 2010) and randomly generated starting trees. We ran the MCMC chain for 40 million generations and sampled every 4,000 generations. A maximum credibility consensus tree was constructed using TREEANNOTATOR with the initial 10% of sampled trees discarded. The result was checked with TRACER for chain convergence and posterior sample size. We required an effective sampling size of more than 200 for all parameters.
To better visualize phylogenetic relationships among F. cancrivora populations sampled in this study, we constructed a haplotype network for each of the three mtDNA regions and for the concatenated sequences (a total length of 1,818 bp) of each sample using the median-joining method in NETWORK version 5.0 (Bandelt et al., 1999).  (Raymond & Rousset, 1995), and for multiple comparisons, p values were adjusted using the Bonferroni correction (Rice, 1989). We quantified the genetic diversity of the microsatellite loci in each sampled population by using the following estimates:
Isolation by distance (IBD) refers to a positive correlation of genetic distance and geographic distance that is often observed in species with spatially limited gene flow (Rousset, 1997;Wright, 1943).
We used Mantel tests implemented in GENALEX version 6.5 (Peakall & Smouse, 2006 to assess the correlation between pairwise Nei's genetic distance (calculated with GENALEX) and geographic distance among sampled individuals. Geographic distance among samples was measured as the natural logarithm of the present coastline distance. As F. cancrivora is closely associated with coastal mangrove habitats, coastline distance is a more biologically realistic measure of dispersal than is straight-line distance (Karns et al., 2000). Present coastline distances between our sampling locations were obtained using Google Earth (http://www.

google.com/earth/) and the open ocean between Hainan Island and
the Leizhou Peninsula was measured as the shortest distance across the Qiongzhou Strait. Statistical significance was evaluated by comparing the estimated r to a distribution of r-scores calculated with 10,000 permutations.
We used two methods to infer population genetic structure based on our samples' microsatellite data. First, we used a Bayesian clustering algorithm in STRUCTURE version 2.3.4 (Pritchard et al., 2000) in which the value of the number of genetic clusters, K, was set to 1-10, with the largest K equal to the number of sampled localities. Ten independent runs were performed for each K value, using the admixture model with correlated allelic frequencies among clusters. Simulations were run for one million MCMC iterations after a burn-in of 100,000. The most likely K value was determined using the ΔK method following Evanno et al. (2005), as implemented in STRUCTURE HARVESTER version 0.6.94 (Earl & vonHoldt, 2012).
To detect the potential hierarchical population structure nested within the first level of genetic clustering, we repeated an iterated STRUCTURE analysis on each identified cluster until the genetic cluster was inferred down to one (Evanno et al., 2005;Janes et al., 2017).
The STRUCTURE results for the most likely value of K were summarized using CLUMPP version 1.1 (Jakobsson & Rosenberg, 2007) with 1,000 permutations and the FullSearch algorithm and plotted using DISTRUCT version 1.1 (Rosenberg, 2004).
We also performed principal coordinate analysis (PCoA) using GENALEX to assess population genetic structure based on genotypic genetic distances between individuals.
Estimates of the effective number of migrants (N m ) among the inferred clusters were assessed using a Bayesian coalescent-based framework under the default setting in MIGRATE-N version 3.6.11 (Beerli, 2006;Beerli & Palczewski, 2010). MIGRATE-N provides estimates of θ (4N e μ; N e = effective population size, μ = mutation rate) and M (m/μ; m = migration rate) from microsatellite data. We ran five independent replications for each genetic cluster under the Brownian stepwise mutation model with initial parameters based on F ST values and constant mutation rates. We conducted searches using 5 million steps for each replication, with 500,000 discarded as burn-in, and sampled the parameter values every 20 iterations. We calculated N m using the equation N m = (θ × M)/4.

F I G U R E 2
Six likely scenarios of population divergence in the DIYABC analysis using Fejervarya cancrivora microsatellite data. Each scenario represents a divergence order among the three genetic clusters inferred from microsatellite data. Time is not strictly to scale. See Figure 1 for map labels

| Demographic history and population divergence
The different mutation rates of mtDNA and microsatellite loci allow inferences of population demographic patterns in different time frames. Slowly evolving mtDNA can reveal ancient population demographic history, whereas the rapid mutation rates of microsatellites enable tests of ancient as well as recent demographic events (Fumey et al., 2018;Zhang et al., 2014). We evaluated the ancient historical demographic changes of the inferred genetic clusters using the combined mtDNA sequence (1,818 bp) data. Because genetic Clusters 2 and 3 each contained few mtDNA haplotypes, we combined the data of those two clusters for the subsequent analyses.
First, we estimated Tajima's D (Tajima, 1989) and Fu's F S statistics (Fu, 1997) in DNASP to seek the genetic signature of demographic expansions. We then performed Bayesian skyline plots with BEAST to assess ancient demographic dynamics in the study populations.
This analysis was conducted using the same settings as above for divergence time estimation. The data failed to generate an effective sample size (n > 200) when Clusters 2 and 3 were analysed both separately and combined, likely because of the limited sample sizes of those populations. Therefore, we conducted this analysis only on data from Cluster 1 and the total population.
We used our microsatellite data to analyse population divergence history through approximate Bayesian computation (ABC) implemented in DIYABC version 2.0.1 (Cornuet et al., 2008(Cornuet et al., , 2014. diverged from Cluster 1 at t 2 and Cluster 3 split from Cluster 2 at t 1 (i.e. a dispersal route from Hainan to Guangdong, then to Guangxi).
In Scenario 5, Cluster 3 diverged first from Cluster 1 at t 2 and then Cluster 2 split from Cluster 3 at t 1 (i.e. a dispersal route from Hainan to Guangxi, then to Guangdong). Finally, in Scenario 6, a population diverged from Cluster 1 at t 2 and subsequently split to Clusters 2 and 3 at t 1 (i.e. a population diverged from the Hainan population and split to the Guangdong and Guangxi populations) (Figure 2a). We set prior values for N e and divergence time estimates with a uniform distribution for all parameters (Table S3, Appendix S1). In all scenarios, we assumed that each of the three clusters had a constant size over time and post-divergence from the ancestral population. The microsatellite mutation rate was set between 1 × 10 −6 and 1 × 10 −4 substitutions/generation (Komaki et al., 2017) with a uniform distribution and under the stepwise mutation model. All loci were set to 4-bp repeats except A33 and A127, which were set to 2-bp repeats (Zheng et al., 2016). We used default settings for the other microsatellite parameters. We simulated six million datasets for each scenario and calculated summary statistics (mean number of alleles per locus, mean genetic diversity and mean Garza-Williamson's M) for each cluster and F ST and the mean classification index between pairs of clusters. We estimated the posterior probabilities of the modelled scenarios using a logistic regression plot of the 1% of simulated datasets closest to the observed data (Cornuet et al., 2008). The most likely scenario was the one with the highest significant posterior probability value and non-overlapping 95% confidence interval (CI).
For this scenario, we estimated type I (probability of rejecting the selected scenario when it is true) and type II (probability of selecting the scenario when it is false) errors by using the "confidence in scenario choice" function in DIYABC with 1,000 independent datasets and logistic regression approaches (Cornuet et al., 2010). We also used principal component analysis in DIYABC to check goodness-offit between the simulated and observed data by using 10,000 simulations of the most likely scenario.
Finally, we used the microsatellite data with BOTTLENECK version 1.2.02 (Cornuet & Luikart, 1996;Piry et al., 1999) to detect the molecular signatures of recent population declines. We conducted the analysis under both the stepwise (SMM) and the two-phased mutation models (TPM), with 95% single-step mutations and 5% multi-step mutations. We tested for the significance of bottleneck signatures using the sign test and Wilcoxon's signed-rank test.

| Genetic diversity and phylogenetic relationships inferred from mtDNA
Three mtDNA fragments (Cytb, 16S and the D-loop region) were successfully sequenced for all 83 samples. Six, five and four unique haplotypes were recovered from the Cytb, 16S and D-loop sequences, respectively, and combining the three fragments from each individual yielded nine haplotypes (ChMT-1-ChMT-9). The number of haplotypes and diversity estimates for each sampling locality are shown in Table 2.
We constructed phylogenetic trees using Cytb and 16S sequences from our samples and from F. cancrivora samples from other Southeast Asian countries (Figure 3). Both ML and BI trees showed similar topologies ( Figure S1, Appendix S1). All haplotypes from our study grouped with the mangrove-type haplotypes from the Philippines, Thailand, India and Bangladesh. We estimated the divergence event between the mangrove-type and the clade of the large-type from Malaysia and Indonesia and the P/S-type from Indonesia to be 9.95 Mya (95% highest posterior density [HPD]: 7.40-13.70) and the large-type and the P/S-type clades split 5.99 Mya (95% HPD: 4.23-8.44). Sequence divergences estimated by p-distances between our samples and mangrove-type haplotypes from other Southeast Asian countries ranged from 0.2% to 3.9% (mean 1.8% ± SD 1.4%) for Cytb and 0.0% to 1.5% (mean 0.5% ± SD 0.5%) for 16S haplotypes. In comparison, Cytb haplotypes from our study showed 13.7%-15.7% (mean 14.5% ± SD 0.5%) sequence divergences from the large-type haplotypes and 15.3%-16.5% (mean 15.8% ± SD 0.4%) from the P/S-type haplotypes, respectively. While our samples' 16S haplotypes showed 8.7%-9.4% (mean 9.0% ± SD 0.2%) sequence divergences from the large types and 9.9%-10.5% (mean 10.2% ± SD 0.2%) from the P/S types, respectively.

| Genetic diversity and population structure inferred from microsatellite data
Samples from 83 frogs collected at 10 localities were successfully PCR-amplified and genotyped. Seventy-three (88%) samples were genotyped at 12-13 microsatellite loci, while 10 lacked genotype data at two to three loci. The number of alleles varied from 5 to 27, with a mean of 13 (SD 6.3) alleles per locus (Table S4, Appendix S1).
Genetic diversity estimates for each sampled population are shown in Table 2. Three, two and zero loci showed significant departure (p < .05) from Hardy-Weinberg equilibrium in the inferred genetic Clusters 1, 2 and 3 (see below), respectively, while no linkage disequilibrium was found between any loci. Two loci showed evidence of null alleles in Cluster 1, but estimated percentages of null alleles were low (< 8%). A greater number of loci showed null alleles in Clusters 2 and 3, but this pattern was likely due to small sample sizes so we retained those loci in the following analyses.

TA B L E 2
Summary of genetic diversity parameters of the mtDNA (combined sequences of 1,818 bp) and 13 microsatellite loci of Fejervarya cancrivora sampled from 10 localities in China and northern Vietnam The Mantel test of microsatellite data revealed a significant positive correlation between pairwise genetic distance and geographic distance (coastline distance) among sampled individuals (Figure 4) MIGRATE analysis showed migration rates (N m ) ranging from 0.1 (from Cluster 3 to 1) to 0.9 (from Cluster 3 to 2) individuals per generation (Table 3), thus indicating low levels of gene flow between all cluster pairs.

| Demographic history and population divergence
We inferred the ancient population demographic history using the  Table S5, Appendix S1). This suggests a historical lack of large population expansion or contraction.
The Bayesian skyline plots showed overall stable historical effective population sizes for Cluster 1 and for the total population ( Figure S4, Appendix S1).
The most supported scenario of F. cancrivora population divergence history, using DIYABC analysis, was Scenario 6 ( Figure 2 and Figure S5, Appendix S1), in which a first historical divergence event occurred between Cluster 1 and the other populations, and a second event occurred with the divergence of Clusters 2 and 3 from a common ancestral population (logistic approach PP = 0.605, 95% CI = 0.584-0.627). The posterior probability of Scenario 6 was considerably higher than that of Scenario 4 (PP = 0.250, 95% CI = 0.218-0.283), the second highest ranked scenario and that of the other scenarios, which varied between 0.000 and 0.140 (Table S6, Appendix S1). Scenario 6 showed high type I (0.525) but generally low type II errors (mean = 0.138; Table S7, Appendix S1), and a good match between observed and simulated data ( Figure S6, Appendix S1), thus indicating the posterior distributions of this scenario were well corroborated by the observed data. The posterior distribution of N e of the three clusters had reasonable estimates, and the median values of the divergence time, t 2 and t 1 , were 34,600 and 3,250 generations ago, respectively (Table S3, Appendix S1).
BOTTLENECK analysis using microsatellite data to infer recent population bottlenecks detected population decline signatures in Clusters 1 and 2 (under both SMM and TPM) and in the total population (SMM only)and was supported by significant results of both sign tests and Wilcoxon's signed-rank tests (Table 4). Cluster 3 showed non-significant results, but possibly lacked power due to its particularly small size (8 samples). Overall, these results suggest range-wide population declines across Chinese coasts in recent times.

| Phylogenetic relationships and colonization history of Chinese Fejervarya cancrivora
The phylogenetic relationships and taxonomic structure among species of the genus Fejervarya have long been a topic of confusion.
Species identification based on morphological characters is often problematic in this clade, owing to minor interspecific differences and intraspecific morphological variations, and because morphologically close groups may show considerable genetic divergence and reproductive isolation (Islam et al., 2008;Kotaki et al., 2010). Molecular phylogenetic analyses have revealed many cryptic species and erroneous species identifications Sumida et al., 2007;Toda et al., 1998). This all points to a need for careful examination and validation of phylogenetic relationships in the currently recognized genus Fejervarya. Our mtDNA analysis of the F. cancrivora sampled in China and Vietnam provides the first molecular evidence of the phylogenetic position of these populations within the species complex.
Our results placed our samples within the mangrove-type also found in Bangladesh, Thailand and the Philippines, which diverged from the large-and P/S-type clades. The deep divergence history and large mtDNA sequence divergences among the mangrove-, large-and P/S types support that F. cancrivora is a complex of cryptic species and the three types can be considered three distinct species (Kurniawan et al., ,2010(Kurniawan et al., , , 2011. The limited sequence divergences and phylogenetic relationships between our samples and those from the Philippines, Thailand and Bangladesh suggest that the expansion of this species from Southeast Asia to China occurred relatively late. It is known that repeated glaciations during the Pleistocene caused sea levels to drop below the present level for extended periods of time (Voris, 2000;Yao et al., 2009). At low sea levels, the Beibu Gulf area was exposed above sea level and southern Vietnam and the Hainan area were directly connected by a short coastline (Voris, 2000) (Figure 6).
However, the high type I error for the most supported population divergence scenario suggests limited confidence of the ABC analysis, and we cannot definitively prove our dispersal route hypothesis due to the limited geographic scope of our samples. Mangrove-type F. cancrivora samples from other Southeast Asian countries (e.g. India, Bangladesh, Thailand, southern Vietnam and the Philippines) will help resolve the intraspecific relationships and evolutionary history of this species.
An interesting question remains regarding the origin and phylogenetic position of the F. cancrivora population on the island of Taiwan. The 16S rRNA gene sequence of a specimen from southern Taiwan grouped with the large-type haplotypes from Indonesia and Malaysia, but was separate from the mangrove-type (Kurniawan et al., 2010). Samples from the landmass close to Taiwan, including mainland China to its west and the Philippines to its south, were all shown to be the mangrove type (Kurniawan et al., 2010;this study (2000) in southern China. Similarly, an introduced F. cancrivora population, transported from East Asia, was also reported in Guam (Christy et al., 2007).

| Population genetic structure of Fejervarya cancrivora in China
A strong IBD pattern of genetic differentiation was detected among all sampled populations in southern China, thus supporting the likelihood that the species may have a limited capacity for long-distance movement and that most dispersal may occur in small incremental steps between adjacent populations. Similar patterns of IBD were also reported in a coastal mangrove-associated snake (Karns et al., 2000), thus suggesting a commonality among mangrovedependent species in their dispersal modes, which are possibly shaped by the patchily distributed habitat along coastlines.
Genetic analysis of the population structure using the microsatellite data revealed three genetically diverged groups among our sam- N m estimates were still less than one, suggesting the presence of dispersal barriers between the two clusters, despite relatively short coastline distances. Interestingly, the one Vietnam population was genetically grouped with Cluster 2 despite the long geographic distance, thus suggesting frequent inter-population gene flow that may be maintained by either ocean currents in the Beibu Gulf (Daryabor et al., 2016;Gao et al., 2013) or human-mediated transport from Guangxi to Vietnam.

| Conservation implications for mangrovedependent biota
Southeast Asia holds most of the world's mangrove forests (Giri et al., 2011), yet mangroves are rapidly disappearing due to that region's urban and economic developmental pressures (Duke et al., 2007;Lovelock et al., 2015). Mangrove forests provide essential habitats for highly diverse terrestrial and marine biota, and the rapid disappearance of this critical habitat is inevitably threatening its dependent biodiversity (Daru et al., 2013;Polidoro et al., 2010). However, very few studies have investigated the impact of mangrove loss on the population genetic structure of associated species. The strong genetic differentiations among population clusters of F. cancrivora may have roots in historical climatic and geographic processes, but the recent and ongoing areal declines of mangrove forests are likely aggravating the genetic isolation among local populations. Further population surveys and genetic monitoring of this and other mangrove-associated species should be carried out to evaluate the effects of mangrove deforestation and restoration on dependent biota.
Based on the results of our genetic analyses and field survey, we recommend a twofold strategy for mangrove protection and management. First, large forest patches should be vigorously protected because large habitats have the capacity to support rich biodiversity and large populations, generally show higher resilience against environmental disturbances and have better prospects for species persistence in the long term (Ferraz et al., 2003;Ribeiro et al., 2009). This is corroborated by observations that F. cancrivora is most readily found in larger (>100 hm 2 ) mangrove forests (Hong et al., 2011;our observation). However, maintaining large but spatially distant habitat fragments may create genetic isolation among populations, thus hampering the associated species' metapopulation processes and evolutionary potential (Garant et al., 2007;Thomas, 2000). The strong genetic differentiation between population clusters in this study provides an example of such isolation. Therefore, our second suggestion is to preserve and restore a number of smaller (10-100 hm 2 ) and geographically proximate mangrove patches along the coastlines. Smaller patches constitute a considerable fraction of remnant mangroves, often neglected in current conservation policy. Yet, they can play an important role in improving population viability in fragmented landscapes, either by functioning as movement corridors and stepping stones to facilitate gene flow between larger patches (Tewksbury et al., 2002;Uezu et al., 2008) or by forming networks of functionally connected habitat fragments, which would allow species such as F. cancrivora to persist in disturbed landscapes (Martensen et al., 2008;Ribeiro et al., 2009). University to M.Y. and J.Z.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Mitochondrial DNA sequences are available from GenBank (accession nos. MK396081-MK396095). Data used in this study are available from the authors upon request.