Habitat fragmentation is one of the major threats to amphibian species. In a previous study, population genetic analyses of the Japanese brown frog Rana japonica were conducted using a mitochondrial DNA (mtDNA) marker in a typical Japanese agricultural landscape (known as satoyama) in Chiba, Japan. This previous study revealed that gene flow was restricted by the roads and cement-walled urban river that divide this site. In the present study, we reanalyzed the genetic structure of the same meta-population using microsatellite markers in comparison with the mtDNA results and elucidated fine-scale gene flow. The genetic structure derived from the microsatellite clustering analysis was almost identical to that of the mtDNA results, although some important details differed. We recognized boundaries of genetic structure are consistent with the major roads and cement-walled river, however, we also detected gene flow across those artificial barriers. We concluded that the current genetic structure was formed in the past when gene flow was strongly restricted. Gene flow among breeding populations is now being restored by the maintenance of breeding sites, although it is not sufficient to erase the signature of historical isolation.
Introduction
Habitat fragmentation increases the vulnerability of amphibian populations (Stuart et al., 2004; Beebee, 2005; Hamer and Mcdonnell, 2008; Allentoft and O'Brien, 2010). Amphibian species often form metapopulation structures because of their limited breeding habitats (e.g., ponds, rivers, and wetlands) (Marsh and Trenham, 2001; Beebee, 2005; Smith and Green, 2005). Anthropogenic activities such as highway construction and deforestation may divide breeding habitats by forming barriers to migration and gene flow (Arens et al., 2008; Safner et al., 2011; Mikulíček and Pišút, 2012; Hale et al., 2013; van Strien et al., 2014; Kakehashi et al., 2014; Dias et al., 2015). For example, Safner et al. (2011) reported that not only highways but also national roads restrict the migration and gene flow of the European brown frog Rana temporaria, and Arens et al. (2008) reported that urbanization and habitat loss significantly affect the genetic structure and gene flow of the Moor frog R. arvaris.
In a previous study, using a control region of mitochondrial DNA (mtDNA), we found that the genetic structure of the Japanese brown frog Rana japonica was correlated with the locations of major roads (10 m wide) and a cement-walled river. This genetic structure implies that gene flow among breeding sites is impeded by these barriers (Kobayashi et al., 2013). We suggested that agricultural modernization, such as the construction of cement walls in rivers and road improvement, degrades habitat quality for this representative species of a common Japanese agricultural landscape, the satoyama. However, nuclear markers which are more sensitive to fine-scale gene flow were not used in that study.
Nuclear markers such as microsatellite markers can be used to show the fine-scale genetic structure and extensive gene flow of local populations (Rogers and Peacock, 2012; Iwai and Shoda-Kagaya, 2012; Dias et al., 2015). Understanding gene flow and genetic structure using nuclear markers is critical for developing conservation approaches in local populations, because such detailed information provides an improved understanding of the effects of barriers or other strong restrictions on movement among breeding sites and allow for better management (Cushman, 2006; Purrenhage et al., 2009; Safner et al., 2011; Sotiropoulos et al., 2013). Therefore, in this study, we first measured the genetic structure of R. japonica on the western side of Inba Lake, Chiba, Japan using microsatellite markers and compared it with the genetic structure measured using a mitochondrial DNA marker. We also elucidated fine-scale population structure and gene flow across the study site, and assess the implications of these results for effective conservation measures.
Materials and Methods
Study area
Our study was conducted on the western side of Inba Lake, approximately 40 km northeast of Tokyo. We conducted our analyses in a small study area (approximately 3×3 km; Fig. 1i), where we found 13 major breeding sites in March 2010. In the study area, there are major roads (dashed lines, Fig. 1i) and a 3-m-wide urban river with a 2-m-deep concrete embankment built during the 1980s (solid line, Fig. 1i). These structures divide the forest and paddy fields into six parts. We previously reported that these roads and urban river likely restrict gene flow based on mtDNA data (Kobayashi et al., 2013) (Fig. 1ii). The environment surrounding the study site is varied. To the east, there are a number of roads, a lake, and an urban area. The north side of the study area consists of a golf course, a large railway, paddy fields, and an urban area with a large road. To the west of the study area, there is a large golf course. We also conducted census surveys in those surrounding areas (approximately 500-m ranges), but we could not find a major R. japonica breeding site. Therefore, we concluded that the 13 sites in Fig. 1i represent the major breeding sites in this region.
Population sampling
Similar to many other ranid frog species, R. japonica females spawn one egg mass containing 500–3,000 eggs during each breeding attempt. Many individuals congregate to mate in wet paddy fields, ponds, and wetlands during the reproductive season, which usually occurs from February to March (Kaneko and Matsui, 2004). We collected R. japonica eggs from 13 breeding sites within the study area (Fig. 1i) during February and March from 2010 to 2012. The breeding sites were mostly in paddy fields, small drains beside the paddy fields, and wetlands on abandoned paddy fields. We collected 20 egg masses each year at each breeding site (except site 1 in 2011, when we missed a window for egg mass collection, according to the presence of some tadpoles) and approximately three eggs per each egg mass were brought to the laboratory. We reared the eggs to the tail-bud stage at room temperature in the laboratory, at which time one individual per egg mass was euthanized for DNA extraction. The numbers of samples collected and the numbers of egg masses counted in each year are presented in Table 1. A few eggs per breeding site did not reach the tail-bud stage, so we excluded those samples.
Fig. 1.
Geographic distribution of analyzed sample sites and genetic clusters on the western side of Inba Lake. (i) Geographic overview of the study site. The white solid line indicates a cement-walled urban river (approximately 3 m wide), and the dashed lines indicate roads (approximately 10 m wide). The white circles numbered 1–13 indicate the major breeding sites of Rana japonica sampled. (ii) Genetic groups estimated using mitochondrial DNA in a previous study (adapted from Kobayashi et al., 2013). (iii) Genetic structure results from the clustering analysis based on microsatellite markers.
Laboratory protocols
Total DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen Inc., Valencia, CA, USA). We evaluated 10 microsatellite loci developed for R. japonica (Koizumi et al., 2009: Raja01–04, 07, 10–12, 18, and 19) and checked for stuttering and null alleles using Micro-Checker 2.2.3 (Van Oosterhout et al., 2004). Finally, we used six loci for which no bias caused by null alleles was detected. We used a pair of tailed non-labeled oligos and fluorochrome-labeled oligos with different dyes attached to the 5′ end for GeneScan Analysis (Applied Biosystems, Foster City, CA, USA). These markers were amplified in two multiplex and one single-locus polymerase chain reaction (PCR): multiplex-1, Raja01 (NED), Raja03 (6-FAM), Raja04 (VIC); multiplex-2, Raja02 (6-FAM), Raja10 (VIC); and s ingle-locus, R aja19 ( PET). W e c onducted PCR for different microsatellite loci under the same PCR conditions, described as follows: template DNA was added to 10 µl PCR reaction mixture composed of 20 mM Tris—HCl (pH 8.0), 100 mM KCl, 20 mM MgCl2, dNTPs (2.5 mM each), 1 mg bovine serum albumin, 4% dimethyl sulfoxide, primers (20 mM each), and 0.5 U Takara Ex Taq DNA Polymerase Hot-Start Version (Takara Bio, Shiga, Japan). PCR amplification consisted of an initial denaturation step at 98 C for 1 min, followed by 35 cycles of denaturation at 98 C for 10 s, annealing at 60 C for 60 s, and an extension at 72 C for 30 s, followed by a final extension at 60 C for 30 min. Amplified fragments from the microsatellite loci were genotyped using the ABI PRISM 3130xl Genetic Analyzer and GeneMapper 4.0 software (Applied Biosystems). PCR products were mixed with Hi-Di Formamide and 500 LIZ Size Standard (Applied Biosystems) and the mixed solution were used for capillary electrophoresis. R aja19 w as m ixed w ith multiplex-1 and they were electrophoresed at the same time.
Table 1.
Genetic indices for each breeding site and the numbers of samples and the number of egg masses counted each year.
Genetic diversity of each breeding site and genetic structure using clustering analysis
To determine the basic genetic indices for each breeding site, we counted the number of alleles (Na), calculated the observed and expected heterozygosities (Ho and He, respectively), and tested for Hardy—Weinberg equilibrium (HWE) for each locus at each breeding site using GenAlEx 6.5.1 (Peakall and Smouse, 2012). We also corrected the HWE p-values for each breeding site using Bonferroni correction to obtain more robust results. The inbreeding coefficient, which is generally represented by the fixation index of a subdivided population (FIS), was calculated using INEst 2.0 as the average of individual samples using the individual inbreeding model (Chybicki and Burczyk, 2009). Additionally, the allelic richness (AR) calculation and linkage disequilibrium test were performed in FSTAT 2.9.3.2 (Goudet, 2001).
Genetic structure was estimated using two approaches. First, we used Bayesian Analysis of Population Structure (BAPS) ver. 6.0 with the ‘spatial clustering of groups' models followed by admixture analysis (Corander and Marttinen, 2006; Corander et al., 2008a, b; Cheng et al., 2013). BAPS uses both approaches to determine the most appropriate clustering (default mode) and best groups in a fixed number of clusters (fixed-K mode). We used the fixed-K mode, adjusting K from 2 to 10, and ran the analysis 1,000 times for each value of K.
We then used STRUCTURE 2.2.3 to implement an admixture model using the LOCPRIOR option, which uses the sampling location as prior information to assist clustering (Pritchard et al., 2000; Falush et al., 2003; Hubisz et al., 2009). The parameters were as follows: burn-in length, 2,000,000; running length, 2,000,000; and K, 2 to 10, with 20 runs for each value of K. The results were summarized using Structure Harvester 0.6.94 (Earl and von Holdt, 2011) using a table of summary statistics and graphs of each parameter such as log-likelihood. We determined genetic groups based on the consensus of these two clustering analyses.
Gene flow
To elucidate fine-scale gene flow more analytically, we estimated gene flow based on recent migration rates (m) among the genetic groups using an assignment test in BayesASS 3.0.3 (Wilson and Rannala, 2003). BayesASS estimates the percentage of individuals in a population that are assigned to each population in the study area. It should be noted that BayesASS does not work well when the tested population includes multiple gene pools (Wilson and Rannala, 2003). Therefore, we decided to estimate m among the estimated genetic groups. We performed four separate runs using two tuning parameters, delta M and delta F, to change the acceptance-rate fitting following the BayesASS manual instructions. We tested values for delta M, a mixing parameter for m, of 0.10 and 0.20. We tested values for delta F, a mixing parameter for inbreeding coefficients, of 0.40 and 0.50. We conducted this twice for each parameter combination with different seeds (10 and 3099 respectively), thus obtaining results for eight runs. The settings of the other parameters were the same in each run. We obtained the average of 10,000,000 Markov chain Monte Carlo iterations, sampling the chain every 100 iterations; the first 1,000,000 iterations were discarded as burn-in. On the other hand, the sample sizes were sufficiently different among the genetic groups, thus requiring adjustment. Consequently, we estimated the number of migrating individuals per genetic group by multiplying m with the census size (Nc). To determine Nc for each genetic group, we first doubled the egg mass number counted each year (we assumed this was the number of breeding individuals for that year) and then calculated the average number of breeding individuals over the 3 years. We averaged the m values obtained from eight runs in BayesASS and used the product of m and Nc to estimate the number of migrating individuals. Then, we evaluated the present state of gene flow among these 13 breeding sites and discussed effective conservation measures.
Results
Genetic diversity at each breeding site
The calculated basic genetic indices for each breeding site are shown in Table 1; there was no noticeable decrease in AR among sites. The FIS values for each site ranged from 0.009 to 0.081. Some loci deviated from HWE at various breeding sites (Raja01 at site 9, Raja03 at site 2, Raja04 at site 12, and Raja19 at sites 1, 4, 9, and 12), but no loci deviated from HWE at more than four sites. There was also no significant linkage disequilibrium among loci.
Genetic structure according to clustering analysis
The BAPS results indicated an optimal K value of 5 (Fig. 2). The highest mean log likelihood of K over all reps for that K (LnP[K]) in the STRUCTURE analyses was found at K=5 with the admixture model (Fig. 3). The genetic structures in the BAPS and STRUCTURE analyses were nearly identical (Figs. 2 and 3), and the likelihood had a single peak at K=5 (Fig. 3). However, there were some differences in the clustering results between the two analyses. In the STRUCTURE analysis, all individuals from sites 6 and 7 did not consist of one cluster, but rather they included three clusters: cluster I (also at sites 1–5), cluster II (also at site 8), and cluster III (also at sites 11–13). Cluster II was found only at sites 6, 7, and 8, whereas in BAPS, most individuals from sites 6, 7, and 8 were clearly grouped into one cluster (cluster II). Therefore, we defined five genetic groups (A–E; Fig. 1iii) based on the result of BAPS in the next gene flow analysis step. Moreover, the microsatellite clustering results differed from the mtDNA genetic structure at breeding sites 9–13 (Fig. 1).
Gene flow
We estimated m using BayesASS as an indicator of gene flow among the genetic groups. The BayesASS results using different parameters were shown in Table 2. The results of duplicated calculations using the same parameter sets with different seeds were the almost same, as summarized in Table 2. The m values using different parameter sets were almost identical except for those from group C to C and from E to C, which ranged from 0.778 to 0.835 and 0.117 to 0.175, respectively (Table 2). Migration rate (m) were relatively higher from group A to B, E to B, E to C and E to D than other directions. Taking 95% confidential intervals into account, the m-values for the within site migration rate never overlapped with zero. Conversely between site values for only two site pairs did not overlap with 0 (specifically, from A to B, and from E to D). The values of the product of m and Nc, which represents the number of migrating individuals, are shown in Table 3. The estimated numbers of migrating individuals and the value of m were high in the original genetic groups. Excluding the original genetic group, both m and the estimated number of migrating individuals from group A to B, E to B, and E to D were relatively high.
Discussion
Genetic diversity and genetic structure
In a previous study at the same sites Kobayashi et al. (2013) detected much lower genetic diversity indices at sites 6, 7, and 9 when compared to other sites, and based on mtDNA. Only two mtDNA haplotypes were found from site 9 and only four haplotypes were detected at sites 6 and 7, respectively. These indices were lower than the average (7.7) among the other 10 sites (Kobayashi et al., 2013). The authors interpreted lower mtDNA haplotype diversities at these sites as a result of lost genetic diversity caused by habitat fragmentation (Kobayashi et al., 2013). However, in this study, the genetic diversity calculated based on microsatellite markers did not show such obvious lower values at these sites. There are a few possible explanations for this discrepancy. Strong male-biased dispersal may explain the difference between the markers: this has been observed on some mammalian species (e.g., Nyakaana and Arctander, 1999; Escorza-Treviño and Dizon, 2000). However, the home range of R. japonica may not be significantly different between sexes (Osawa and Katsuno, 2001) and there is little information about whether dispersal is sex biased in R. japonica. Therefore, we found strong malebiased dispersal to not be a persuasive hypothesis. An alternative explanation may be the different timescale over which different genetic markers accumulate mutations. For instance, gene flow may be high enough to maintain diversity of the microsatellite markers, which have the quadruple effective population size of the mtDNA marker. It is difficult to reject or accept each of these hypotheses without any previous demographical information.
Table 2.
BayesASS results showing the migration rate (m) with 95% confidence intervals (CI) in parentheses and average m for each parameter set.
Table 3.
The estimated numbers of migrating individuals calculated as the product of the average migration rate (m) and egg mass census (Nc).
The genetic structure derived from the clustering analysis based on microsatellite markers was nearly identical to that based on mtDNA (Fig. 1). The results of microsatellite markers suggest that large genetic groups, such as A, B, and E, comprise the main gene pool in this area. The range of these three groups include several stable breeding sites at which we found >50 egg masses every year (e.g., sites 2, 5, 7, 11, and 12; Table 1). The borders of these genetic groups coincided with the location of the main roads and a cement-walled river (Fig. 1iii). Sites 6–8 were recognized as one group, with borders delineated by the groups as well as by the main roads and cement-walled river. However, in the mtDNA results, sites 1–5 and 11–13 were recognized as one large genetic group, whereas they were divided into two clusters in the microsatellite results. It is generally believed that microsatellite markers reveal more fine-scale population structures because of the greater number of alleles (e.g., Zink et al., 2008). This explanation could be also adapted to the differences of groupings at site 9 and 10 between mtDNA and microsatellite markers. These results implied that there would be some factors restricting gene flow for enough time to form different genetic structuring detectable by microsatellite markers. Therefore, while the recent construction of major roads and cement-walled river were possibly one of the main causes of gene flow restriction, it is possible that some genetic structure in these populations predates these developments.
Gene flow
In the BayesASS results, migration rate (m) values were higher within each genetic group than those between different genetic groups. On the other hand, estimated migration rate were relatively higher from group A to B, E to B, E to C and E to D than other directions (Table 2). This trend was also found for the product of m and Nc (Table 3). However, the 95% confidential intervals of migration rate from group E to B and E to C were overlapped with zero, and both migration rates were estimated with relatively large standard errors (Table 2). Therefore, the amount of migration from groups E to B and E to C seems to be uncertain. The migration rates from group A to B and E to D implied recent gene flow across major roads has occurred even after traffic became increasingly heavy, because the BayesASS results are considered to reflect migrations within the past several generations (Wilson and Ranala, 2003). The traffic increased on the roads after rapid economic growth in the 1960s. The reproduction age of R. japonica is not clear, but a previous study suggested that reproducing individuals are mainly aged 1 to 2 years (Marunouchi et al., 2002). This implies that at least 20–30 generations have passed since the traffic became busy as nowadays. Therefore we concluded that group C was isolated from other genetic groups, because group C shows lower m values with other genetic groups in both directions (except for from group E which we discussed and assumed as uncertain result as mentioned before). On the other hand, the directional gene flow still exists from A to B and E to D.
It was difficult to interpret these complicated results, especially considering that the genetic groups from our analyses of genetic structure implied restricted gene flow along the major roads and cement-walled river, whereas the present results implied recent migrations from group A to B and E to D despite major road has become busy. However, if the number of migrating individuals was sufficient to maintain gene flow, the genetic structure should have recombined into one large gene pool, thus appearing as one cluster such as combined cluster I and II in STRUCTURE in the clustering analysis (Frankham, 2002). Therefore, although we considered directional migration from group A to B and E to D were existed in recent years, it was not enough to maintain or to recover one large gene pool. Some local farmers told us that sites 6 and 7 were partially maintained for brown frogs and satoyama landscape recovery on several years ago. We also checked the past satellite images on the site 10, then we confirmed it seemed paddy field would be modified and a new large puddle for rice seedlings has been made on site 10 in the 2000s. These local activities would possibly help to explain the partial recovery of recent migration from group A to B and E to D suggested by BayesASS analysis.
Implications for conservation research and management
In this research, we have revealed the current state of gene flow and genetic structure in a population of the frog in a satoyama landscape in Japan, and discussed historical changes that may have generated these patterns. Based on these results, we should change conservation measures according to the circumstances, depending on whether gene flow has been restored or is still restricted (Lesbarrères e t a l., 2 014). P artial m aintenance at breeding sites 6 and 7 seems to work well but is not sufficient. F or example, creating a new breeding site would be a considerable conservation measure to mitigate the restricted gene flow. However, to determine the location of a new breeding site, we must first understand whether roads and rivers really act as barriers to gene flow and, if so, we should clarify which of the two has stronger influence than another. Landscape genetic techniques, such as least-cost path analysis, are useful tools that can be used to verify how the road and/or cement-walled river strongly affect to restriction of gene flow (Cushman, 2006; Wang et al., 2009; Anderson et al., 2010; McCartney-Melstad and Shaffer, 2015).
Acknowledgments
We thank Ms. Ishiyama and Mr. Furukawa for their help with the sampling and technical assistance. We especially thank Dr. Nakano for his helpful and constructive comments.
