Landscape genetics of northern crested newt Triturus cristatus populations in a contrasting natural and human-impacted boreal forest

Among vertebrates, amphibians currently have the highest proportion of threatened species worldwide, mainly through loss of habitat, leading to increased population isolation. Smaller amphibian populations may lose more genetic diversity, and become more dependent on immigration for survival. Investigations of landscape factors and patterns mediating migration and population genetic differentiation are fundamental for knowledge-based conservation. The pond-breeding northern crested newt (Triturus cristatus) populations are decreasing throughout Europe, and are a conservation concern. Using microsatellites, we studied the genetic structure of the northern crested newt in a boreal forest ecosystem containing two contrasting landscapes, one subject to recent change and habitat loss by clear-cutting and roadbuilding, and one with little anthropogenic disturbance. Newts from 12 breeding ponds were analyzed for 13 microsatellites and 7 landscape and spatial variables. With a Maximum-likelihood population-effects model we investigated important landscape factors potentially explaining genetic patterns. Results indicate that intervening landscape factors between breeding ponds, explain the genetic differentiation in addition to an isolation-by-distance effect. Geographic distance, gravel roads, and south/south-west facing slopes reduced landscape permeability and increased genetic differentiation for these newts. The effect was opposite for streams, presumably being more favorable for newt dispersal. Populations within or bordering on old growth forest had a higher allelic richness than populations in managed forest outside these areas. Old growth forest areas may be important source habitats in the conservation of northern crested newt populations.


Introduction
Much of the widespread decline in global biodiversity (IPBES 2018) is driven by negative impacts from anthropogenic habitat changes and loss. Among all vertebrates, amphibians is the group that currently has the highest proportion of threatened species, mainly because of habitat loss, caused by agricultural land use, logging, and the changing of freshwater systems (Stuart et al. 2004;Becker et al. 2007;Baillie et al. 2010).
Habitat loss is a negative direct effect, as smaller habitat patches can sustain fewer individuals (Fahrig 2003;Cushman 2006;Fischer and Lindenmayer 2007), i.e., population reduction. However, the same process, typically also leads to habitat changes and split, i.e., an increase in the amount of inhospitable environments between populations, which can reduce landscape permeability and migration between suitable habitat patches (Wiegand et al. 2005;Becker et al. 2007; 1 3 Fischer and Lindenmayer 2007). Connectivity among habitats is suggested to play a key role in preserving amphibian populations (Lehtinen et al. 1999;Cushman 2006;Becker et al. 2007;Coster et al. 2015a).
The potential for migration to counteract habitat isolation effects is connected to the balance between a given species ability and its propensity to move and the geographical distance to travel [the isolation-by-distance effect (Wright 1943;Hutchison and Templeton 1999;Jenkins et al. 2010)]. However, it is also connected to the species ability to traverse the in-between (matrix) habitat (i.e., landscape resistance Jenkins et al. 2010;Richardson 2012;Balkenhol et al. 2016), depending on the species vagility, plus the ability to cross different types of environments. Amphibians are considered unlikely candidates for good dispersal, with their generally low vagility (Bowne and Bowers 2004;Alex Smith and Green 2005), small body sizes and high water loss rate under hot and dry conditions (Oke 1987;Wells 2007). Therefore, landscape factors determining matrix habitat likely may be important in determining the level of connectivity between habitats, and particularly for amphibians (Compton et al. 2007;Todd et al. 2009;Goldberg and Waits 2010).
The general concept and framework of landscape genetics combines population genetics with features from landscape ecology, offering tools to investigate relative effects of landscape composition traits and the putatively associated configuration upon gene flow and genetic drift (Balkenhol et al. 2016). Factors reported to increase genetic differentiation in an amphibian context, are roads (Richardson 2012;Sotiropoulos et al. 2013), rivers (Peter et al. 2009;Richardson 2012), topography (Spear and Storfer 2008;Kershenbaum et al. 2014), urban areas (Emaresi et al. 2011) and open fields (Greenwald et al. 2009), although effects can be species specific or landscape specific. Because most of these studies focus on landscapes affected by agriculture and urban development, an important observation is that forest cover appears to decrease population differentiation (Greenwald et al. 2009;Richardson 2012). However, although forests may still constitute extensive and important habitats for amphibians (Corn and Bury 1989;Gibbs 1998;Semlitsch and Bodie 2003;Cushman 2006;Coster et al. 2015b), few studies focus on contrasting forest ecosystems, which may be disturbed by human impacts, or remain natural.
Here, we focus on a boreal northern forest ecosystem without major human impacts like urban areas, agriculture or major roads. However, the study system consists of boreal forest with different histories of generally low, but quite variable and contrasting impacts from forestry practices. One part of the study area has been affected by clear-cutting and harvest for several decades, whereas another part of the study area consists of old natural Pinus sylvestris dominated forest with long forest continuity (trees 200-300 years old, rejuvenated by natural forest fires) and no record of recent logging (Reiso 2018).
We suggest that such forest habitat changes by clear-cutting likely is negative for amphibian genetic connectivity. However, we also hypothesize that additional landscape features, both natural and infrastructure associated with clearcutting, primarily roads, may play an important role in these generally less human-impacted landscapes. This requires more detailed data and analysis, but is important particularly in proactive conservation. Here, factors like aspect and vegetation cover may be important considering the distribution of cool and humid microclimates, and thereby the accessibility of suitable habitat for amphibians (Oke 2002;Peterman and Semlitsch 2013). Streams could function as humid dispersal corridors (Emel and Storfer 2015), steep slopes as barriers or partial barriers by evoking avoidance behavior or increasing energy cost (Lowe et al. 2006). Low soil productivity and the removal of forest canopy could result in a lack of prey, as invertebrate abundance may be affected negatively by clear-cuts and low soil pH (Stuen and Spidsø 1988;Wareborn 1992;Atlegrim and Sjöberg 1996). Further, forest gravel roads could act as barriers due to a drier microclimate from removal of vegetation and canopy (Marsh and Beckman 2004), or because of steep roadside verges that are too difficult to traverse (Marsh et al. 2005).
Here we use the northern newt as an amphibian model to study landscape genetic relationships and identify what landscape features may affect the genetic differentiation and pattern of northern crested newt populations in a northern boreal forest breeding pond system with contrasting (near-) natural vs. forestry impacted forest landscapes. Northern crested newt is a pond breeding amphibian. They can become at least up to 16 years in the wild (Miaud and Castanet 2011). A female can lay about 200 eggs annually (Arntzen and Hedlund 1990), but because of lethal homozygotes there is a 50% lethality during egg development (Wallace 1987). Sexual maturity is reached at age 3-5 years (Dolmen 1983). Known maximum dispersal distance of the species is 1 km or more (e.g., 860 m Kupfer and Kneitz 2000, 1290m Kupfer 1998. However, long distance dispersal events could be difficult to observe and generalize because of the rarity of such events, and the likely dependence on factors such as corridors and matrix between ponds.

Study area
Thirteen known, fishless northern crested newt breeding ponds were included for field sampling in the study area ( Fig. 1), located within a forested land area of 10.5 × 3.5 km (N59° 37′, E9° 19′) in south-east Norway. The sampled ponds are all in the southern part of a likely network of breeding ponds (based on estimated distances). No intervening breeding ponds were found between sampled ponds or to the south. Mean inter-pond distance is 3841 m (SD ± 2137 m, min-max 677-8717 m), and mean pond size is 1399 m 2 (SD ± 895 m 2 , min-max 78-2712 m 2 ).
Two nature reserves (2.89 km 2 and 0.23 km 2 , established in 1967 and 2014) within the study area preserve old natural forest (Fig. 1). A third area (about 3 km 2 ) also contributes old natural forest with documented high conservation values (Fig. 1), including many rare species connected to old natural boreal forest (Reiso 2018). Contrasting forestry related human impacts in the remaining study area ( Fig. 1) mainly consists of clearcutting and associated construction of forest gravel roads, a main power line, and a few scattered cabins and trails. In 1992 a small forest fire burnt an area of 2.25 km 2 ( Fig. 1) (Slettemo 2008).

Sampling
Capture of adult and juvenile newts for DNA sampling was conducted between 20 May and 17 July 2017, using funnel traps (Dervo et al. 2013). Tail clips were stored in 96% ETOH at − 18 °C. The mean sample size from each pond was 31.4 individuals (SD ± 10.3; range 4-39). Permits for capture and sampling of newts and including ethical considerations, were acquired from the County Governor of Telemark (20.02.2017), Norwegian Environment Agency (20.03.2017) and Norwegian Food Safety Authority (08.08.2016, license #9118).

Molecular methods
Genomic DNA was extracted with the Qiagen blood and Tissue kit, following the manufacturer's instructions (Qiagen 2006). Microsatellite markers (Tables SI1.1, SI1.2) were developed from Håland (2017) and Drechsler et al. (2013). For microsatellites from Håland (2017), different primer and universal primer concentrations were tested to develop an optimum combination for Polymerase Chain Reaction (PCR) amplification (Table SI1.1). All PCRs were run with a final volume of 20 µl on an Applied Biosystems 3130xl Genetic Analyzer (https ://www.therm ofish er.com). The microsatellites from Drechsler et al. (2013) were run with the temperature profile prescribed in the Qiagen kit (Qiagen Multiplex PCR Master Mix), and PCR products were diluted with 100 µl dH 2 O before visualization. The locus Tcri46 primer sequences, as described by Drechsler et al. (2013) had to be corrected, because of an apparent mix-up of forward and reverse primers in the original article (Table SI1.2).
Twelve loci were dropped after initial testing due low polymorphism and no or uninterpretable amplification. Locus Tcri29 was amplified for all samples, but later dropped because of difficulties with defining the alleles, leaving in all 14 loci which were used for subsequent analysis (Table SI1.2). All results were analyzed in GeneMapper v5 (AppliedBiosystems). Error rate was estimated by re-amplifying 10% of the samples (arbitrarily selected) and found to be 1.78%.
The 14 amplified microsatellite loci were tested for departure from Hardy Weinberg equilibrium and for linkage disequilibrium, within pond samples, using Genepop v.4.7.0 with 90,000 and 600,000 iterations respectively (Rousset 2008), and significance assessed after sequential Bonferroni correction (Holm 1979). Micro-Checker v2.2.3 was used to test for null alleles, scoring errors and large allele dropouts using 10,000 iterations and α = 0.05 (Van Oosterhout et al. 2004). The presence of candidate loci under natural selection was investigated using BayeScan v2.1 after 5,000,000 iterations following 500,000 burn-ins (Foll and Gaggiotti 2008), and locus Tc50 dropped from further analysis. The sample from pond F was considered too small to provide reliable estimates (n = 4) and omitted from subsequent analysis. This resulted in a total sample size of 404 unique individuals (range 12-39 per population) with 13 markers.

Population structure analysis
We used GenAlEx v.6.503 (Peakall andSmouse 2006, 2012) to calculate number of alleles (N A ), observed heterozygosity (H O ), expected heterozygosity (H E ) and expected heterozygosity corrected for small samples (uH E ). Allelic diversity was calculated as allelic richness (A R ) with the package "diveRsity" in R (Keenan et al. 2013). Allelic richness was calculated with rarefaction, and with and without the smallest of the remaining samples (pond G), giving a per pond sample size of 12 and 32, respectively. Pairwise genetic differentiation was estimated by Weir and Cockerham's F ST (1984) calculated in SPAGeDI v1.5 (Hardy and Vekemans 2002), and significance of results were evaluated after 15,000 permutations with 95% confidence intervals generated by jack-knifing over loci. We also calculated Chord Distance D C (Cavallisforza and Edwards 1967) using FreeNA (Chapuis and Estoup 2007) ( without INA-correction for null alleles), and 95% bootstrap confidence intervals were from 15,000 replicates. Population structure was inferred using STRU CTU RE v2.3.4 (Pritchard et al. 2000) based on variation in allele frequencies, and the minimization of within-population departure from Hardy Weinberg proportions, and linkage disequilibrium (Pritchard et al. 2009). STRU CTU RE has been criticized for not finding the correct population structure, when samples are unbalanced (Kalinowski 2010;Wang 2017). Including the smallest sample (pond G, n = 12, Table 1) in the material, introduced the problem of unbalanced sampling. To consider the above-mentioned issues, we used these recommended settings (Wang 2017): (1) alternative prior (α inferred for each source population), (2) initial α = 1/K = 1/12, and (3) the uncorrelated allele frequency model. STRU CTU RE was then run with 10 replicates for each possible number of clusters (K), using the admixture model, the above settings and 200,000 replications of burn-in and 500,000 MCMC replicates. K was set to range from 1 to 12. The optimum number of clusters was estimated with both: (1) the mean likelihood of the data (mean Ln P(D) (Pritchard et al. 2009), and (2) the ∆K method (Evanno et al. 2005).
The mean Ln P(D) method may be the best method when working with unbalanced samples and the above recommended settings (Wang 2017). The ∆K method is capable of finding the uppermost level when populations are hierarchically structured (Evanno et al. 2005), but primarily works well with balanced samples (Wang 2017). Therefore, STRU CTU RE was also run without the smallest sample (Pond G), with admixture model, correlated allele frequency model, initial alpha = 1.0, 10 replicates per K [1-11] (the smallest sample G with n = 12 omitted to have balanced samples) and 200,000 replications of burn-in and additional 500,000 MCMC replicates. Both methods for optimum K estimation were implemented for both STRU CTU RE runs using STRU CTU RE-Selector (Li and Liu 2018). Optimal alignment of replicates for the same K was obtained for the most relevant K values, and performed in the software Clumpp, with the Greedy algorithm (2000 repeats). Bar graphs of the aligned individual assignments were generated using Clumpak (Kopelman et al. 2015). Emaresi et al. (2011) developed a simple and flexible landscape genetics approach that identify relevant landscape variables influencing population structure, without a priori and potentially unrealistic assumptions about dispersal. By using strips of varying widths that defined a dispersal corridor between pairs of populations of the newt species (Mesotriton alpestris), Emaresi et al. (2011) were able to identify land-uses that acted as dispersal barriers (i.e. urban areas) and corridors (i.e. forests). Using this method, we quantified proportions of relevant landscape variables in strips between populations. A principal advantage of this strip-based method, is that it does not depend on the parameterization of cost values based on a priori assumptions about dispersal strategies or abilities. Emaresi et al. (2011) also tested performance of 11 different fixed (110-510 m) vs. ratio (1:1-1:9) width strip models. The best model, i.e. with width to length ratio of 1:3, is used here. We also tested two contrasting strip width models, i.e. narrower 1:2 and wider 1:5, by comparing mean marginal R 2 across models, using the package "piecewiseSEM" in R (Lefcheck 2016). In accordance with Emaresi et al. (2011), strip width ratio 1:3 performed best, although only marginally better than 1:2, and 1:3 results are reported here. In addition to distance, six more landscape variables were included in the analysis (Table SI3). Variables that could lead to low prey abundance were assumed to affect propensity to move through an area and/or mortality. These were (1) low soil productivity (PROD), because of its negative correlation with newt reproductive success and prey abundance (Wareborn 1992;Vuorio et al. 2013), and/or (2) removal of forest-canopy (OPEN), leading to less prey through a drier microclimate (Stuen and Spidsø 1988;Atlegrim and Sjöberg 1996). Variables that could lead to a drier microclimate may decrease landscape permeability for newts, while moist microclimates may increase it. Landscape variables considered proxies for a drier microclimate with more solar radiation and evaporation, were (3) aspect (ASP; south/south-west facing), particularly in non-forested areas due to the loss of the canopy's shadowing effect (Oke 1987), and also (4) forest gravel roads (ROAD), because of the creation of forest edges and the "edge effect" (Marsh and Beckman 2004). Streams (5) (STRM) may function as humid dispersal corridors (Spear and Storfer 2008). A steeper topography (6) SLOP) per se likely increases the cost of moving through the landscape and may even evoke avoidance behavior. Slopes have thus been found to work as amphibian dispersal barriers (Marsh et al. 2005; Richards-Zawacki 2009). There is little knowledge as to how steep, and it will depend on local landscape features. We tentatively set steep slopes as 20 resp. 30° and steeper. Because they correlated strongly (r = 0.88), we only included the latter. Steeper inclines than 30° are rare in this landscape (2.7%), and were therefore aggregated. All landscape variables were quantified using ArcMap v10.4.1 (ESRI 2015) with 3 × 3 m cell size, and based on Lidar-data, aerial photos and existing land cover data (Kartverket 2008(Kartverket , 2010a(Kartverket , b, 2015. Land cover data (from 2010) were used to create rasters delimiting areas of low soil productivity, gravel roads and streams. Lidar-data (from 2008) were used to create a Digital Terrain Model (DTM) and Digital Surface Model (DSM) of the area. These models were used to calculate slope and delimit areas with no forest cover. The most recently available aerial photos (2015) were used to update information about forest cover (https ://www.nibio .no/tema/ skog/kart-over-skogr essur ser/satsk og). Forest road data were also updated from the recent topographic map (norgeskart. no), but with virtually no change. Spearman Rank Correlations were calculated across all pairs of quantified landscape variables.

Landscape resistance and permeability
Relationships between genetic structure and landscape features, including geographical distance, were explored with a maximum likelihood population effects linear mixed effects model using the package "ResistanceGA" in R (Peterman 2018). Geographical distance and landscape variables were incorporated as fixed effects, with dependency of observations caused by the pairwise study design as random effect (Clarke et al. 2002). Two simple models: one random pattern (only intercept) and one isolation by distance (IBD) only (geographical distance), were tested together with the 16 models combining landscape factors, which were all run with and without distance as fixed effect (Table 2). All models were run with either pairwise F ST (Weir and Cockerham 1984) or D C (Cavalli-sforza and Edwards 1967) to model genetic differentiation as the response variable. All predictor variables were standardized before the models were run (Schielzeth 2010).
Models were compared using an information-theoretic model approach with a correction for finite sample size (AICc) (Akaike 1974;Burnham 2002), exercising caution in interpretation if models include additional parameters, but are within 2 AIC units of the top ranking model (Arnold 2010). AICc, ∆AICc, Akaike's weights (= w) and log likelihood were calculated using the package "AICcmodavg" in R (Mazerolle 2017). Models run with the response variables Fst resp. Dc were compared separately. According to Burnham (2002) models with a ∆AICc < 4-7 are plausible models. Relative support of models with ∆AICc < 4 were evaluated using the pairwise evidence ratio = w i /w j (Burnham, 2002).
The strength of the relationship between predictors and response variable, i.e. effect size, was expressed using standardized estimated regression coefficients for the highest ranked models. Standardized regression coefficients can be used to compare effect size across alternative models, and even between studies (Schielzeth 2010). The uncertainty of the estimated regression coefficients was assessed with approximate 95% confidence intervals, attained by the ruleof-thumb formula: ± 2 × Standard Error.

Impact of old forest on genetic and allelic diversity
Potential population differences in estimated genetic diversity (observed heterozygosity H O ) and allelic diversity (allelic richness A R ) between ponds in areas of old forest compared to those in managed forests, were tested with a two-sided permutation test implemented in FSTAT v2.9.3.2 (15,000 permutations, Goudet 2001). Ponds located entirely within (A, G, K, L) or bordering on old forest (< 50 m; B, H) were considered 'old forest' (group 1). Areas of old forest were the two nature reserves and the area of old natural forest with documented high conservation values (Fig. 1). Ponds in managed forest (group 2) were C, D, E, I, J and M, all with a distance greater than 390 m from old forest (Fig. 1). Allelic richness was estimated with and without pond G, as its outlying small sample size might bias the subsample size used across all samples to estimate allelic richness. The subsample size used to calculate allelic richness was therefore 12 resp. 32. Effect sizes were evaluated using Cohen's d (Li 2010).

Genetic variation
Levels of genetic variation for the microsatellite loci varied among the sampled populations ( Table 1). The number of alleles (N A ) showed a total mean of 3.9, and ranged from 3.2 (pond G) to 5.2 (pond K). Unbiased expected heterozygosity (uH E ) ranged from 0.408 (pond I) to 0.539 (pond L). Pond I also had the lowest allelic richness (A R = 2.9/3.2), and pond K the highest (3.7/4.5; Table 1). No significant deviation from Hardy Weinberg equilibrium was found within the samples (P > 0.0036; 11 of 168 tests were significant. before correction, affecting 6 loci and 7 samples). One locus (Tc50) appeared to be under balancing or purifying selection-(P (α ≠ 0) = 1.00). Also, one loci pair Tcri36-Tc50 in pond M was found to be in linkage disequilibrium (P = 0.000003, 67 of 1183 tests were significant before correction, affecting all loci and samples). Thus, Tc50 was dropped from subsequent analyses.

Population structure
Both F ST and D C indicated structure across all sampled pools (Table 2, significant pairwise tests), STRU CTU RE suggested an overarching structure with number of clusters substantially lower than number of sampled ponds. The optimum number of clusters estimated from STRU CTU RE with the ∆K procedure was four clusters (Fig. 2), and for both balanced and unbalanced samples. The alternative mean Ln P(K) method predicted a higher number of clusters both with balanced and unbalanced sampling. This number was somewhat lower with the uncorrelated than with the correlated allele frequency model (optimum K = 8 and 9-10, respectively).
Based on the lower and likely more accurate ∆K result (Pritchard et al. 2000(Pritchard et al. , 2004Bergl and Vigilant 2007), three main clusters were located in the best model, one in the south around the nature reserve (A, B, C, E, G, H, (SD ± 1071 m), Fig. 3), one in the east/north-east (I, M (dist. = 1543 m)), and one in the north/north-west (J, K, L, (SD ± 448 m)) inside or close to the area with old natural forest (Figs. 1, 3). The pond at the forest fire site (D) constituted a fourth cluster on its own. The geographical distance between neighbor ponds within clusters (mean distance 2048 m, 1543 m and 1744 m, respectively), was sometimes longer than between neighbor

Landscape resistance and permeability
Of the 34 models tested (Table 3), landscape variables associated with microclimate affected permeability for newts the most (Tables 4, 5), and with very similar results with either F ST or D C as response variables, For both F ST and D C two top models clustered within ΔAICc less than 2 and four models within ΔAICc less than 4 (Table 4). All top models included geographic distance, aspect and gravel roads as important landscape factors (Tables 4, 5). Adding streams appeared to improve the top ranked model somewhat, suggesting that streams may perhaps increase landscape permeability for Northern newts (Tables 4, 5). The effect of OPEN was included in two of the four top models, but did not improve performance of the two top models (Table 4). The topographic factor slope did not contribute to any of the best models (Tables 3, 4. The only notable difference between F ST and D C as response variable, might be that models ranked 1 and 2 were closer for F ST (Table 3), but the range is less than 2 and likely is insignificant. As expected, some of the landscape factors showed some correlation (Table 6). Notably, aspect and gravel roads were moderately correlated (ρ = 0.45), as were aspect and streams (ρ = 0.54) and gravel roads and streams (ρ = 0.60) ( Table 6). The simple random and IBD models both ranked considerably lower than the more complex landscape variable models, indicating the importance of landscape features. However, notably the isolation by distance model (Table 3; model 02), performed substantially better than the random model (Table 3; model 01)) ( Table 4).
Relationships among landscape predictors and population genetic differentiation (F ST and D C ), as evaluated by the standardized (z-transformed) regression coefficients from the top ranked models, all indicated a positive relationship with genetic differentiation, but negative for streams (Tables 4,5). Geographic distance, aspect, and gravel roads were all important with substantial effect sizes. For example, for the best model (M3d) with F ST as response, an increase from 0 to 2.6% (the most extreme value observed in our study; ~ 6 SD) area of gravel roads, could increase F ST from 0.100 to somewhere between 0.152 and 0.269, i.e., from moderate to great, or very great genetic differentiation (Hartl and Clark 2007). A similar increase in F ST by geographical distance would require a distance of 16.5 km (= 7.7 SD).

Impact of old forest on genetic and allelic diversity
Old forest ponds had a significantly higher allelic richness compared to the other ponds in managed forest (average A R = 3.63 and 3.24, respectively, two-sided permutation test, P = 0.030), and with a substantial effect size (Cohen's d = 1.15) (Sullivan and Feinn 2012). However, the difference in observed heterozygosity between old forest ponds and the other ponds, was not significantly different (average H O = 0.509 and 0.488, respectively, two-sided permutation test, P = 0.372). Excluding pond G from the analysis did not affect the result significantly (average H O = 0.503 and 0.491, respectively, two-sided permutation test, P = 0.686).

Discussion
Habitat loss and changed composition reduce, and represent a substantial threat to, many amphibian populations (Cushman 2006). Conservation efforts focus on the importance of avoiding the negative genetic effects of small population sizes, e.g. genetic drift and inbreeding (Frankham et al. 2010). For proactive management and conservation, it is important to know what may enhance, or reduce these effects. This could be specific for the type of organism, the species and certainly the type of habitat the populations inhabit (Sih et al. 2000;Keyghobadi 2007). In our study on the Northern crested newt populations in the boreal forest ecosystem exhibited strong genetic structure with substantial diversity. The landscape features aspect (south-southwest facing) and gravel roads, and likely also streams, together with, geographical distance best explained genetic differentiation. In addition, ponds connected to old forest showed a higher allelic diversity compared to ponds in areas more affected by forestry.

Landscape effects on genetic differentiation
In the analysis of landscape effects on genetic differentiation, isolation by distance hypothesis was used as a simple model. However, this model performed poorly relative to models which also included landscape variables. Clearly, the intervening landscape between newt breeding ponds plays an important role in shaping the pattern of genetic differentiation. This was also consistent with the results from STRU CTU RE. Here, geographical distance between ponds in the same clusters was sometimes much larger than between ponds from different clusters. All the top models, however, included both landscape factors and geographical distance, indicating their additive effects.
Forest gravel roads emerged as an important landscape factor in addition to distance. Roads as such, have been found to be a barrier for amphibians, e.g. the smooth newt (Lissotriton vulgaris) (Sotiropoulos et al. 2013) and the spotted salamander (Ambystoma maculatum) (Richardson 2012). Roads also directly increase mortality by vehicles killing or injuring the amphibians (Mitchell et al. 2008). This mortality depends on amount of traffic in relation to activity periods and behaviors of the relevant amphibian species (Hels and Buchwald 2001). The forest gravel roads in Fig. 3 Results from STRU CTU RE, run with the uncorrelated allele frequency model and unbalanced sample sizes (see Fig. 1 for symbol legend). The pie graphs indicate STRU CTU RE results by sampled pond, and for 4 clusters (Fig. 2a) and 8 clusters (Fig. 2b), respectively. The lower bar graph shows individual assignment probabilities grouped by sampled pond. Each horizontal bar chart represents number of clusters (K), from 4 through 8 Table 3 A priori models and hypothesis' of the relationship between predictor landscape variables and geographical distance, and response genetic differentiation measured as F ST or D C Expected negative effects are denoted with (−), and positive effects with (+). The seven modelled predictors were geographical distance = DIST, low soil productivity = PROD, non-forested areas 2015 = OPEN, gravel roads = ROAD, aspect (south/south-west) = ASP, streams = STRM, and slope 30° or steeper = SLOP the study area provide access to outdoors activities (fishing, hunting, hiking), and sporadically for logging. Thus traffic is limited and mostly daytime, while the northern crested newt is mostly active during night (Dervo and Kraabøl 2010). Therefore, it is unlikely that traffic-related mortality caused the relationship of gravel roads on genetic differentiation in our study. A more likely explanation comes from the microclimatic effects of forest. Construction of roads entails forest removal and creation of forest edges. Forest edges lead to increased evaporation; when air from hot and dry surfaces meets a more humid vegetated surface (Oke 1987). Marsh and Beckman (2004) found that in Virginia, USA, forest redbacked salamanders (Plethodon cinereus) were observed less frequently in the edges around forest roads (up to 20 m from the roads), due to decreased soil moisture. In an experimental study, Marsh et al. (2005) found that when they displaced red-backed salamander and measured their return rates, forest roads (width 5-8 m) reduced the return rate by 51%. The landscape factor south/south-west facing slopes, moderately correlated with gravel roads, explained much of the increased genetic differentiation. Since many amphibians, including the genus Triturus, have a poor ability to regulate water loss (Wells 2007), drier microclimates likely pose a problem, and a general preference for moist and cool microclimate amongst amphibians (Harper and Guynn 1999). High solar radiation load has also been found to affect gene flow negatively, e.g., in the southern torrent salamander (Rhyacotriton variegatus) in California (Emel and Storfer 2015) and the coastal tailed frog (Ascaphus truei) in Washington (Spear and Storfer 2008). Conversely, the stream variable presumably providing more suitable dispersal habitat for newts, appeared to be related to a decrease in genetic differentiation. The model including the stream variable had 2.7 times more support (D C as response variable), compared to the same model without stream. Although the results indicate that streams can have an effect, but the amount is uncertain due to large confidence intervals, and needs further study. For the southern torrent salamander (Rhyacotriton variegatus) streams have been shown to increase gene flow (Emel and Storfer 2015). On the other hand, larger waterways such as rivers have been found to limit gene flow for both the spotted salamander and the northern crested newt (Peter et al. 2009;Richardson 2012). Rivers and streams constitutes a dynamic continuum between barriers and corridors (Puth and Wilson 2001). Whether the water way is experienced as barriers or corridors, depend on species mobility (vagility and ability) relative to width and length of the water way and the amount and gradient of water flow (Puth and Wilson 2001). In addition, in a managed forest in Idaho, USA, the Rocky Mountain tailed frog (Ascaphus montanus) used streams as corridors. This was not found for frog populations in a control area less impacted by forestry (Spear and Storfer 2008). In other words, streams may have different effect on gene flow in amphibians depending on stream characteristics or the surrounding environment. In our study, the effect of streams correlated with gravel roads, likely because they follow natural valleys in a rather hilly landscape. This likely confounds any effect of streams per se.
Unexpectedly, open non-forested areas, i.e., clear-cuts, did not turn up as an important landscape variable affecting genetic structure in our study. Clear-cutting seems to affect many amphibian species negatively (Semlitsch et al. 2009;DeMaynadier and Hunter 2011;Tilghman et al. 2012). Northern crested newts tend to avoid clear-cuts and other non-forested areas in their terrestrial habitat around the breeding pond (Kupfer and Kneitz 2000;Vuorio et al. 2015). Forest cover has also been found to be positively related to gene flow for several salamander species (Greenwald et al. 2009;Emaresi et al. 2011;Richardson 2012;Emel and Storfer 2015). However, most of these studies have tested the effect of forest in contrast to other landscape types, such as agricultural areas, open fields, etc. This would represent more stable landscape patterns than the dynamic of nonforested areas caused by forest harvest, and ensuing diverse successional processes at different rates (Sih et al. 2000;Keyghobadi 2007). Besides, an area could be used for dispersal for years, but after logging suddenly become a partial barrier. The effect would, however, not be immediately detectable in genetic data because of time lag. For example, the time lag between forest harvest and genetic response (measured as G'st) was 20-40 years for the coastal tailed frog (Ascaphus truei) (Spear and Storfer 2008).

The impact of old forest
In our study, allelic diversity, measured as allelic richness, was significantly higher in the ponds grouped as old forest ponds. Areas of old forest represent less fragmented habitat, in time and space, compared to the rest of the study area. This is also congruent with the results from STRU CTU RE, where populations connected to old forest remained clustered, while populations not connected to such areas constituted more or less separated populations (8 clusters). It is also supported by the fact that the pond most affected by the loss of canopy cover (pond D in the 1992 forest burn), was singled out as one single cluster in STRU CTU RE at the level of just 4 clusters. Although elevated levels of genetic diversity and similar STRU CTU RE clustering patterns are not necessarily related, it does not appear unlikely that similarly high diversity among close sites stem from shared alleles.
Several studies have found lower genetic and allelic diversity in amphibian populations located in fragmented habitats, e.g. Cosentino et al. (2012); Hitchings and Beebee (1997); Johansson et al. (2005). However, these studies focused on the effect of habitat changes caused by urbanization or agricultural land use, and not forest as in this study. On the other hand, a study in British Colombia, Canada, on the coastal giant salamander (Curtis and Taylor 2004) found that allelic richness and heterozygosity was positively correlated with the age of the forest stands and allelic richness showed a higher correlation (r 2 = 0.59) compared to heterozygosity (r 2 = 0.37). Old growth forest showed the highest amount of both measures of genetic variability. The results were attributed to the negative effects of forestry harvest on population sizes, leading to higher impact of genetic drift (Curtis and Taylor 2004).
We found a significant relationship between old growth forest and allelic richness. However, genetic diversity expressed as observed heterozygosity was not significantly different between old forest and more impacted habitat. This discrepancy may be an artifact of time elapsed since disturbance. Forestry activities may have led to reduced population sizes and/or migration rates, making the populations more exposed to the effect of genetic drift. Alleles are lost more rapidly than heterozygosity when population sizes are reduced (Hedrick 2011), and it takes less time before the effect of some disturbance affects the number of alleles, compared to heterozygosity (Lloyd et al. 2013). Another possible scenario is that populations outside of old forest areas more frequently experience events of low population size, but grow fast enough so that heterozygosity is not substantially affected (Nei et al. 1975;Allendorf et al. 2013). Loss of heterozygosity depends not only on population reduction, but also on the duration of the period of low population size, whereas allelic richness is more connected to the size of the reduced population (Hedrick, 2011). Both scenarios, though, would lead to a loss of alleles and thus evolutionary potential (Caballero and García-Dorado 2013).

Conclusion
We found that landscape composition factors associated with microclimate and geographical distance between ponds affected the genetic differentiation of northern crested newt populations in a northern boreal forest ecosystem Moreover, we found that populations located within or near old forest exhibited a significantly higher allelic diversity, relative to populations in more managed forest, i.e., more affected by forestry harvest and road constructions.
Both habitat loss and habitat changes can be a challenge for many animal populations (Cushman 2006). A prerequisite for sustainable populations in affected areas is the maintenance of gene flow. This requires pro-active knowledge about the effects of the intervening landscape on dispersal, which we provide here for amphibians in boreal forest ecosystems. Furthermore, such baseline knowledge can be used to evaluate the future likely viability of discrete northern crested newt populations, or in a meta-population framework, by empirically parameterizing model parameters. Such, data can potentially be used in virtual population models projecting and evaluating the viability of genetic population units into future given current and future environmental challenges.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.