INTRODUCTION
A significant role in the distribution of species is attributed to historical paleoclimate changes and the associated periodic existence of appropriate living conditions, tectonic movements, and processes of terrain evolution (e.g., Hrbek et al., Reference Hrbek, Stölting, Bardakci, Küçük and Wildekamp2004; Lister and Stuart, Reference Lister and Stuart2008; Markova and Kolfschoten, Reference Markova and van Kolfschoten2008; Fløjgaard et al., Reference Fløjgaard, Normand, Skov and Svenning2009; Conroy et al., Reference Conroy, Patton, Lim, Phuong, Parmenter and Höhna2016; Gottscho, Reference Gottscho2016; Baca et al., Reference Baca, Nadachowski, Lipecki, Mackiewicz, Marciszak, Popović, Socha, Stefaniak and Wojtal2017; Stojak et al., Reference Stojak, Borowik, Górny, McDevitt and Wójcik2019). Geologic processes may affect ecological (climatic, soil, and habitat) conditions of different areas by shaping their contemporary topography and lithology. This favors the formation or disappearance of habitats, preferring (or not) the occurrence of species with specific requirements. Such interconnections have been analyzed in recent years, and molecular analyses of genetic variability in populations and the genetic diversity of species have become tools for assessing landscape differentiation. Owing to such analyses, it is possible to identify even indistinct, “hidden” diversity within a population of the same species (Hrbek et al., Reference Hrbek, Stölting, Bardakci, Küçük and Wildekamp2004).
Today, such relationships have been documented for numerous plant and animal species in areas of strong tectonic and erosional activity, including North America (along the San Andreas Fault), Central America, Tenerife, Turkey, and New Zealand (Hrbek et al., Reference Hrbek, Küçük, Frickey, Stölting, Wildekamp and Meyer2002; Daza et al., Reference Daza, Castoe and Parkinson2010; Gutiérrez-García and Vázquez-Domínguez, Reference Gutiérrez-García and Vázquez-Domínguez2013; Waters et al., Reference Waters, Wallis, Burridge and Craw2015; Gottscho, Reference Gottscho2016; Puppo et al., Reference Puppo, Curto and Meimberg2016). Molecular methods indicate that highly diverse relief forms—ranging from river valleys to mountain chains, which developed because of tectonic movements or erosional processes, such as river captures or karst processes—are barriers for the migration of individuals of species living in those areas (Waters et al., Reference Waters, Craw, Youngson and Wallis2001; Hrbek et al., Reference Hrbek, Stölting, Bardakci, Küçük and Wildekamp2004, Reference Hrbek, da Silva, Dutra, Gravena, Martin and Farias2014; Burridge et al., Reference Burridge, Craw, Fletcher and Waters2008; Craw et al., Reference Craw, Burridge, Morris and Waters2008; White et al., Reference White, Davidson and Paquin2009). However, less work has focused on the relationship of genetic diversity and genetic differentiation of populations with geomorphological processes and changes occurring in the Pleistocene in repeatedly glaciated lowland areas.
The greatest role in shaping the distribution of contemporary species and their genetic diversity in northern and central Europe and Asia in the Pleistocene is generally attributed to the last glaciation (e.g., Jaarola et al., Reference Jaarola, Tegelström and Fredga1999; Brunhoff et al., Reference Brunhoff, Galbreath, Fedorov, Cook and Jaarola2003, Reference Brunhoff, Yoccoz, Ims and Jaarola2006; Kotlík et al., Reference Kotlík, Deffontaine, Mascheretti, Zima, Michaux and Searle2006; Wójcik et al., Reference Wójcik, Kawałko, Marková, Searle and Kotlík2010; Haring et al., Reference Haring, Sheremetyeva and Kryukov2011; Herman et al., Reference Herman, McDevitt, Kawałko, Jaarola, Wójcik and Searle2014; García-Vázquez et al., Reference García-Vázquez, Bilton, Foster and Ribera2017). Postglacial colonization by species and their genetic variability are considered to have been primarily controlled by climatic processes that took place in the Pleistocene and Holocene. These studies focused on several related species of voles, including the root vole (Microtus oeconomus) (Brunhoff et al., Reference Brunhoff, Galbreath, Fedorov, Cook and Jaarola2003, Reference Brunhoff, Yoccoz, Ims and Jaarola2006), the bank vole (Myodes glareolus) (e.g., Wójcik et al., Reference Wójcik, Kawałko, Marková, Searle and Kotlík2010), the common vole (Microtus arvalis) (Ratkiewicz and Borkowska, Reference Ratkiewicz and Borkowska2006; Stojak et al., Reference Stojak, McDevitt, Herman, Searle and Wójcik2015), and the field vole (Microtus agrestis) (e.g., Jaarola and Searle, Reference Jaarola and Searle2002). In this context, however, the impact of previous glacial episodes on the development of various habitat types should also be considered.
In the middle and late Quaternary, successive glaciations in the Northern Hemisphere covered various and generally successively smaller areas of Eurasia (Różycki, Reference Różycki1972; Baltrūnas et al., Reference Baltrūnas, Švedas and Pukelytė2007; Astakhov, Reference Astakhov2013). The resulting landscape shows zonation, well marked in Germany, Poland, Belarus, Lithuania, Latvia, Estonia, and part of Russia (Lindner, Reference Lindner1987; Marks, Reference Marks, Ehlers, Gibbard and Hughes2011; Astakhov, Reference Astakhov2013).
In central and eastern Europe, a good example of postglacial landscape diversity is eastern Poland (Fig. 1), where different types of postglacial landscape are present over a relatively short distance of several hundred kilometers with a latitudinal arrangement in accordance with the ranges of successive Scandinavian ice sheets (Różycki, Reference Różycki1972; Marks et al., Reference Marks, Dzierżek, Janiszewski, Kaczorowski, Lindner, Majecka, Makos, Szymanek, Tołoczko-Pasek and Woronko2016) (Fig. 1). Thus, the landscape maturity level, understood as transformation and even obliteration of glacial landscape elements, increases from north to south. This maturity was also influenced by the advance of denudation processes.
Location of the research sites in eastern Poland on the background of Pleistocene glaciations limits and the characteristics of surface relief shown on a digital elevation model. Limits of the main Pleistocene glaciations after Marks (Reference Marks2005, Reference Marks, Ehlers, Gibbard and Hughes2011). Southern range limit of the root vole in Poland after Sałata-Piłacińska (Reference Sałata-Piłacińska1990). B, Biebrza River; Nw, Narew River; Nk, Narewka River; Nc, Nurzec River; T, Tanew River; S, San River. *MIS 6 according to Marks et al. (Reference Marks, Dzierżek, Janiszewski, Kaczorowski, Lindner, Majecka, Makos, Szymanek, Tołoczko-Pasek and Woronko2016)
Considering the glaciation history and landscape diversity, related to the modern distribution of many plant and animal species, it can be assumed that in Europe and Asia: (1) glaciation history that varies across latitudes contributes to (a) different numbers and characters (e.g., different degrees of connectivity) of favorable habitats and (b) different colonization times of the postglacial landscape by species because of the presence of small local refugia in the ice-marginal zone; (2) different numbers of favorable habitats and different degrees of their communication/isolation results in different courses of in situ evolution, which should be reflected in the genetic structure; and (3) different connectivity of favorable habitats, their variability in attractiveness, and different availabilities contribute to different directions and intensities of contemporary species migration.
In postglacial areas, moisture-loving species are good for testing the impact of geologic history on their genetic differentiation because of specific habitat requirements. Postglacial humid (aqueous and boggy) environments undergo relatively fast evolution in a temperate climate, which in turn forces quick reaction of the inhabiting species to the changing environment.
The aim of this research was to determine whether and how the landscape shaped during subsequent Pleistocene glaciations affects the diversity and persistence of populations, expressed by the genetic differentiation of species. We focused our research on the root vole (tundra vole) Microtus oeconomus (Arvicolinae, Rodentia). We know that (1) this model species shows clear preferences for inhabiting wetlands (e.g., Tast, Reference Tast1966; Gliwicz and Jancewicz, Reference Gliwicz, Jancewicz, Jędrzejewska and Wójcik2004), (2) the southern limit of its contemporary natural range runs across southern Poland (Sałata-Piłacińska, Reference Sałata-Piłacińska1990; van Appeldorn, Reference van Apeldoorn, Mitchell-Jones, Amori, Bogdanowicz, Kryštufek, Reijnders, Spitzenberger, Stubbe, Thissen, Vohralik and Zima1999), and (3) the history of postglacial colonization of the species in this area is complex (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015).
GEOMORPHOLOGICAL CHARACTERISTICS OF THE STUDY AREA
The geologic history of eastern Poland has resulted in differences in the features and character of the land relief. From north to south, the following features clearly change: lake density, acreage and percentage of wetlands, the character of rivers, land surface gradients, and relief variations. Postglacial landscape zones showing different geomorphological and landscape features can be distinguished, consistent with the limits of the individual glaciations: Vistulian zone (V-Z), Wartanian zone (Wa-Z), Odranian zone (O-Z), and Sanian 1 and Sanian 2 zones (herein referred to as Sanian zone S-Z) (Marks, Reference Marks2005) (Fig. 1, Table 1).
Geomorphological characteristics of the landscape zones analyzed. Data sources: aMojski (Reference Mojski2005), bHarasimiuk et al. (Reference Harasimiuk, Dobrowolski, Rodzik and Radwan2002), cbased on various authors, after Lindner (Reference Lindner1987) and Mojski (Reference Mojski2005), dDembek et al. (Reference Dembek, Piórkowski and Rycharski2000), eChoiński (Reference Choiński2007), fMajdanowski (Reference Majdanowski1954). Drawn by T. Falkowski.
Vistulian/Weichselian zone (V-Z)
This is a young glacial area of the East Baltic Lake District (Kondracki, Reference Kondracki2011). The “fresh landscape” (Marks, Reference Marks2005) of this area is characterized by highly varied terrain relief, complicated geologic structure, and the presence of numerous lakes and boglands (Table 1). The surface relief of this area was ultimately formed during the Vistulian glaciation (115–11.7 ka BP; Marine Oxygen Isotope Stage [MIS] 2–5d) (e.g., Słowański, Reference Słowański1981; Szumański and Laskowski, Reference Szumański and Laskowski1993; Lisicki and Rychel, Reference Lisicki and Rychel2006; Marks et al., Reference Marks, Dzierżek, Janiszewski, Kaczorowski, Lindner, Majecka, Makos, Szymanek, Tołoczko-Pasek and Woronko2016). Young glacial relief landforms, such as strings of frontal moraine hills and consecutive glaciofluvial accumulation plains (e.g., Ber, Reference Ber2006), document the consecutive deglaciation stages in the individual phases of this glaciation.
A characteristic feature of this zone is the occurrence of extensive and numerous small melt-out depressions. Their surface is occupied by lakes and peat plains that developed from lakes being overgrown. In this zone, there are also elongated subglacial tunnel valleys or meltwater valleys, narrow and deeply incised in the upland areas. Nowadays, there are lakes or wetlands in these tunnel valleys (Szumański and Laskowski, Reference Szumański and Laskowski1993; Lisicki and Rychel, Reference Lisicki and Rychel2006). They are bounded by steep slopes attaining a height of 20–70 m. Such landforms often exist as isolated elements of the landscape. The relief diversity of the East Baltic Lake District, which results from glacial morphogenesis, and the poor degree of denudation processes are the reasons for the presence of many habitats preferred by the root vole, which are potential sites of occurrence of this rodent species. These are mainly lakeshores (still not overgrown by vegetation) or partially boggy and peaty melt-out depressions of various sizes. Many of these depressions represent isolated landscape patches. This area is subject to rapid evolution: lakes are being overgrown and turned into peat bogs. A measure of these changes is the decrease in lake density from 1.1% in 1954 to 0.9% in 1991 (Dembek et al., Reference Dembek, Piórkowski and Rycharski2000; Choiński, Reference Choiński2007).
Wartanian/Warthe zone (Wa-Z)
This area developed mainly during the Wartanian glaciation and then was partly denuded during the Eemian interglacial and the Vistulian glaciation. It comprises the Podlasie-Belarus Uplands and the eastern part of the Middle Polish Lowland (Kondracki, Reference Kondracki2011). In some publications, the Warta glaciation is termed the Odra glaciation stade (Marks et al., Reference Marks, Dzierżek, Janiszewski, Kaczorowski, Lindner, Majecka, Makos, Szymanek, Tołoczko-Pasek and Woronko2016). Despite the differences in the stratigraphic interpretation, this area is distinguished by a different landscape relief than other zones discussed here.
The area is characterized by a higher level of landscape maturity in relation to the northern zone: relief variations are smaller, and the geomorphological boundaries are gentler and less easily visible (Table 1, Wa-Z). Denudation processes are more advanced here than in the Vistulian zone. For this reason, the level of communication between humid landscape elements is higher. This area has developed natural habitats preferred by the vole. These are mainly peaty melt-out depressions, which were originally lakes, then included in the runoff system of surface waters at later stages of landscape development, and used by rivers as polygenetic sections of their valleys (Falkowska, Reference Falkowska2001; Falkowski, Reference Falkowski2003). For this reason, the valleys of the Biebrza (B), Narew (Nw), Narewka (Nk), and Nurzec (Nc) Rivers form a network of connected wetlands (Fig. 1, Table 1, Wa-Z). Melt-out landforms that are not included in the surface runoff system are now isolated boggy and peaty depressions. A feature that distinguishes this part of the area of Poland analyzed from the V-Z zone is the lack of lakes (Choiński, Reference Choiński2007) (Table 1).
In the area adjacent to melt-out depressions and river valleys of the glaciogenically modified glacial upland, there are also small depressions after minor dead-ice blocks. These are often filled with deluvial material deposited here mainly in periglacial conditions. Currently, peat deposition takes place within them. Such landforms are usually the routes of surface water runoff to rivers that appear periodically during precipitation events (Falkowska, Reference Falkowska2009).
The relief of the area is enriched in landforms whose development was related to the mode of deglaciation. These are kames and eskers occurring within the melt-out areas and uplands, kame terraces on the upland slopes, extensive levels of glaciofluvial plains, and ablation covers, as well as strongly denuded moraine hills and dunes formed in periglacial conditions (Bałuk, Reference Bałuk1991; Musiał, Reference Musiał1992; Falkowska, Reference Falkowska2009; Marks et al., Reference Marks, Makos, Szymanek, Woronko, Dzierżek and Majecka2019).
Odranian/Drenthe zone (O-Z)
The relief of the Polish part of Polesie has been shaped mainly by glacial processes (deposition and glaciotectonics) during the Odranian glaciation (Fig. 1, O-Z) (Lindner and Marks, Reference Lindner and Marks1994; Marks, Reference Marks2005) and subsequent periglacial processes, as well as denudation processes associated with warming periods of various stratigraphic ranks. For this reason, Polesie shows a higher level of landscape maturity in relation to both the northern zone and the boggy river valleys in the central zone (Fig. 1, Table 1) of the study area. The land surface relief is characterized by low diversity. Generally, this is a flat lowland area with relief variations reaching only a few meters and average gradients smaller than 1° (occasionally up to 3°). Larger hills are formed by outcrops of both older basement deposits and glacial forms: kames, eskers, moraine uplands, and single moraine hills. The most characteristic feature of this area is the presence of vast peat plains and bogs (as in Wa-Z) and numerous lakes (unlike in Wa-Z). The formation of these landforms was associated with glacial and thermokarst processes (Wojtanowicz, Reference Wojtanowicz1994) that took place in conditions of the shallow occurrence of Cretaceous poorly permeable marls and limestones under the glacial deposits. The watershed zones between these depressions form strings of eskers and kames. Large areas in this zone are occupied by flat, sandy, glacial (outwash) plains (Harasimiuk et al., Reference Harasimiuk, Dobrowolski, Rodzik and Radwan2002).
South Polish glaciations zone, Sanian1/Elsterian (Elstera1), Sanian2/Elsterian (Elstera2) zones (S1-Z, S2-Z; herein referred to as Sanian zone S-Z)
The South Polish glaciation zone of the study area is located within the Sandomierz Valley (Kondracki, Reference Kondracki2011) and is characterized by the highest level of landscape maturity (Fig. 1, Table 1, S-Z). The landscape shows signs of denudation here. The major features of the relief have developed because of the activity of the Sanian 1 and Sanian 2 ice sheets (Kwapisz, Reference Kwapisz1998; Wieczorek, Reference Wieczorek2006). The surface relief was also shaped by processes occurring in the periglacial zone during the later glaciations and by river erosion and accumulation processes that operated from the Mazovian interglacial to the present day (Starkel, Reference Starkel1988; Gębica, Reference Gębica2004). During this long period, denudation of the upland areas was accompanied by the filling of melt-out depressions and shaping of valleys. In many cases, signs of valley polygenesis, easily visible in the northern and central parts of the entire study area, have been completely obliterated in the terrain morphology here. Deposition of the largest loess succession and calm fluvial sedimentation, which led to the deposition of sandy alluvial sediments of considerable thickness, occurred in this area during the Vistulian glaciation (Kwapisz, Reference Kwapisz1998). Hydrogenic habitats, suitable for the root vole (Tast, Reference Tast1966), are associated within the S-Z zone mainly with river valleys (these occur predominantly at oxbow lakes or in near-slope boggy areas) (Falkowska and Falkowski, Reference Falkowska and Falkowski2015; Falkowska et al., Reference Falkowska, Falkowski, Tatur and Kałmykow-Piwińska2016) and with deflation zones on Pleistocene terrace surfaces overlain by eolian sands. Their formation might have been related to changes in climate humidity and a regional groundwater-level rise in the Holocene (e.g., in the Tanew River valley [T], Fig. 1) (Kwapisz, Reference Kwapisz1998). The moraine upland that occurs in the southern zone is characterized by a poorly diverse relief and the presence of various Pleistocene glaciofluvial landforms and dunes. Its surface zone is composed of glacial deposits and loess (Kwapisz, Reference Kwapisz1998). There are a few boggy areas (Table 1, S-Z).
MATERIALS AND METHODS
Animal trapping
In the period of 2008–2010, intense trapping of root voles was carried out at 33 sites in eastern Poland in four zones located within the limits of the individual Pleistocene glaciations (Fig. 1). Environmental conditions favorable for the occurrence of this rodent species were selected as research sites, including Carex wetlands, local boggy areas, lakeshores, and meadows in river valleys. The number of selected sampling sites in each zone was adequate to the proportion of habitats appropriate for the vole. The geographic distance between the trapping sites was 3.5–405 km.
Animal trapping was conducted from June to September or October—that is, in the growing season and breeding season, when the density of herbivorous rodents (including root voles) is the highest (Gliwicz and Jancewicz, Reference Gliwicz, Jancewicz, Jędrzejewska and Wójcik2004, Reference Gliwicz and Jancewicz2016). To catch animals, wooden single-capture live traps for rodents (8 × 9.5 × 16.5 cm) with bait and food (oat grains and pieces of juicy cucumber) were used. The traps were spaced every 10 m in lines (Jancewicz and Gliwicz, Reference Jancewicz and Gliwicz2017) at the most likely sites of vole occurrence. They were checked once a day. At a given site (20 traps), 7–23 individuals of M. oeconomus were captured. From each newly captured M. oeconomus individual, tissues were collected for examination (from the ear, as recommended by the ethics committee), and the animal was released at the point of capture. Recaptured individuals and those of other species were released at the trapping site immediately. The tissue samples were preserved in 70% alcohol and kept at –20°C until the analysis started. A total of 549 samples were collected from various M. oeconomus individuals from 33 sites.
All procedures of capture of root voles and sampling tissues were in accordance with Polish law and were accepted by the Third Local Ethical Commission in Warsaw (permission no. 45/2007). Trapping site locations were approved by the administrators of the area. All animals were released at the point of capture.
Laboratory analysis and analysis of genetic data
DNA was extracted from tissue samples of 549 individuals of M. oeconomus with the Genomic Mini Kit (A&A Biotechnology, Gdynia, Poland). For all individuals, 12 microsatellite loci were established according to the method described by Czajkowska et al. (Reference Czajkowska, Borkowska, Wieczorek and Zub2010). FSTAT (Goudet, Reference Goudet1995) was used to estimate the number of alleles per locus and allelic size range. Observed (Ho) heterozygosity was calculated in ARLEQUIN 3.11 (Excoffier et al., Reference Excoffier, Laval and Schneider2005). Departures from Hardy-Weinberg equilibrium (HWE) were calculated using Markov chain Monte Carlo with 1 million steps and 100,000 dememorization steps (ARLEQUIN 3.11).
STRUCTURE 2.3.4 was used to infer population structure and assign individuals to subpopulations (clusters) (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000; Falush et al., Reference Falush, Stephens and Pritchard2003). The number of genetic groups (K) was tested from 1 to 10. The results were summarized using the STRUCTURE HARVESTER v. 06.94 (Earl and Vonholdt, Reference Earl and Vonholdt2012). The same program was used to calculate the posterior probability of the data [LnP(D)] and delta K (ΔK) (Evanno et al., Reference Evanno, Regnaut and Goudet2005).
Mitochondrial DNA cytochrome b (mtDNA cyt b), based on the 908 bp fragment of the gene, was established for 427 individuals using the method described by Jancewicz et al. (Reference Jancewicz, Falkowska and Ratkiewicz2015). The relationships between the mtDNA haplotypes (GenBank accession number KP684101-KP684121; Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015) were illustrated by the minimum spanning network using ARLEQUIN 3.11.
ARLEQUIN 3.11 was also used to estimate the number of haplotypes (Nh), haplotype diversity (h), and nucleotide diversity (π). Genetic differentiation between populations was assessed as pairwise Fst for both markers, and statistically tested using this program (1000 permutations).
In GENALEX v. 6.0 (Peakall and Smouse, Reference Peakall and Smouse2006), principal component analysis (PCA) was performed for microsatellite DNA based on Fst and mtDNA Fst data.
All genetic parameters were averaged for four the geomorphological zones analyzed. Genetic parameters/genetic data/genetic variability of M. oeconomus—the number of alleles per locus (Na) and allelic size range, expected heterozygosity (He), number of haplotypes (Nh), haplotype diversity (h) and nucleotide diversity (π), genetic differentiation (Fst), and minimum spanning network of cytochrome b gene sequences—from eastern Poland were analyzed taking into account the zonal variability of the landscape, especially the zonal variability of the distribution and character of wetland habitats.
Statistical relations between geographic coordinates of the population and both the haplotype diversity (h) and nucleotide diversity (π) in these populations were determined by means of a regression curve in STATISTICA 10 (StatSoft Inc., 2012).
RESULTS
Genetic variability within geomorphological zones
Microsatellite DNA
Average allelic richness for different landscape zones and average allelic size range and average heterogeneity were determined. The lowest values for all parameters calculated were found in the Sanian zone (S-Z), and the highest values in the Wartanian zone (Wa-Z) (Table 2). The average values of the parameters calculated for a single population in the Wartanian zone (Wa-Z) and the Vistulian zone (V-Z) were always significantly higher, and the average values for the Sanian zone (S-Z) were always significantly lower than the average for the whole study area. In the case of the Odranian zone (O-Z), only the average allelic size range was higher than the average value for the whole study area. The remaining parameters did not exceed the average values for the whole study area. All populations were in HWE.
Genetic variability parameters for the root vole (Microtus oeconomus) within the geomorphological zones analyzed, genotyped at 12 microsatellite loci (division into zones and places of research on the Fig. 1). Because the populations were in Hardy-Weinberg equilibrium (HWE), only He values have been given. Fst, genetic differentiation value for zones; He, expected heterozygosity; Na, number of different alleles per locus; Nind, number of individuals analyzed. Fst according to the Wright scale: 0.000–0.05, little differentiation; 0.05–0.15, moderate differentiation; 0.15–0.25, great differentiation; >0.25, very great differentiation (Hartl and Clark, Reference Hartl and Clark1997).
According to Wright's scale (Hartl and Clark, Reference Hartl and Clark1997), the average pairwise Fst between the populations in the landscape zones for nuclear DNA was low in the Vistulian zone (V-Z) and moderate in the remaining landscape zones. The average Fst values increase across zones from north to south (Table 2).
Because of reaching a constant value by LnP(D) at K = 5, STRUCTURE analysis was carried out for divisions up to K = 5 (Fig. 2). ΔK showed the highest support for K = 3. For this value of K, STRUCTURE identified 3 subpopulations: (1) Vistulian zone, (2) Wartanian and Odranian zone, and (3) Sanian zone. The Odranian zone was distinctly separated only for K = 5. Division at K = 2 was also very clear and indicated clusters: the Sanian zone and the rest of the populations analyzed. In mtDNA, these two clusters correspond to the central European phylogroup (CE) and South Polish (CE-PL S) group described by Jancewicz et al. (Reference Jancewicz, Falkowska and Ratkiewicz2015) (Table 3). For K from 2 to 4, the Odranian zone is only slightly marked, and that shows that it may be a mixing zone for voles from the Sanian and Wartanian zones.
(color online) Results of STRUCTURE analysis of root vole data based on 12 microsatellite loci assuming K = 2 to 5 genetic clusters. Vertical lines separate voles from different geomorphological zones.
Genetic variability parameters of mitochondrial DNA cytochrome b of the root vole (Microtus oeconomus) within the geomorphological zones analyzed (see Fig. 1). CE, central European haplotype group; CE-PL S, South Polish group–distinct haplotype group within the CE phylogroup (see Fig. 4); Fst, genetic differentiation value for zones; h, haplotype diversity; Nh, number of haplotypes within the zones; nh, number of haplotypes at a single study site; Nind, number of individuals analyzed; π, nucleotide diversity. Scale for Fst as in Table 2.
mtDNA cytochrome b
The haplotype diversity h and nucleotide diversity π were highly variable within the zones. The highest values of both parameters were recorded in the Wartanian zone (Wa-Z, Table 3, Fig. 3). In the Sanian zone (S-Z), determination of Fst was impossible because of the fixation of one haplotype (for more information, see Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015). We found correlations (P < 0.0001) between the latitude and the values of both parameters (h and π) here. A correlation between these values and the longitude, and an increase along the west-to-east axis (P < 0.001), were found only for the nucleotide diversity (Fig. 3, right side).
The relationship between haplotype diversity (h; A, C) and nucleotide diversity (π; B, D) for mitochondrial DNA of the root vole (Microtus oeconomus) and the geographic location of 33 research sites (λ, latitude [A, B]; φ, longitude [C, D]). O-Z, Odranian zone; S-Z, Sanian zone; V-Z, Vistulian zone; Wa-Z, Wartanian zone.
The average pairwise Fst between the populations within the landscape zones was high in the Wartanian zone (Wa-Z) and very high in the Vistulian zone (V-Z) and the Odranian zone (O-Z) (Table 3). The average Fst was greater than the average value for the whole study area only in the Odranian zone.
The glacial landscape zones differed in the number of recorded haplotypes: from 1 in the Sanian zone (S-Z) to 17 in the Wartanian zone (Wa-Z) (with an average of 1 to 4.2 for each locality, and an overall average of 2.9, as reduced to a single locality) (Table 3). In the Wartanian zone (Wa-Z) and Odranian zone (O-Z), both haplotypes were recorded, which belong to the central European clade CE (Table 3, Fig. 4). The haplotypes of these groups show different distribution patterns. In the most southerly Sanian zone (S-Z), there is only the dominant haplotype from the South Polish group CE-PL S. Haplotypes of this group are absent from the peripheral, most northerly Vistulian zone (V-Z). In the remaining landscape zones—the Wartanian zone (Wa-Z) and the Odranian zone (O-Z)—haplotypes from both phylogenetic groups are recorded (Table 3). Out of 21 haplotypes described (for details, see Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015), only 5 occur in at least two distinguished landscape zones: 4 haplotypes were found in two zones, and only 1 haplotype in three zones (Fig. 4). The distribution of the remaining 16 haplotypes, commonly of low and very low frequency, is limited to a single landscape zone. Widely distributed haplotypes were usually found in landscape zones neighboring each other. Only haplotype PL-9 (Fig. 4), dominant in the South Polish group CE-PL S, had a different distribution pattern: its presence was recorded in the Wartanian zone (Wa-Z) and the Sanian zone (S-Z), but it was absent from the intervening Odranian zone (O-Z).
Minimum spanning network of 21 cytochrome b gene sequences of the root vole (Microtus oeconomus) in eastern Poland. Phylogenetic groups (central European haplotype group [CE] and South Polish group–distinct haplotype group within the CE phylogroup [CE-PL S]) after Jancewicz et al. (Reference Jancewicz, Falkowska and Ratkiewicz2015). The size of the circles corresponds with the frequency of haplotypes in the entire sample of 427 individuals analyzed. The presence of haplotypes in the geomorphological zones analyzed (Fig. 1) is marked. O-Z, Odranian zone; S-Z, Sanian zone; V-Z, Vistulian zone; Wa-Z, Wartanian zone.).
Genetic differentiation between the geomorphological zones
The average pairwise Fst between the zones was very high (Wright scale) in all mtDNA pairwise comparisons (Table 4, above the diagonal). The highest values were between the Sanian zone (S-Z) and the remaining zones. For microsatellite DNA, the pairwise Fst was low between the neighboring Vistulian zone (V-Z) and the Wartanian zone (Wa-Z), and large between the Sanian zone (S-Z) and the Wartanian zone (Wa-Z) and between the Sanian zone (S-Z) and the Odranian zone (O-Z). In the remaining three comparisons, the pairwise Fst was moderate (Table 4, below the diagonal).
Pairwise genetic differentiation Fst (average ± standard deviation) for voles computed between geomorphological zones: below the diagonal, for microsatellite markers; above the diagonal, for mitochondrial DNA (mtDNA). Fst according to the Wright's scale as in Table 2.
PCA performed on the Fst values reflects the subdivision into the four postglacial landscape zones but with overlap between some populations, particularly those from the Wartanian and Odranian zones (Fig. 5). This subdivision is more clearly pronounced for the nuclear marker (Fig. 5, upper image), although one population from the Odranian zone (O-Z) shows similarity to the population from the Wartanian zone (Wa-Z). The PCA based on Fst values for mtDNA is slightly different: the zones are clearly distinct, but the populations of the Odranian zone are contained within the Wartanian zone. In addition, two populations from the Vistulian zone were grouped with populations from the neighboring Wartanian zone, and one population from the Wartanian zone was grouped within the Sanian zone (Fig. 5, lower graph).
Principal component analysis of the root vole (Microtus oeconomus) representing the geomorphological zones analyzed for Fst for microsatellite DNA (A) and mitochondrial DNA cytochrome b (mtDNA cyt b; B).
DISCUSSION
The linking of geologic processes with the occurrence of species, unequivocally postulated relatively recently by Conroy et al. (Reference Conroy, Patton, Lim, Phuong, Parmenter and Höhna2016) and Gottscho (Reference Gottscho2016), seems to be necessary to explain genetic variability in populations of contemporary animal species. The results presented in our work suggest that—apart from tectonic and erosional processes—Pleistocene glacial processes can also affect the genetic differentiation of a population. Most phylogeographic studies make general reference to the effects of glaciations and do not consider geomorphology in any detail (e.g., Brunhoff et al., Reference Brunhoff, Galbreath, Fedorov, Cook and Jaarola2003; Haring et al., Reference Haring, Sheremetyeva and Kryukov2011; Hulejová Sládkovičová et al., Reference Hulejová Sládkovičová, Dąbrowski, Žak, Miklós, Gubányi, La Haye and Bekker2018). The results of our research indicate a relationship between the genetic structure of M. oeconomus and the diversity of landscape relief, geologic structures, and geologic history of the glacial area. Genetic differences between populations inhabiting different glacial zones are visible despite the fact that the voles studied are included by mtDNA into a single phylogenetic group—the central European phylogroup CE (Brunhoff et al., Reference Brunhoff, Galbreath, Fedorov, Cook and Jaarola2003; Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015), in which the voles from the southern zone differentiate as a separate branch, CE-PL S, that survived the last glaciation in the local northern latitude refugium situated in Poland (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015). The separate origin of voles from the southern part of Poland living close to the range of occurrence of the species is also supported by results of the microsatellite analysis that include determination of basic genetic variability parameters: number of alleles per locus, allelic size range, expected heterozygosity (Table 2), PCA grouping analysis (Fig. 5, upper figure), and clustering for two clusters (for K = 2) in STRUCTURE (Fig. 2).
There are many studies showing the influence of contemporary environmental and landscape fragmentation and the impact of various types (geographic, topographic, geologic, environmental, and climatic) of barriers on the genetic structure of plant and animal populations (e.g., Gerlach and Musolf, Reference Gerlach and Musolf2000; Ratkiewicz and Borkowska, Reference Ratkiewicz and Borkowska2006; Banaszek et al., Reference Banaszek, Jadwiszczak and Ziomek2011; Huck et al., Reference Huck, Budel and Schmitt2012; Smissen et al., Reference Smissen, Melville, Sumner and Jessor2013; Fietz et al., Reference Fietz, Tomiuk, Loeschcke, Weis-Dootz and Segelbacher2014; Czarnomska et al., Reference Czarnomska, Niedziałkowska, Borowik and Jędrzejewska2018). However, there have been no studies so far that would focus on explaining the role of postglacial landscape differentiation on a larger temporal and spatial scale. The research presented here discusses this problem based on a very large sample (549 individuals) and a large amount of data collected from a morphologically diverse area, shaped over several Pleistocene glaciations and intervening warming periods.
The study area has been shaped by geologic processes occurring during these glaciations and by periglacial and interglacial processes (Lindner and Marks, Reference Lindner and Marks1994, Reference Lindner and Marks1999; Marks et al., Reference Marks, Makos, Szymanek, Woronko, Dzierżek and Majecka2019). The clear subdivision of the area of eastern Poland into four zones, consistent with the extents of consecutive glaciations (Fig. 1, Table 1), allowed testing the genetic-landscape relations in differently developed environments. The impact of modern climate, which is very important in biogeographic research, was less significant here, because the study area lies within one climate zone—the temperate climate zone. The diversity of climatic conditions is not high enough to affect the contemporary distribution of the boreal root vole. In terms of climatic conditions, an exception is the Polesie region (within the Odranian zone), in which a slight influence of continental climate is observed (Kaszewski, Reference Kaszewski and Radwan2002).
The PCA (Fig. 5) enabled grouping the populations based on Fst between the populations and divided them into groups largely according to the division of the landscape into zones. This division is very clearly visible in analyses of microsatellite DNA (Fig. 5, upper graph), which is a marker inherited from both parents and reflects the spatial distribution of the voles. However, this distribution is not so clearly visible in mtDNA, which is inherited only in the maternal line (Fig. 4, lower graph). In PCA grouping of mtDNA, the lack of a clear separation of populations of the Odranian zone from the Wartanian zone may result from the fact that the landscapes of these zones could have been shaped in a similar period and under similar conditions. The stratigraphic position of these zones is not unambiguous because for some authors (Marks et al., Reference Marks, Dzierżek, Janiszewski, Kaczorowski, Lindner, Majecka, Makos, Szymanek, Tołoczko-Pasek and Woronko2016) the Wartanian glaciation is considered to be the stade of the Odranian glaciation (MIS 6). This pattern and the way of mtDNA grouping in PCA can be partly explained by the haplotype network and by the presence or absence of these haplotypes in populations within the zones analyzed (Fig. 4). Despite such a result of haplo-grouping of mtDNA in PCA, the pairwise comparison of the zones shows very large differences in Fst (Table 4), which suggests limited flow of haplotypes between the zones. This is also visible in the haplotypes network (Fig. 4), in which the occurrence of haplotypes in different landscape zones is additionally marked.
STRUCTURE (Bayesian clustering) shows very similar results. In this analysis, the highest fit (probability) was found for division of the populations into three clusters in accordance with the Vistulian, Wartanian and Odranian, and Sanian zones (Fig. 2). However, populations from the Odranian zone form a separate cluster for K = 5 (Fig. 2). This could be caused by the presence of a mixed zone there. In addition, the distinctness of the Odranian zone is visible in parameters from Table 2. The number of alleles per locus and allelic size range for voles from this zone reveal a very wide range of values (minimum to maximum): from low values characteristic for the Sanian zone to very high ones that were found in the Vistulian zone.
As we expected, the vole finds its best conditions for living in the area of vast boggy valleys of the Biebrza, Narew, Narewka, and Nurzec Rivers (Fig. 1) in the Wartanian zone. High genetic diversity and usually low Fst values between populations in this area indicate free mating of large numbers of individuals and free flow of genes within a zone. The Wartanian zone landscape (Wa-Z, Table 1) began to be shaped earlier than the Vistulian zone landscape (V-Z, Table 1); it developed during the Odranian and Wartanian glaciations and then was modified by denudation processes in the Eemian interglacial and the Vistula glaciation. These processes resulted in the “levelling” of the terrain relief. The lakes that remained in melt-out depressions were later partly inundated with sediments from interglacial rivers and then overgrown, forming extensive peat bogs. Most of these depressions have later become river valleys (Falkowski, Reference Falkowski1975; Falkowska, Reference Falkowska2001, Reference Falkowska2009). On the other hand, the hills were largely denuded, which has limited the number of natural barriers for the migration of individuals of this small rodent. The landscape of this zone favors the presence of moisture-loving species, a free spread of individuals, stability of spatial relations between them, and free mixing of genes. In the case of the root vole, this is indicated mainly by the highest number of alleles per locus, the range of allelic size range, and high heterozygosity (Table 2).
The nature and distribution of habitats (including extensive boggy river valleys) in the landscape of this zone translate into the continuity of habitats optimal for voles and thus promote the free spread of individuals, stability of spatial relations between them, and free mixing of genes between populations. The lack of permanent isolation in the landscape is the result of its maturity level. This type of environment favors free migration of the root vole, although this species has very low dispersal capability (Steen, Reference Steen1994). Relations between haplotypes, visible in the haplotype network, and the distribution of haplotypes closely related to each other within the zones (Fig. 4) reflect the sedentariness of females, characteristic of this species (e.g., Gliwicz and Jancewicz, Reference Gliwicz, Jancewicz, Jędrzejewska and Wójcik2004; Pilot el al., Reference Pilot, Dąbrowski, Jancewicz, Schtickzelle and Gliwicz2010). In stable and heterogeneous environments, which are suitable for this rodent, daughters settle near mothers and other closely related females, which causes an increase in the number of haplotypes (because of fixed mutations and a relatively large number of individuals) in a small space. This is more possible in the extensive favorable environments of the Wartanian zone (Wa-Z), resulting in an increase in the gene diversity. This population parameter, seemingly insignificant to landscape researchers, is important for understanding both the mechanism of functioning of the population and the effects of environmental changes for the population in the past and in the future.
Because of this, genetic diversity in stable and rich environments of the Wartanian zone (Wa-Z), expressed in the diversity of haplotypes, nucleotides, and alleles, is high (Tables 2 and 3). It is higher than in the other zones where other types of landscape predominate, and the environments have different features. It can be assumed that this diversity will continue to be high, even if there were periodic decreases in the population size, or rapid reductions, in the past (Pilot et al., Reference Pilot, Dąbrowski, Jancewicz, Schtickzelle and Gliwicz2010), because only rich environments guarantee the appropriate level of genetic differentiation. Similarly high allelic diversity was observed in other single populations of root vole within the Wartanian zone—in the Narewka valley (Pilot et al., Reference Pilot, Dąbrowski, Jancewicz, Schtickzelle and Gliwicz2010; Hulejová Sládkovičová et al., Reference Hulejová Sládkovičová, Dąbrowski, Žak, Miklós, Gubányi, La Haye and Bekker2018) and the Biebrza valley (Czajkowska et al., Reference Czajkowska, Borkowska, Wieczorek and Zub2010).
Because of a different geologic history of the Odranian zone and the impermeable basement (Cretaceous marls) to the glacial deposits, bogs and lakes occupy a large part of the area (O-Z, Table 1). High lake density makes this area similar to the northern Vistulian zone (V-Z), but because of the extensive areas of bogs and much lower relief gradients than in the north, this area is similar to the Wartanian zone (Wa-Z) in terms of landscape features. The Odranian zone (O-Z) was shaped in a specific way (thermokarst process) during the Odranian glaciation, on the basement composed of marls (Harasimiuk et al., Reference Harasimiuk, Dobrowolski, Rodzik and Radwan2002). This difference can also be seen in the mtDNA genetic structure of the local root vole population. No PL-9 haplotype from the separate South Polish branch of the central European clade, CE-PL S, was found here (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015), which is dominant in eastern Poland (Fig. 4) and fixed in the immediately neighboring Sanian zone (S-Z). However, the presence of two other haplotypes, which have not been found in the other zones, was recorded only in the O-Z: haplotype PL-7, directly originating from PL-9, and PL-8, a rare (present in this area only) haplotype from the central European group (Fig. 4). Perhaps the PL-9 haplotype, although most numerous among the 21 haplotypes found in this area, has not been revealed during this study, yet still exists here. A similar situation was previously observed in the root vole population studied for many years in eastern Poland (Dąbrowski et al., Reference Dąbrowski, Pomorski and Gliwicz2013). In spite of some geologic differences of the Polesie region, spatial analysis of molecular variance (SAMOVA) of mtDNA includes the root vole population from this area to the common group of populations from the Vistulian, Wartanian, and Odranian zones (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015). To some extent, this is also suggested by the mtDNA PCA analysis (Fig. 5, lower graph). However, based on analyses of microsatellite DNA transmitted to progeny by mothers and fathers, we can distinguish the vole populations from the Odranian zone as a separate population group (Fig. 5, upper graph). This population was also distinguished in STRUCTURE analysis for K = 5 (Fig. 2).
In the most northerly and southerly zones (i.e., in the Vistulian zone [V-Z] and the Sanian zone [S-Z]), moist habitats are scattered and the fragmentation is more pronounced. However, the lack of habitat continuity differs in character between both these zones. In the oldest landscape of the southern zone, the habitat discontinuity is ordered and linear (S-Z, Fig. 1, Table 1). Moist habitats are concentrated along river valleys: in oxbow lakes and near-scarp boggy areas, and in deflation depressions on upper river terraces overlain by eolian sands. In the northerly located Vistulian zone, the terrain relief gives an image of random and chaotic fragmentation, which results from the young age of the landscape and the short period of its evolution. Here, in topographic lows, moist places are more frequent, but their mutual arrangement seems to be relatively disordered (V-Z, Fig. 1, Table 1). However, a system of connections exists between the moist terrain depressions that can be convenient migration routes despite considerable relief variations and gradients. The evolution of this young glacial area, already observed nowadays, consists of denudation of uplands, shallowing of lakes and their transformation into bogs (a change of lake density by 0.2% over 37 yr; Dembek et al., Reference Dembek, Piórkowski and Rycharski2000; Choiński, Reference Choiński2007), and their incorporation into the surface outflow system. This will increase both the area occupied by hydrogenic habitats and the communication between them. It can be expected that these phenomena will increase the genetic diversity of the root vole, which has already been observed to some extent (six haplotypes and heterozygosity of 0.8021; Tables 2 and 3). If this is the direction of its evolution, it is possible that, in the future, the system of habitats in the Vistulian zone will be similar to the system observed currently in the Wartanian zone.
The present-day differences in the terrain relief and the system of habitats in the Vistulian zone (V-Z) and the Sanian zone (S-Z) give both peripheral zones different characters in terms of root vole requirements. The southern zone of the landscape studied (S-Z) is special also because it is currently inhabited by voles that probably originate from a small northern glacial refugium located in the territory of Poland (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015). This different genesis, previously described/postulated on the basis of analysis of mtDNA (inherited only in the maternal line), is supported by the results of the analysis of the other genetic marker provided herein—the nuclear (microsatellite) DNA. The allelic size range, which is clearly lower than in other zones (Table 2), indicates the different origin of voles currently inhabiting the zone located close to the limit of its present-day occurrence (Sałata-Piłacińska, Reference Sałata-Piłacińska1990; Linzey et al., Reference Linzey, Shar, Lkhagvasuren, Juškaitis, Sheftel, Meinig, Amori and Henttonen2016; see also Fig. 1). In addition, this thesis is supported by the high pairwise Fst between the S-Z and O-Z and between the S-Z and Wa-Z zones for microsatellite DNA (Table 4), clearly higher than between the S-Z and V-Z.
The voles studied from the south of Poland (Fig. 1, localities 28–33) do not show mtDNA genetic variability (Table 3). All individuals have descended from a common (pra)female. The loss of genetic variability by the fixation of one haplotype indicates a strong reduction in abundance in the past, which is undoubtedly because of the limited availability of suitable habitats and the lack of, or difficult, communication between them in the Pleistocene. This may be because of the origin of the local refugium (Jancewicz et al., Reference Jancewicz, Falkowska and Ratkiewicz2015). Paleozoological research provided evidence for the presence of rodents in southern Poland during the Vistulian. On the Kraków-Częstochowa Upland, in the Holy Cross Mountains and in the Carpathians, fossil remains of root voles, common voles, field voles, bank voles, and other mammals were found (e.g., Madeyska, Reference Madeyska1981; Nadachowski, Reference Nadachowski1989; Nadachowski et al., Reference Nadachowski, Harrison, Szyndlar, Tomek and Wolsan1993; Sommer and Nadachowski, Reference Sommer and Nadachowski2006; Żarski et al., Reference Żarski, Winter, Nadachowski, Urbanowski, Socha, Kenig and Marcinkowski2017). In small and isolated populations, rapid decreases in abundance are dangerous, because these populations quickly lose their genetic variability because of the fixation of alleles (Leijs et al., Reference Leijs, van Apeldoorn and Bijlsma1999; van de Zande, Reference van de Zande, van Apeldoorn, Blijdenstein, de Jong, van Delden and Bijlsma2000; Báldi et al., Reference Báldi, Csorba and Korsós2001). Heterozygosity in the Sanian zone (S-Z) is relatively low (0.687, Table 2) and is comparable to that found in the endangered M. oeconomus populations in isolated and fragmented environments of the Netherlands and Hungary (van de Zande et al., Reference van de Zande, van Apeldoorn, Blijdenstein, de Jong, van Delden and Bijlsma2000; Hulejová Sládkovičová et al., Reference Hulejová Sládkovičová, Dąbrowski, Žak, Miklós, Gubányi, La Haye and Bekker2018). The contemporary Sanian zone environment, although severe in terms of the existence of moisture-loving species, provides the possibility of migration and exchange of genes between the vole populations. This is reflected in the genetic differentiation of nuclear DNA inherited from both parents.
In the Sanian zone (S-Z), hygrophilous species have a limited number of places to settle. This “dilution” of the distribution, as we move from the center of the species range to its borders, is a universal phenomenon related to the decreasing number of habitable environments (Brown, Reference Brown1984). From the point of view of M. oeconomus, the mature river valleys in a highly denuded landscape do not provide enough places and space to live, and the limited number of potential places to colonize increases the isolation between populations. It also impedes the migration of individuals, thus restricting the free mixing of genes. This phenomenon overlaps with the universal model of distribution and functioning of species at the range margins (Hampe and Petit, Reference Hampe and Petit2005), because the range limit of M. oeconomus runs across this zone. It seems, therefore, that the lack of mtDNA variability in the Sanian zone (Table 3) and relatively low heterozygosity (Table 2) are a result of several factors—the lack of appropriate habitats and a strong reduction in abundance in the past, which may also result from the lack of habitats in a highly mature and transformed landscape. The result is probably the current course of the range limit of this boreal species. Paleozoological data show a completely different (wider) range of the root vole in Europe (Madeyska, Reference Madeyska1981; Nadachowski, Reference Nadachowski1989; Kowalski, Reference Kowalski2001; Falkowska, E., Jancewicz, E., unpublished manuscript), which was undoubtedly related to the distribution of favorable habitats.
CONCLUSIONS
The landscape of eastern Poland, shaped during consecutive Pleistocene glaciations, forms a specific system of hydrogenic habitats in a zonal pattern. This geographic patterning is reflected in the genetic structure of the moisture-loving root vole.
The genetic structure of the root vole shows that the best living conditions for the root vole persist in the extensive habitats of the Wartanian zone with a high degree of connectivity—in boggy river valley systems. Such conditions promote free migration, random mixing of genes, and increase of genetic diversity. The lowest genetic diversity indicates that the most severe conditions for the vole occur in the Sanian zone, which offers a completely different pattern of habitats preferred by the vole. Its original landscape developed during the oldest glaciations. Changes in both the nature and connectivity of environments, resulting from the geologic evolution of postglacial areas, may be the beginning of a change in species ranges. This is supported by the results of molecular research.
The relationships found as a result of the studies may become to some extent an indicator useful in explaining the geologic history of the postglacial region. This issue requires further analysis concerning other species with different habitat preferences. Thus, it is appropriate to combine geologic knowledge with biological data in order to understand fully the biogeographic history and contemporary distribution of a species in postglacial areas.
ACKNOWLEDGMENTS
This study was financed by research grant no. N304 232035 from the Polish National Science Centre NCN to E. Jancewicz, and the budget of the Faculty of Forestry of the Warsaw University of Life Sciences – SGGW and the Faculty of Geology of the University of Warsaw. The authors thank Mirosław Ratkiewicz and Magdalena Czajkowska from the University of Bialystok for conducting molecular analyses and Mirosław Ratkiewicz for valuable comments and discussion. The authors are grateful to Joanna Gliwicz from the Museum and Institute of Zoology Polish Academy of Sciences and Leszek Marks from University of Warsaw for the constructive comments with respect to the manuscript. The authors would like to express their gratitude to the reviewers and the associate editor for their time and effort dedicated to improving the manuscript for publication in Quaternary Research. The authors also thank Tomasz Falkowski for drawing diagrams of geomorphological zones.
