Regional and local patterns of genetic variation and structure in yellow‐necked mice ‐ the roles of geographic distance, population abundance, and winter severity

Abstract The goal of this study, conducted in seven large woodlands and three areas with small woodlots in northeastern Poland in 2004–2008, was to infer genetic structure in yellow‐necked mouse Apodemus flavicollis population and to evaluate the roles of environmental and population ecology variables in shaping the spatial pattern of genetic variation using 768 samples genotyped at 13 microsatellite loci. Genetic variation was very high in all studied regions. The primal genetic subdivision was observed between the northern and the southern parts of the study area, which harbored two major clusters and the intermediate area of highly admixed individuals. The probability of assignment of individual mice to the northern cluster increased significantly with lower temperatures of January and July and declined in regions with higher proportion of deciduous and mixed forests. Despite the detected structure, genetic differentiation among regions was very low. Fine‐scale structure was shaped by the population density, whereas higher level structure was mainly shaped by geographic distance. Genetic similarity indices were highly influenced by mouse abundance (which positively correlated with the share of deciduous forests in the studied regions) and exhibited the greatest change between 0 and 1 km in the forests, 0 and 5 km in small woodlots. Isolation by distance pattern, calculated among regions, was highly significant but such relationship between genetic and geographic distance was much weaker, and held the linearity at very fine scale (~1.5 km), when analyses were conducted at individual level.

time. Species that undergo cyclic density fluctuations are extreme cases of demographic instability. Computer simulations showed that IBD pattern depends strongly on how populations differ in size and how they fluctuate, but it does not depend on the gene flow frequency or pattern (Björklund, Bergek, Ranta, & Kaitala, 2010).
In the cyclic species which undergo demographic fluctuations over shorter periods of time (i.e., less than 10 generations), the spatial genetic structure at a particular time may be shaped by both presenttime and recent-past demographic events. During the phase of low density, genetic drift is expected to reduce genetic diversity and substantially increase genetic differentiation among isolated subpopulations (Leblois, Rousset, & Estoup, 2004). In a low phase of population fluctuations, spatial aggregations of small mammals form temporary metapopulation structures (e.g., Lima, Marquet, & Jaksic, 1996). Dispersal of individuals among suitable patches of habitat is an essential ingredient of metapopulation dynamics. It was also reported by Szacki and Liro (1991) and Szacki, Babińska-Werka, and Liro (1993) that rodents move longer distances in heterogeneous forest habitat and suburban habitat compared with more homogeneous forest. Increased mobility was suggested as a behavioral adaptation to fragmented landscape.
Studies conducted by Ehrich, Yoccoz, and Ims (2009) on two numerically dominant vole species of northern Fennoscandia, red voles Clethrionomys rutilus, and gray-sided voles C. rufocanus, showed that high-amplitude density fluctuations increased gene flow and genetic variability in vole population. It was also noted that spring densities, when genetic bottleneck could be observed, had no effect on differentiation or local genetic diversity. In certain cases, subdivided population with frequent local extinctions may retain more genetic diversity because founders of newly colonized patches tend to be a mixture from different subpopulations rather than from one larger source population, which increase effective population size (N e ). The fraction of genetic variation lost during the bottleneck is a function of the population growth rate (Nei, 1975). Populations, which recover quickly after the bottleneck, lose little genetic variation even after severe density crash. Moreover, overlapping generations can minimize the effect of environmental fluctuations on population size and thus on the level of maintained genetic variability (Gaggiotti, 2003). Adult individuals reproduce several times throughout their lives; therefore, the genetic variation observed in a given cohort is more likely to be transferred to future generations than in case of organisms with discrete generations (Ray, 2001).
Yellow-necked mouse ( Figure 1) is one of the most numerous forest rodent species in Central European temperate forest zone (Niedziałkowska, Kończak, Czarnomska, & Jędrzejewska, 2010). It prefers deciduous forest (especially oak-lime-hornbeam stands) and generally feeds on tree seeds with slight supplementation of its diet with invertebrates (Butet & Delettre, 2011;Gębczyńska, 1976;Hansson, 1985). Recent isotope analyses of rodent have shown that yellow-necked mice are strongly associated with mast of deciduous trees (>80% of diet; Selva, Hobson, Cortés-Avizanda, Zalewski, & Donázar, 2012). Population dynamics of yellow-necked mice is largely shaped by seed crops of deciduous trees. In eastern Poland, heavy crops (usually synchronous in common oak Quercus robur, hornbeam Carpinus betulus, and maple Acer platanoides) occur at 6-to 9-year intervals and trigger winter breeding in mice, their outbreak in the following year, and crash in the next year (Pucek, Jędrzejewski, Jędrzejewska, & Pucek, 1993). Two years of outbreakcrash are followed by 4-7 years of moderate densities with winter mortality averaging 86% of autumn numbers of mice (Pucek et al., 1993).
Genetic studies on yellow-necked mice are scarce and mostly restricted to phylogeography (Michaux, Libois, & Filippucci, 2005;Michaux, Libois, Paradis, & Filippucci, 2004). The only studies that aimed at including landscape features in explaining the genetic structure of yellow-necked mice had been focused on evaluating the direct effect of natural and anthropogenic barriers on defined a priori population and covered a small geographical scale (Gortat et al., 2010;Kozakiewicz et al., 2009;Rico, Kindlmann, & Sedláček, 2009).
The aims of this study were to infer genetic structure in yellownecked mice (Figure 1) population and to evaluate the role of environmental and population ecology variables in shaping the observed pattern of genetic pattern of the species in northeastern Poland.
To address these questions, we combined microsatellite, ecological (population density and habitat structure in study sites), and environmental (geographic distance, climatic conditions) data. Given the ecological characteristic of the species, we were particularly interested in testing the hypotheses that: (a) yellow-necked mouse population is characterized by high level of genetic variation structured in a complex hierarchical way, (b) mouse abundance has profound effect on local population structure, and (c) genetic structure is related to ecological and environmental conditions in the study sites.

| Study area
The study was conducted in seven large woodlands located in the lowlands of northeastern Poland: Augustów (AUG), Białowieża The landscape of the region has been shaped by glaciations (mainly the Riss, 310,000 to 130,000 yBP, and the Würm, 70,000 to 10,000 yBP). It is mainly a plain, with some belts of frontal and moraine hills and numerous postglacial lakes. The altitude is between 25 and 312 m a.s.l. Forest cover is slightly higher (29.9%) than the mean for the country (28.7%) (Statistical Yearbook of the Republic of Poland, 2005). The climate is transitional between continental and Atlantic types. Annual precipitation is 550-700 mm. Snow cover lasts for 76-90 days, and the growing season lasts for 180-200 days (Kondracki, 1972). During the study period, the mean temperature in January was −7°C and that in July 22°C (Supporting Information Table S1).
Most of the sampling sites were located in managed forests administrated by the Olsztyn and Białystok Regional Forest Directorates. Minor parts of Białowieża and Augustów Forests are protected as national parks, and from two to 24 small nature reserves are located in each of the seven woodlands. The dominant tree species were Scots pine (Pinus silvestris), Norway spruce (Picea abies), common oak, silver birch (Betula pendula), white birch (B. pubescens), and common alder (Alnus glutinosa). The woodlands differed in size (from 126 to 1,600 km 2 ), productivity, and mean age of tree stands.
The Białowieża Forest is the oldest and the best preserved woodland in the whole study area (Sokołowski, 2006  At each site, two types of livetraps were placed: small wooden traps (Dziekanów type, 16.5 × 9.5 × 8 cm) for rodents, large wooden boxes (11 × 11 × 24.5 cm) for dormice (applicable also for mice). Six to 10 small traps were placed on the ground, five large boxes were located on shelves stuck by sticky tape to tree trunks at a height of about 2 m. The trapping and sampling was conducted according to the protocol described by Niedziałkowska et al. (2010). In total, 820 individuals of yellow-necked mouse were captured at 176 trapping sites within a period of 2004-2008. Tissue samples for DNA analysis (tail or ear clipping) were placed in a tube with 96% alcohol and stored in −20°C upon DNA extraction. All capture procedures were in accordance with Polish law and were accepted by the Local Ethical Commission in Białystok (permissions nos 07/2004, 15/2006, 65/2007). Trapping site locations were approved by the administrators of the area. All animals were released at the place of capturing.

| Habitat structure and climatic data of trapping sites as analyzed with GIS tools
Precise location of each trapping site was identified using GPS GARMIN60CSx. Complete botanical description of forest floor vegetation, undergrowth and tree stand was collected for each trapping site. Based on this description, each plot was classified according to forest habitat type as coniferous, mixed or deciduous forest. In program ArcGIS (ESRI), a buffer was defined around each trapping site with a radius equal to the mean distance between all trapping sites (~1 km). Next, Corine Landcover 2006 (CLC2006) maps were used to measure the percentage cover of each land use category (coniferous, mixed and deciduous forests, meadows and arable lands, waters and other) within a defined buffer zone (Supporting Information Table   S1). Moderate Resolution Imaging Spectroradiometer (MODIS) was used to collect land surface temperature (LST, Wan, Zhang, Zhang, & Li, 2003) with high spatial (grid 1 × 1 km) and temporal (four measurements per month for each trapping plot for 2004-2008) resolution. Collinearity between the explanatory variables was assessed by calculating pairwise Pearson's correlation coefficients and only uncorrelated variables (R < |0.5|) were selected for further steps of analysis: mean temperature of January (T JAN ), mean temperature of July (T JUL ), and share of deciduous-mixed forest (DMF, the most suitable habitat for yellow-necked mouse) in 1-km buffer zone around the trapping sites.

| Laboratory analysis
Genomic DNA was extracted using DNeasy Blood & Tissue Kit (QIAGEN) according to the manufacturer's protocol. Quantity and quality control was performed using spectrophotometer NanoDrop ND-1000. We selected 15 markers that were tested previously on yellow-necked mouse by other researchers and therefore guaranteed amplification success: AF246520, AF246522, AF246523 (Harr, Musolf, & Gerlach, 2000) annealing at 54°C, 58°C or 60°C for 30 s (following the protocol in Makova et al., 1998) and extension at 72°C for 60 s with initial step of denaturation for 15 min at 94°C and final elongation for 3 min at 72°C. A mixture containing 1 μl of diluted PCR product, 0.2 μl of GeneScan 400 ROX size standard (Applied Biosystem) and 10 μl of Hi-Di Formamide (Applied Biosystem) was separated on Applied Biosystems 3100 DNA Analyser. Genotypes were determined with GeneMarker ver. 1.75 (SoftGenetics LLC, 2007). All PCR reactions were carried out using a DNA Engine Dyad Peltier Thermal Cycler (BIO RAD).
A final data set consisted of 768 samples genotyped at 13 microsatellite loci (GTTA1A and GTTD8S were excluded due to poor amplification and detected monomorphism). Detailed information on the number of analyzed samples from each woodland, transect, and year is presented in Supporting Information Table S2.

| Statistical analysis of genetic data
Microsatellite variability statistics, calculated per locus and per geographic region (number of alleles per locus, number of private alleles, observed and expected heterozygosity) were calculated in GenAlEx (Peakall & Smouse, 2006). Allelic richness and inbreeding coefficient (F IS ) were calculated using FSTAT (Goudet, 1995 (Weir & Cockerham, 1984) was calculated in Arlequin ver. 3.11 (Excoffier, Laval, & Schneider, 2005). Mean values of Queller and Goodnight relatedness coefficient in each woodland and the transect were calculated in GenAlEx. We estimated effective population size N e using a bias-corrected version of the linkage disequilibrium method of Waples and Do (2008) with the option of random mating selected using NeEstimator V2.1 (Do et al., 2014). STRUCTURE 2.3.4 was used to infer population structure and assign individuals to subpopulations (clusters) based on individual multilocus genotypes without any prior spatial information (Pritchard, Stephens, & Donnelly, 2000). STRUCTURE was run using 50 5 MCMC iterations, following a burn-in period of 50 4 iterations with admixture model and correlated allele frequencies (Falush, Stephens, & Pritchard, 2003). Number of genetic groups (K) was tested from 1 to 10 with 10 runs completed per each K. It is also known that sampling scheme might significantly affect the assessment of population subdivision (Oyler-McCance, Fedy, & Landguth, 2013;Schwartz & McKelvey, 2009). Therefore, additionally to the whole data set, we ran STRUCTURE analyses solely for mice collected in seven woodlands with model settings of no admixture (it assumes that each individual belongs to a single cluster, data not shown) and admixture (each individual draws some fraction of its genome from each of the K populations) and separately on samples collected in each year to check for the consistency of results among years. Program STRUCTURE HARVESTER v. 06.94 (Earl & Vonholdt, 2012) was used to summarize the results and produce the graph presenting the posterior probability of the data [LnP(D)] and delta K (Evanno, Regnaut, & Goudet, 2005). Subsequently, Principal Component Analyses (PCA) were performed on individual genotypes from 10 geographical regions using the adegenet package (Jombart, 2008) in R (R Development Core Team, 2016). The PCA does not assume Hardy-Weinberg equilibrium (HWE) and is highly efficient in revealing genetic structure in the form of clines, which is more difficult to detect than clusters (Jombart, Pontier, & Durfour, 2009).
We also evaluated the genetic structure results by applying spatially explicit algorithms GENELAND ver. 3.3 (Guillot, Mortier, & Estoup, 2005) and TESS ver. 2.3.1 (Chen, Durand, Forbes, & Francois, 2007), Bayesian clustering programs that incorporate spatial coordinates into analyses. Additionaly, TESS allows to include "dummy points" when samples are not evenly distributed. Such points are not included in the analysis but help to avoid ghost populations. Inferences in GENELAND were performed in a single step as recommended by the author (Guillot, 2008). This approach makes inferences faster and avoids the issue of ghost populations. Thus, the MCMC were ran 50 times with settings as follows: delta.coord 0.05, 200,000 iterations, thinning 200, uncorrelated allele frequency, and 1-15 population. TESS was run using the MCMC algorithm with the admixture model. At first we used three initial starting values for the interaction parameter, which represents spatial interactions (w = 0.6, 0.75 and 0.9) to select the optimal interaction parameter using likelihood (w = 0.6). Next, 10 independent simulations were performed for each K value (2-12) using 50,000 iterations and a burn-in period of 10,000 iterations to identify, which K values produced the highest likelihood runs (K max ). Finally, we conducted 100 simulations for each K max and calculated membership probabilities from the 20 highest likelihood simulations for each K max (Chen et al., 2007).
To test for genetic differentiation between defined spatial groups, pairwise F ST was calculated in Arlequin ver. 3.11 (Excoffier et al., 2005). The levels of significance for multiple tests were adjusted by the sequential Bonferroni method (k = 6) (Rice, 1989).

| Association of genetic variation with ecological factors
Isolation by distance was assessed using linear Mantel test on defined regions and individual level in R package vegan (Oksanen et al., 2011). First, F ST coefficient values calculated among 10 geographic regions were plotted against the mean geographic distances.
Significance was determined using 999 permutations.
Fine-scale genetic structure and dispersal distance were analyzed using spatial autocorrelation methods implemented in GenAlEx, performed separately for individuals from each region (woodland or transect). Pairwise genetic distance and geographical distance matrices were used to calculate autocorrelation coefficient (r) for each distance class presented as a correlogram. The geographical distances were calculated as Euclidean distances between samples and the analyses were performed with an equal distance classes set at 1 km. Numbers of the calculated distance classes differed between regions due to differences in the spatial distribution of individuals in each region. For each analysis, we used 999 permutations to test the hypothesis of no spatial genetic structure, and 1,000 bootstraps to estimate 95% confidence intervals for r in a given geographical distance (Peakall, Ruibal, & Lindenmayer, 2003). Subsequently, we evaluated the association between environmental variables (mean temperature of January, mean temperature of July, and percentage cover of deciduous-mixed forests in 1-km buffer zone around trapping sites) and genetic variation measured as probability of each individual (n = 752) to be assigned to Cluster 1 (based on STRUCTURE with K = 2, see Section 3). We applied the generalized linear mixed models (GLMM) for binomial data using the "glmer" function implemented in the lme4 package (Bates, Maechler, Bolker, & Walker, 2015). We used mixed-model framework with observation-level random effect because our data expressed excess variation which had led to overdispersion problem. An information theoretic approach with a second-order correction for small sample size (AIC c ) (Burnham & Anderson, 2002) implemented in R package MuMIn (Bartoń, 2015) was applied to select the most parsimonious model, which best explained the observed spatial pattern in genetic variation of mice. Statistical analyses were performed in R (R Development Core Team, 2016).

| Genetic variability and structure of the yellownecked mouse population
All 13 loci included in statistical analyses were highly polymorphic with number of alleles per locus ranging from 15 to 32 in the whole population of yellow-necked mice in northeastern Poland (Supporting Information Table S3). Allelic richness, a parameter that controls for differences in a sample size, was similarly high in all regions ranging from 9.38 in Pisz Forest to 10.55 in Białowieża-Mielnik Transect (Table 1)  Inferences in program STRUCTURE indicated major division of the yellow-necked mouse population into two genetic clusters followed by reasonable divisions up to K = 5 (Figure 3). Although LnP(D) continued to increase until it reached a plateau at K = 5, ΔK (Evanno et al., 2005) showed highest support for K = 2 with a slight increase at K = 4 and K = 5 (Supporting Information Table S4). Subdivision at K = 2 was characterized by genetic gradient from north to south, with Cluster 1 (green color in Figure 3 Information Table S5).

| Genetic similarity among mice in relation to geographic distance and population abundance
The population abundance index of yellow-necked mice varied from 1.94 to 9.37 per 100 trapnights (Supporting Information   Table S1). Such a wide range suggests huge differences in population density. Jędrzejewski, Jędrzejewska, and Szymura (1995) have shown that linear increase in abundance index between two and nine yellow-necked mice per 100 trapnights reflects the exponential increase in absolute population density between 5-10 TA B L E 2 Pairwise genetic differentiation F ST among yellow-necked mice in 10 geographical regions (see Figure 2). All values of F ST were statistically significant at p < 0.001 after Bonferroni correction  Table S1) and was positively correlated to the index of mice abundance Forest showed significant autocorrelation within 1-km distance class only, whereas mice from Białowieża Forest exhibited significant rvalues up to 12-km distance classes (Supporting Information Figure   S2). The maximum distance, at which significant autocorrelation was observed, was 14 km for individuals from Augustów-Knyszyn Transect. There was no clear spatial pattern found with high variation observed among geographical regions (Supporting Information Figure S2).  Figure S3). Mean relatedness coefficients and genetic distances calculated for mice captured along transects had moderate values, within the range found in forests. Relatedness level in the first distance class ranged from 0.034 to 0.093 and genetic distance from 0.029 to 0.052 (Supporting Information Figure S5).
In the six forests (PIS excluded due to small sample size), individual genetic distances were significantly positively associated with an abundance index of mice at the first distance class (R 2 = 0.70, p = 0.02), the relationship got weaker at 1-km distance, and disappeared at 2 km ( Figure 7). Dependence of pairwise relatedness on abundance of mice was even higher, with significant relationship at the first and second distance classes (Figure 7). Interestingly, high level of genetic similarity was negatively associated with increasing abundance of mice. The relationship between abundance of mice and their pairwise relatedness/genetic distance was not visible for data collected along three transects (Supporting Information Figure   S5).
In our testing the association between environmental variables and genetic assignment of yellow-necked mice to Cluster 1, the AIC ranking indicated that-among all the combinations of considered submodels-the model with all explanatory variables had the lowest AIC c scores. Therefore, we selected this submodel as a single best model (Table 3). The GLMM showed a negative effect of growing temperatures (both January and July) on the probability of genetic assignment of mice to Cluster 1 (p < 0.001; Table 4, Figure 8). In addition, the prevalence of Cluster 1 decreased with the growing share of deciduous-mixed forests (p < 0.001). Winter severity (approximated by mean temperature of January) had the strongest impact on Cluster 1 spatial distribution (see z-values in Table 4). In general, mice belonging to Cluster 1 prevailed in colder regions and in the boreal-type, coniferous forests.

| Genetic variation and differentiation among local populations
Genetic variation of yellow-necked mice in NE Poland, estimated using 13 microsatellite markers, was very high in all studied regions.  Figure S2 and Table S1 F I G U R E 7 Mean values of Rousset's genetic distances (a rupper panel) and Queller and Goodnight relatedness coefficients (lower panel) calculated for individuals from six forests (PIS excluded due to small sample size) in three distance classes in relation to abundance of yellow-necked mice. Sources of data in Supporting Information Figures S3 and S4 and Table S1 et Rico et al., 2009) and in several vole species (e.g., Aars et al., 2006;Berthier et al., 2005Berthier et al., , 2006Ehrich & Jorde, 2005;Ehrich et al., 2009;Gauffre, Estoup, Bretagnolle, & Cosson, 2008;Gauffre et al., 2014;Plante, Boag, & Bradley, 1989;Redeker et al., 2006;Rikalainen et al., 2012;Vuorinen & Eskelinen, 2005 Greenwood, 1980) one might expect negative F IS values as a result of mating between individuals originating from different source populations (Prout, 1981). Such phenomenon is a common feature in fragmented vole populations (e.g., Ratkiewicz & Borkowska, 2006) where high dispersal rate is negatively correlated with population size. Lack of negative values in our data might indicate that there is no sex bias in dispersal rate in the yellownecked mouse, low number of immigrants are able to reproduce in favorable habitat (e.g., Białowieża Forest; F IS = 0.101) or individuals from adjacent areas that successfully reproduce are genetically similar. It was found in voles trapped in Białowieża meadow that recruitment of individuals from outside did not substantially affect the TA B L E 3 Model selection (based on the AIC c criteria) for the considered GLMMs. The models aimed at assessing the effects of mean temperature of January (T JAN ), mean temperature of July (T JUL ), and the percentage of deciduous-mixed forests (DMF) on the probability of genetic assignment of individual mice to Cluster 1 (see Figure 3; green color) F I G U R E 8 Relationships between mean individual assignment of individual yellow-necked mice to Cluster 1 (based on Structure results) and mean temperature of January, mean temperature of July, and percentage of mixed and deciduous forests in the trapping sites-results of the most parsimonous model presented in Table 4 genetic results as they exhibited highly similar genetic profiles (Pilot et al., 2010).
Deviation from Hardy-Weinberg equilibrium observed in our data is reasonable and results from the assumption of the theoretical model that is rarely met in natural populations. The geographical scale of the study area greatly exceeded the potential for unlimited dispersal of individuals. However, when HWE was tested in a smaller scale, substantially lower number of loci exhibited significant excess of homozygotes. Similar results were found in the yellow-necked mouse population investigated by Kozakiewicz et al. (2009) and Gortat et al. (2010). In their analyses, all samples from Pisz Forest pooled together resulted in two loci deviating significantly from HWE, whereas such deviation was not found when analyses were conducted at the very fine scale.
Genetic differentiation among the studied regions was low but significant. Hierarchical population structure was found. Primal genetic subdivision was observed between the northern and the southern parts of the study area. The obtained results on population genetic structure were highly concordant among applied nonspatial (STRUCTURE, PCA) and spatial (GENELAND, TESS) methods.
Clear genetic distinction of northern subpopulation (Cluster 1, G1) was supported by all programs. The remaining samples with more continuous spatial distribution were less genetically structured, and formed a gradient but with several detected areas of disruption in gene flow, which resulted in a more subtle but still visible differentiation.
Low values of F ST (between geographical regions and genetic subpopulations) observed in this study are in a range of values found in similarly abundant and highly variable small mammal species. Nagylaki (1998) and Hedrick (1999)

| Factors explaining genetic variation pattern
Identification of ecological variables that contribute to population genetic divergence is a primal goal in landscape genetics (Manel, Schwartz, Luikart, & Taberlet, 2003). Although the impacts of various features are not mutually exclusive and spatial genetic pattern usually results from combination of several factors, we made the attempt to evaluate the role of each category of factors separately.
At the larger scale, geographic distance was found to be the main factor shaping a nuclear genetic variation. It was supported by F ST values, PCA plot structure and the gradient of genetic differentiation. Nevertheless, such relationship between genetic and geographic distance was much weaker, and held the linearity at very fine scale (~1.5 km), when analyses were conducted at individual level. Such incoherent results are in high concordance with published data on small mammals (e.g., Trizio et al., 2005). Despite unrealistic assumptions of the theoretical model introduced by Wright (1943), isolation by distance (IBD) is still commonly used by population genetists to describe spatial pattern of genetic variation and such correlation is indeed often found (Jenkins et al., 2010 and reference therein). However, genetic IBD patterns may not necessarily be observed even if dispersal is actually restricted in space (Guillot, Leblois, Coulon, & Frantz, 2009;Leblois et al., 2004;Rousset, 1997). Lack of genetic IBD pattern is usually explained as results of substantial, long-distance dispersal Redeker et al., 2006;Schweizer et al., 2007), dispersal barriers (Ratkiewicz & Borkowska, 2006), large effective population size (Guivier et al., 2011) or migration-drift disequilibrium Francl, Glenn, Castleberry, & Ford, 2008).
Overall high genetic variation of mice in NE Poland and their consistent genetic structure across consecutive years indicated that fluctuations in population size had no serious impact on genetic parameters. Most probably it is due to a large effective population size, which can buffer population response to evolutionary forces. The differences in habitat composition (the share of deciduous forest) among the sampled areas had a profound impact on population density indices, what subsequently influenced kinship and spatial autocorrelation results. The decrease of relatedness in dense population was observed also in the root vole Microtus oeconomus (Pilot et al., 2010) and it was suggested that with increasing density more kin groups occupy the same space. On a contrary, in years with sparse distribution of individuals, one family group can span over larger area what results in higher estimates of mean relatedness. Positive autocorrelation was also found to be highly associated with population density, but it is unclear whether the stronger population structure resulted from dispersal rate or favorable habitat share.
Winter severity was found to be the main factor that explained the spatial distribution of the two main genetic clusters of mice in NE Poland. Although the pattern revealed by microsatellites is believed to reflect the neutral genetic variability, it can be highly consistent with the adaptive variability as has been shown for European wolves Canis lupus (Stronen et al., 2013(Stronen et al., , 2015. It was already found by Wójcik (1993) that-in a population inhabiting Białowieża Forest-mice surviving winter were not a random representation of population in respect of polymorphism in the transferrin locus.
Heterozygotes had a selective advantage in winter survival over homozygotes.

| CON CLUS IONS
In NE Poland, the yellow-necked mouse population exhibited a major division into two genetic clusters. Their spatial distribution was correlated to environmental factors, mainly winter severity. Fine-scale genetic structure was highly influenced by mouse population abundance, which was determined by the share of optimal habitats (deciduous forests).

ACK N OWLED G M ENTS
We thank the MRI PAS technicians, students, and volunteers for help in the field work and Jakub Bubnicki for the assistance with GIS data collection. The project was funded by the Ministry of Science and Higher Education (grant nos: N304 0522 33 and 2 P06L 006) and the budget of the Mammal Research Institute PAS.

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

AUTH O R S' CO NTR I B UTI O N S
BJ, MN, and SDC conceived the ideas and applied for the financial support, SDC and MN collected the samples, SDC, BJ, and TB analyzed the data, and SDC, MN, and BJ wrote the paper.

DATA ACCE SS I B I LIT Y
Sampling locations, ecological and climatic data and microsatellite genotypes supporting the results are deposited in Dryad repository: https://doi.org/10.5061/dryad.f8q8qf5.