Role of individual dispersal in genetic resilience in fluctuating populations of the gray‐sided vole Myodes rufocanus

Abstract Population densities of the gray‐sided vole Myodes rufocanus fluctuate greatly within and across years in Japan. Here, to investigate the role of individual dispersal in maintaining population genetic diversity, we examined how genetic diversity varied during fluctuations in density by analyzing eight microsatellite loci in voles sampled three times per year for 5 years, using two fixed trapping grids (approximately 0.5 ha each). At each trapping session, all captured voles at each trapping grid were removed. The STRUCTURE program was used to analyze serially collected samples to examine how population crashes were related to temporal variability, based on local‐scale genetic compositions in each population. In total, 461 and 527 voles were captured at each trapping grid during this study. The number of voles captured during each trapping session (i.e., vole density) varied considerably at both grids. Although patterns in fluctuations were not synchronized between grids, the peak densities were similar. At both grids, the mean allele number recorded at each trapping session was strongly, positively, and nonlinearly correlated with density. STRUCTURE analyses revealed that the proportions of cluster compositions among individuals at each grid differed markedly before and after the crash phase, implying the long‐distance dispersal of voles from remote areas at periods of low density. The present results suggest that, in gray‐sided vole populations, genetic diversity varies with density largely at the local scale; in contrast, genetic variation in a metapopulation is well‐preserved at the regional scale due to the density‐dependent dispersal behaviors of individuals. By influencing the dispersal patterns of individuals, fluctuations in density affect metapopulation structure spatially and temporally, while the levels of genetic diversity are preserved in a metapopulation.

reaching peak density, populations are almost exterminated, or crash, within a few months. However, true extinction does not occur because those populations recover to reach peak density again in several years. With such fluctuations, population genetic diversity may decrease due to repeated bottlenecks (Frankham et al., 2017).
However, high levels of population genetic diversity are frequently reported (e.g., Berthier et al., 2005Berthier et al., , 2006Ehrich et al., 2009;García-Navas et al., 2015;Gauffre et al., 2014;Pilot et al., 2010). In their review of northern population cycles, Norén and Angerbjörn (2014) outlined nine predictions covering the direct and genetic feedback consequences of population cycles on genetic variation and population structure and verified them with empirical evidence. They concluded that although genetic variation in northern cyclic populations is generally high and the geographic distribution and amount of diversity are usually suggested to be determined by various forms of context-and density-dependent dispersal exceeding the impact of genetic drift, the signatures of microevolutionary processes such as genetic drift and selection are weaker and obscured by densitydependent dispersal. In fluctuating vole populations, it is thought that genetic variation is maintained through a combination of ecological and evolutionary processes. These include (1) differences in dispersal patterns throughout fluctuations, (2) crash phases with a short duration that are accompanied by weak genetic drift, and (3) the rapid accumulation of new alleles through mutation or immigration (Gauffre et al., 2014;Norén & Angerbjörn, 2014). Berthier et al. (2005) suggested that in the fossorial water vole Arvicola terrestris, dispersal rates are low among isolated remaining populations during periods of low density, resulting in substantial differentiation even among nearby patches. Subsequently, as the populations increase to peak-density levels, the populations expand spatially, facilitating the exchange of individuals between neighboring populations; thus, populations become more similar. Additionally, Berthier et al. (2006) showed that genetic drift, which has a strong influence during periods of low density, and migration, which mainly occurs as population size increases, interact closely to maintain high genetic diversity within a metapopulation (i.e., a population of populations) (Hanski & Gilpin, 1997) at the regional scale (also see Gauffre et al., 2014). Although there is clarity regarding the mechanism by which genetic variation is preserved within a metapopulation, there still is a lack of information regarding how genetic diversity varies with density at the local scale, and how this process is influenced by density-dependent dispersal.
Before and after a population crash, there may be substantial changes in the genetic compositions of individuals that inhabit a specific area. Considering the mechanism by which genetic variation is maintained, and because of breeding among new immigrants in populations with few or no remaining previous residents, local populations that have recovered from a crash would likely acquire different combinations of alleles across loci, compared to a precrash population. Recently, Bayesian inference analyses have been performed in many studies to spatially delineate populations (e.g., Ouchene-Khelifi et al., 2018;Pelletier et al., 2011;Priadka et al., 2019;Viengkone et al., 2016;Wang et al., 2017). Bayesian inference analysis of samples that are collected serially at a fixed site will be likely to show drastic changes in the cluster compositions of individuals, reflecting density crashes.
Individuals of both sexes disperse from their natal sites when they reach sexual maturity (Ishibashi & Saitoh, 2008b;Saitoh, 1995).
Male-biased dispersal is well-known. Within a 3-ha enclosure, approximately 70% of females reproduced within two homerange lengths of their natal sites, whereas more than 80% of the males bred more than 33 m away from their natal sites within their birth year, with a median distance of 103.5 m (interquartile range, 56.9-135.0 m) recorded (Ishibashi & Saitoh, 2008b). During the non-breeding season (i.e., in winter), voles had smaller ranges than during the breeding season; most seemed to nest communally (Ishibashi et al., 1998;Saitoh, 1989). Vole abundance usually reaches a maximum at the end of the autumn breeding season.
Peak densities in areas with the highest vole abundances can reach 100-200 or more individuals per hectare . A time-series analysis of long-term census data in Hokkaido grouped local populations into four types based on their fluctuation patterns, with the most common type showing high fluctuation among years, a ratio of maximum to minimum density of >10, and peaks at intervals of 3-4 years (Saitoh, 1987). Cyclicity is obvious among populations located in northeastern Hokkaido .
In this study, we examined the variations in microsatellite loci diversity in voles three times per year for 5 years at two fixed sites, from which all captured voles were removed in each trapping session, to address the role of density-dependent dispersal in the maintenance of genetic diversity within fluctuating vole populations. We also used a Bayesian inference application to analyze serially collected F I G U R E 1 The gray-sided vole Myodes rufocanus. In Japan, the vole occurs only on Hokkaido, the northernmost island of Japan samples, to explore the relationship between population crashes and local-scale temporal variability in the genetic composition of each population.

| Study site and sampling protocol
On the Nemuro Peninsula, the breeding season of the gray-sided vole is from mid-April to early October. To accurately identify changes in population density and collect liver samples during the breeding season, voles were captured three times per year (mid-spring/late May, mid-summer/early August, and late autumn/early October) for 5 years (2002)(2003)(2004)(2005)(2006) at two fixed sites on Nemuro Peninsula-grids A and I (located 9.3 km apart; Figure 2). Grid A was located in a shelterbelt and grid I was located in a low-lying and damp area. To investigate parasitic Echinococcus multilocularis in the liver, one of the authors (KT) has collected voles three times per year at these grids since the late 1980s with the permission of the Hokkaido Government. Tissue samples were collected from these voles during 2002-2006 for DNA extraction. During each trapping session, 50 Sherman-type traps were arranged in a 5 × 10 grid pattern at an interval of 10 m (~0.5 ha) with a handful of oatmeal in each trap for three nights (i.e., a total of 150 trap-nights per session). Under this sampling protocol, almost all weaned gray-sided voles in the grid would be captured, allowing us to estimate vole density per 0.5 ha at the study sites. Upon capture, voles were euthanized by cervical spine dislocation. At the laboratory of the Nemuro Subprefectural Bureau of the Hokkaido Government, a few small pieces of liver per vole were fixed in 99.5% ethanol after biopsy. Liver samples were stored at −80°C in Sapporo, Hokkaido until DNA extraction. To determine whether similar density fluctuation occurs around the grids, we concurrently performed additional trapping in the areas surrounding grids A and I, during which voles were captured using 49 or 50 live traps on three nights; these voles were released shortly after capture.
Because of a cold sea current (i.e., the Kuril Current), the eastern part of Hokkaido is often enveloped in sea fog from early summer to early autumn. The peninsula has been impacted by anthropogenic activities since the late 19th century, with the creation of many pastures.
Dwarf bamboo Sasa nipponica, the main habitat of gray-sided voles in eastern Hokkaido, is distributed in various places, such as around pastures, along streams, and in shelterbelts with deciduous shrubs (e.g., Alnus hirsuta) as undergrowth. Although the middle of the peninsula has been urbanized, most of the peninsula remains vegetated

| Genotyping
Genomic DNA was extracted from each liver sample using the conventional phenol-chloroform method (Sambrook et al., 1989).
Polymerase chain reaction (PCR) amplification was performed at eight F I G U R E 2 (a) Location of Nemuro Peninsula, Hokkaido, Japan. (b) Locations of two study sites. Areas with some vegetation are colored red. This falsecolor composite image was created with QGIS 3.14.0-Pi (QGIS.org, 2020) using Landsat-8 images. (c, d) Composite images of aerial photographs and shaded-relief maps for grids A and I, which were constructed using QGIS. Squares indicate the positions of each trapping grid (0.5 ha per grid). Aerial photographs and shadedrelief maps were obtained from the Geospatial Information Authority of Japan eight-strip tubes were used and always centrifuged before the caps were opened. Six loci (i.e., MSCRB-01, -04, -06, -07, -11, and -13) were amplified as reported previously (Ishibashi & Saitoh, 2008a,b;Ishibashi et al. 1995Ishibashi et al. , 1997 with an ABI 2720 Thermal Cycler system (Applied Genotyping was performed with an ABI PRISM 310 Genetic Analyzer (Applied Biosystems) using GeneScan Analysis 3.1.2 and Genotyper 2.5 software (Applied Biosystems).

| Genetic characteristics
Basic information such as the number of observed alleles and heterozygosity level was calculated using GenAlEx 6.51b2 (Peakall & Smouse, 2006. Levels of genetic differentiation among trapping sessions (F ST ) were also calculated using GenAlEx; data were used from trapping sessions in which five or more voles were captured. Allelic richness, a measure of the number of alleles within samples that considers the differences in sample size (Mousadik & Petit, 1996), was calculated using the hierfstat package (Goudet, 2005) in R 3.6.2 software (R Development Core Team, 2019). To ensure that the heterozygosity of all individuals was measured on an identical scale, we standardized individual heterozygosity by calculating the proportion of heterozygous typed loci relative to the mean heterozygosity of typed loci (Coltman et al., 1999). Using the web-based version of Genepop (Raymond & Rousset, 1995;Rousset, 2008

| Temporal changes in cluster composition
To investigate temporal changes in cluster composition among individuals, Bayesian clustering analysis of serially collected samples was performed using STRUCTURE 2.3.4 (Falush et al., 2003;Hubisz et al., 2009;Pritchard et al., 2000), which was also used to estimate the population membership of each individual. The proportions of clusters of admixed individuals were estimated using the admixture model. The models used in this program version could detect structure at lower levels of divergence, compared to the original models, with the inclusion of information regarding sampling locations (Hubisz et al., 2009). Because the presence of null alleles may affect analytical outcomes, all loci suspected to have null alleles were excluded from our analyses. The number of clusters (K) varied from one to eight. Each run, replicated 100 times, comprised a burn-in period of 10 4 and 2 × 10 4 iterations with an admixture model under default settings, with the exception of option LOCPRIOR, in which the order of trapping sessions was considered to represent sampling location. Using the R package pophelper 2.3.0 (Francis, 2017), the most probable number of clusters was inferred from ΔK plots (Evanno et al., 2005); plots of membership coefficients (q values) were drawn for the inferred K value, for which individual assignments to clusters were determined using the CLUMPP 1.1.2 program (Jakobsson & Rosenberg, 2007).

| Fluctuations in density
The number of voles captured at each session (i.e., vole density) varied considerably at both grids, despite constant trapping efforts

| Genetic diversity
Genotypes were determined based on the eight microsatellite loci for all collected samples (461 and 527 at grids A and I), with the exception of a single animal, which had a missing value for one locus, MSCRB-11 (Table 1; Appendices S1 and S2). At each grid, each locus examined exhibited high variability; allelic diversity was similar between grids. Each locus had 11 or more alleles. The  grid I (Table 1; Appendices S1 and S2); however, the presence of null alleles was not suspected in some trapping sessions for these loci.   F I G U R E 5 Relationship between density and allelic richness, which is a measure of the number of alleles within samples that considers the differences in sample size (Mousadik & Petit, 1996

| Relationship between density and allelic diversity
As density increased, the average number of alleles observed at the eight loci increased, after the number of alleles at loci suspected to have null alleles had been increased by one for each trapping session ( Figure 7). However, as vole density approached peak levels, the rate of increase gradually slowed. At both grids, we did not record all alleles for a locus in a single sample, even at peak density, except in one instance-all 11 alleles of MSCRB-01 were recorded in session 9 (October 2004) at grid I.
To describe the nonlinear relationship between density and allelic diversity, the nls command in R was applied to obtain nonlinear least-squares estimates of parameters for three assumed models after Crawley (2015) and Logan (2010) (Appendix S5). The appropriateness of the models was then assessed by comparing their Akaike Information Criterion (AIC) values (Burnham & Anderson, 2002).
Although the coefficients of determination (R 2 ) were higher than 94% for all models, one model, Y = a + b × ln (X), where Y is the mean number of observed alleles and X is the density, produced both the highest R 2 values and lowest AIC values for both trapping grids To investigate the timing of the appearance of new alleles, we counted alleles that had newly appeared in the first five sessions after a crash, in which each allele was counted only at the first appearance.
Because two crash events occurred in grid I during the study period (i.e., sessions 1 and 10), three time-series data sets were used for this analysis: sessions 11-15 for grid A and sessions 2-6 and 11-15 for grid I. This analysis revealed that more new alleles appeared during the early increasing phase immediately following a crash (Figure 8).

| Temporal changes in cluster composition
First, STRUCTURE analyses were performed for all samples from both trapping grids, to analyze samples as if they had been collected at different sites. The highest ΔK value was observed at K = 2 (Figure 9). In the bar plot representing cluster assignment among individuals (hereafter referred to as a "STRUCTURE plot"), two clusters, each from one grid, were distinctly separated, although genotype data had only been used for three loci (MSCRB-01, -04, and -13; Figure 10). Additionally, at each grid, temporal variability in cluster composition was relatively small among trapping sessions.
When STRUCTURE analyses were performed separately for To clearly demonstrate temporal changes in cluster compositions, the frequency distributions of q values were compared before and after a crash in each grid, such that we used only the q values corresponding to each individual's cluster assignment, that is, its largest q value (Sacks et al., 2005). For grid A, before a crash (sessions 1-6), two of the three clusters accounted for higher proportions, and the frequency of the remaining cluster was very low (Figure 13a).
By contrast, after a crash (sessions 11-15), all individuals in the latter cluster had the largest q values (Figure 13b). Similarly, for grid I, one of two clusters accounted for almost all of the q values before a crash (sessions 2-9) and after a crash (sessions 11-15), the other cluster accounted for a large amount (Figure 13c,d). Although similar critical changes in cluster composition occurred in both grids, the grids exhibited distinct frequency distribution patterns; throughout the study period, the largest class (0.9-1.0) of q values was most frequent at grid I, whereas values in the lower classes (0.5-0.8) were more frequent at grid A.

Mean number of alleles ± SE
Grid A Grid I

| Resilience in genetic diversity
Our results revealed that genetic diversity was maintained even as the population sizes of gray-sided voles fluctuated. Although levels of heterozygosity, allelic richness, and standardized individual heterozygosity were continuously high at both sites (Figures 4-6), allelic diversity was strongly, positively correlated with density (Figure 7). At each grid, if the density was similar, the mean number of alleles observed was also similar, irrespective of trapping sessions. Large fluctuations in density have been observed at both study sites since the late 1980s, and we suspect that such fluctuations have been happening for decades. Gauffre et al. (2014) found that, in the common vole Microtus arvalis, genetic diversity fluctuates at the local scale, while it remains constant at the regional scale. In gray-sided vole populations, genetic diversity may be well-preserved within metapopulations at the regional scale, whereas fluctuations in genetic diversity occur continuously at the local scale. Here, we observed that allele numbers were restored within a short period of time. The mutation rate for microsatellite loci is low (<7 × 10 −3 per locus per gamete per generation; Brinkmann et al., 1998); thus, mutation cannot explain the observed accumulation of new alleles, because there were only a few generations of voles between periods of increasing and peak densities (Berthier et al., 2006). The rapid accumulation of new alleles after periods of low diversity must be due to individual dispersal.
In the present study, significant departures from Hardy-Weinberg equilibrium were observed in some trapping sessions at both grids.
The Hardy-Weinberg equilibrium may not be common in wild populations with high levels of admixture caused by dispersal (Waples & Gaggiotti, 2006). In addition to the presence of null alleles, active individual dispersal might have caused a slight shift from random to nonrandom mating at the study sites. Gauffre et al. (2014) suggested that common voles travel over longer distances when density increases. Therefore, in the populations of gray-sided voles studied here, individual dispersal must be the main contributor to the maintenance of genetic diversity. Although the number of alleles observed was correlated with density, this correlation was not linear; the rate of increase gradually declined with density ( Figure 7). Furthermore, the alleles observed during the 5-year study period did not always appear at the same time. These results imply that new genetic variants appear in populations mainly due to longdistance dispersal during the early increasing phase when densities are very low; as density reaches high levels, the dispersal distance gradually decreases. Short-distance dispersal appears to occur frequently under mid-to high-density conditions, although little is known regarding individual dispersal behavior during transient and settlement stages in arvicoline rodents (Le Galliard et al., 2012). Our observation that more new alleles appeared during the early increasing phase supports this implication (Figure 8).
Notably, the relationship between density and allelic diversity was very similar between trapping grids (Figure 7). This implies F I G U R E 8 Numbers of newly appearing alleles after each crash. Alleles newly appearing in the first five trapping sessions after a crash were counted; alleles were counted only at their first appearance. Mean numbers for eight microsatellite loci are shown. Vertical bars indicate SEs. Three time-series data sets were used for this analysis (sessions 11-15 for grid A; sessions 2-6 and 11-15 for grid I). Null alleles were not considered acted in a similar manner to maintain genetic diversity within each metapopulation. From our human perspective, the trapping grids appeared to be located in habitats that are quite different, such that grid A was located in a shelterbelt and grid I was located in a lowlying and damp area. However, these differences in habitat did not appear to affect vole dispersal behavior. Because we only surveyed genetic diversity at two sites, it remains unknown whether similar relationships can be found at other sites near the grids.

| Dynamic changes in genetic composition
Because of moderate genetic differentiation among samples at each grid, STRUCTURE worked well for inferring the numbers of clusters F I G U R E 1 0 STRUCTURE plot for all data from grids A and I (K = 2). Numbers 1-15 and 16-30 above bars correspond to trapping sessions 1-15 at grid I and trapping sessions 1-15 at grid A, respectively. Genotype data for three microsatellite loci (MSCRB-01, -04, and -13) were used. At grid I, no voles were captured in session 10 (May 2005) F I G U R E 11 ΔK plots for (a) grid A, for which genotype data at three microsatellite loci (MSCRB-01, -04, and -13) were used, and (b) grid I, for which genotype data at six microsatellite loci  were used. For each grid, the loci used for analysis were limited due to the possible presence of null alleles at other loci PEAK (Latch et al., 2006). Our STRUCTURE analyses revealed that patterns of cluster composition differed between grids (Figure 12).   (Figure 2c).

At grid
Accordingly, cluster compositions at this grid varied throughout the study period.
When the analysis was performed using all data from both grids, the cluster compositions exhibited considerable differences. Two clusters were detected (K = 2) and voles were sorted into distinct groups based on the grid in which they were captured, independent of trapping session ( Figure 10). In the STRUCTURE plot, critical changes were not observed at low density for either grid; temporal variability in cluster composition was small among trapping sessions at each grid. Because there is an urbanized area between the two grids (Figure 2b), the rate of gene flow between grids may be low.
Based on differences in the demographic, environmental, and historical processes that have led to the current genetic structure at each site, levels of genetic organization are also likely to differ between sites (Meirmans, 2015). areas, thus helping to preserve genetic variation and counteract differentiation among subregions. However, these are snapshot studies of populations, and genetic diversity is not examined regularly in lemmings (Ehrich et al., 2020). Long-term genetic studies are needed to address the role of density-dependent dispersal in the maintenance of genetic diversity in cyclic lemming populations.

| CON CLUS ION
In this study, we found that density fluctuated at similar levels between the two study sites, although patterns of fluctuations were not synchronized. Allelic diversity was well-preserved at each site, and the relationship between density and allelic diversity was similar between sites. In their review of population cycles, Norén and Angerbjörn (2014) concluded that the signatures of genetic drift and selection on population genetic diversity are weaker and obscured by density-dependent dispersal. The present study adds to the growing body of evidence that dispersal usually overshadows the impact of other microevolutionary processes in cyclic populations. Rikalainen et al. (2012) proposed that in addition to increased individual dispersal and the resulting accumulation of new alleles during peak phases, a constant presence and large effective population size facilitate the maintenance of high genetic diversity within fluctuating populations. In the gray-sided vole, it is unknown how fluctuations in density affect metapopulations in which genetic diversity is totally preserved.
Because a few voles were captured from local populations during lowdensity phases, the effective population size of metapopulations at both study sites must be large. In addition to landscape features, which can affect connectivity among populations, fluctuations in density affect the spatial and temporal structure of metapopulations; genetic diversity is preserved within a metapopulation through density-dependent individual dispersal. Further studies are needed to clarify how fluctuations in density influence metapopulations at the regional scale.

ACK N OWLED G M ENTS
We thank the Nemuro Public Health Office, Department of Health and Welfare of the Hokkaido Government for providing us with laboratory space to dissect voles.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sampling locations and microsatellite genotypes were submitted to Dryad: https://doi.org/10.5061/dryad.mcvdn cjzf.