Past forest cover explains current genetic differentiation in the Carpathian newt (Lissotriton montandoni), but not in the smooth newt (L. vulgaris)

Current genetic variation and differentiation are expected to reflect the effects of past rather than present landscapes due to time lags, that is, the time necessary for genetic diversity to reach equilibrium and reflect demography. Time lags can affect our ability to infer landscape use and model connectivity and also obscure the genetic consequences of recent landscape changes. In this work, we test whether past forest cover better explains contemporary patterns of genetic differentiation in two closely related but ecologically distinct newt species—Lissotriton montandoni and L. vulgaris.


| INTRODUC TI ON
The interpretation of landscape genetic patterns from rapidly changing environments is not always straightforward, as such patterns may reflect past rather than current population parameters, such as effective size or migration rate (hereby time lag; Epps & Keyghobadi, 2015). In dynamic landscapes such time lags may be the rule rather than the exception. Considering this, some landscape genetic studies have taken advantage of spatial data collected at multiple time intervals to assess whether genetic patterns more strongly reflect past rather than contemporary landscapes (Anderson et al., 2010;Epps & Keyghobadi, 2015). Detecting and estimating time lags is crucial not only for understanding the limits of predictions and inferences made from studies based on contemporary landscape data but also for recognizing the effects of recent landscape changes (Epps & Keyghobadi, 2015;Landguth et al., 2010).
Importantly, if such time lags exist, they provide a time window for preserving the existing genetic diversity through rapid restoration activities (Reinula et al., 2021).
Despite the acknowledged need to consider time lags in landscape genetic studies, their prevalence in nature is not yet known (Manel & Holderegger, 2013). Simulation studies suggest that the genetic effects of landscape changes can be rapidly discerned (1-15 generations; Landguth et al., 2010;Prunier et al., 2014), but many factors (e.g. generation times, mutation rates, dispersal rates and distances, nature of landscape changes, effective population sizes, etc; Epps & Keyghobadi, 2015;Keyghobadi et al., 2005;Landguth et al., 2010) may delay their onset and lead to lags of more than 100 generations . Variation in time lags is reflected in empirical studies, with evidence from multiple taxa (e.g. bush cricket, Holzhauer et al., 2006; costal tailed frogs, Spear & Storfer, 2008; copperhead snake, Maigret et al., 2020) ranging from ~2 to 50 generations. Many studies also show that genetic divergence can develop rather quickly in response to strong barriers to gene flow (Ascensão et al., 2016;Holderegger & Di Giulio, 2010;Keller et al., 2004;Mapelli et al., 2020), and commonly, the observed genetic patterns reflect contemporary landscapes rather than historical ones (Clark et al., 2010;Crossley et al., 2019;Winiarski et al., 2020;Zellmer & Knowles, 2009). Overall, there is a clear need for studies investigating time lags in nature and discussing their implications for landscape genetic methods and biodiversity conservation (Manel & Holderegger, 2013).
The Carpathian (Lissotriton montandoni) and the smooth newt (L. vulgaris) are closely related (Pliocene divergence; Pabijan et al., 2017;Zieliński et al., 2016), but ecologically distinct salamandrid species co-distributed along the Carpathian Mountain range (see inset map Figure 1a; Wielstra et al., 2018). Lissotriton montandoni is distributed at higher elevations (up to 2000 m) and is found mainly across humid forests, forest edges and nearby meadows while L. vulgaris is distributed at lower elevations, generally occurring in woodland habitats, but also in meadows and a wide variety of human-modified habitats (e.g. parks, orchards in rural and urban areas). The species have experienced a long history of hybridization and interspecific introgression Zieliński et al., 2013Zieliński et al., , 2019. However, given their relatively well understood and restricted geography of hybridization, it is possible to sample populations of each species that have F I G U R E 1 (a) Sampling localities for Lissotriton montandoni (red) and L. vulgaris (green) across the Carpathian Mountain range. Point size represents the number of genotyped pure (min. 97% ancestry from one species as inferred by ADMIXTURE) individuals per locality and species. Black triangle shows North. Inset map shows species distribution ranges. (b) Population structure (ADMIXTURE K = 2) for L. montandoni. (c) Population structure (ADMIXTURE K = 3) for L. vulgaris. Transparent grey polygons outline the study areas explored in landscape genetic analyses. been largely unaffected by introgression of the nuclear genome for use in a landscape genetics framework (Antunes et al., 2022). Antunes et al. (2022) found that anthropogenic landscapes negatively affect both species, reducing habitat connectivity and increasing genetic differentiation between populations. Still, populations of L. vulgaris tend to show higher genetic differentiation and lower census and effective sizes, suggesting stronger isolation as a possible result of exposure to recent habitat destruction and fragmentation.
Resistance models demonstrated that forest promotes population connectivity in both species, but also revealed differential use of forested habitat, with L. montandoni and L. vulgaris showing the highest population connectivity at forest core and forest edges, respectively (Antunes et al., 2022). However, given the substantial changes in forest cover and fragmentation in the Carpathians over the last five decades (see Figure S1), contemporary genetic diversity may still reflect historical landscapes, and this could confound landscape genetic inference and misinform conservation actions. Importantly, time lags may differ considerably between species. Contrasting time lags may arise not only due to differences in effective population size (lower Ne can reduce time lags), but also due to the potentially different speed, direction and extent of changes in key habitats (Epps & Keyghobadi, 2015).
In this work, we test whether past forest cover better explains contemporary patterns of population genetic differentiation in L. montandoni and L. vulgaris. We explore such time lags using historical forest-cover data from 1963 to 2015 (ca. 13 newt generations). We expect to find evidence for time lags due to the recent forest-cover changes around the Carpathians ( Figure S1). Time lags of different duration are expected for each species due to their demographic disparities, dissimilar ecologies and unequal exposure to land-cover changes. Exact predictions are, however, difficult to make, because different factors can push time lag duration in opposite directions.
Shorter or even absent time lags can be predicted in L. vulgaris, due to its smaller population sizes, which accelerates drift and/or generalist ecology that may make it more tolerant to land-cover changes.
But longer time lags in L. vulgaris can also be hypothesized, due to its potentially higher exposure to land-cover changes, in comparison with L. montandoni. Thus, here, we compare the time lag duration between species, and only later discuss the likeliness of scenarios as the ones described above.

| Sampling, sequencing and filtering
Sampling was done in the Carpathian Mountains and their foothills, covering most of the distribution of L. montandoni but restricted to the areas of species co-occurrence for L. vulgaris (see Figure 1a).
Adult newts were sampled in water during the breeding season (March-May) with most regions visited during a single breeding season but some over the course of several years (2010-2019).
Individual ponds, or sets of adjacent (within a radius of 100 m) puddles, were treated as local populations. Tail-tip biopsies were taken, and animals were released afterwards. Tissues were stored in 96% ethanol and DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega).
Targeted resequencing was done using molecular inversion probes (MIPs), an approach previously optimized and successfully used in Lissotriton newts, following the protocol of Niedzicka et al. (2016). MIPs are single-stranded DNA molecules containing at their end sequences complementary to the two regions flanking a target. After hybridization of MIPs to the target, gap-filling and ligation result in circularized DNA molecules that contain the target sequence together with adaptors. Barcoded primers complementary to the adaptors are used to obtain linear DNA fragments ready for downstream analyses (Niedzicka et al., 2016). MIPs were designed to target 112 bp from protein-coding exons and 3′UTR regions of Lissotriton transcripts . In total, we About half of the samples (1723 from a total of 3544, Table S1) had already been collected for a study of hybrid zones between L. montandoni and L. vulgaris (see Zieliński et al., 2019). However, because of extremely limited hybridization and introgression between species (the width of the average allele-frequency cline is 0.1-5.1 km, and even in the centre of the zone parental genotypes predominate), these data are well suited for landscape genetic analysis (see Antunes et al., 2022). Nonetheless, we adopted a conservative approach and estimated the proportion of L. vulgaris and L. montandoni ancestry per individual using ADMIXTURE (Alexander et al., 2009) with K = 2 and removed from the initial dataset individuals with more than 3% admixture from the other species (Table S1).

| Assessing phylogeographic structure
Three major evolutionary units were identified in previous works: L. montandoni, occurring along the Carpathian Mountain range; and two lineages of its sister species, L. v. ampelensis and L. v. vulgaris, which occur inside and outside the Carpathian Basin, respectively. Here, we assess phylogeographic structure by applying ADMIXTURE analyses at the species level evaluating the number of clusters K up to 10. These results confirmed the presence of the previously defined groups (based on cross-validation error and Fst among the main clusters) and revealed further levels of substructure that need to be accounted for in the landscape genetic analyses (see Section 2.5).

| Genetic differentiation, diversity and isolation by distance
For each evolutionary unit (L. montandoni, L. v. ampelensis and L. v. vulgaris; confirmed by phylogeographic analyses, see results below) we estimated genetic differentiation between populations using pairwise Fst (Nei & Li, 1979) as calculated in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). We also obtained genetic diversity measures that could provide important hints about past and more recent evolutionary history. Specifically, we used the R package 'HierDpart' (Gaggiotti et al., 2018) to calculate three information-based diversity measures based on Hill numbers of different orders that give different weights to common and rare alleles (Jost, 2006): richness (effective number of alleles; q = 0) being completely insensitive to allele frequencies; exponential Shannon entropy (q = 1) weighing the contribution of each allele by its frequency without favouring either common or rare alleles; and a heterozygosity-related measure (q = 2) favouring the most common alleles (Gaggiotti et al., 2018).
Results from genetic differentiation (pairwise Fst) were used to test for isolation by distance (IBD; Wright, 1943) in each evolutionary unit, including only pairs of populations at distances from 1 to 50 km (threshold used in landscape genetic analyses, see below), but also for pairs of populations at lower distances (0-5, 0-2.5 and 0-1 km), because such analyses can provide key information about species dispersal. This was done using the maximum likelihood population effects model (MLPE; Clarke et al., 2002). For implementation, we used the R package 'corMLPE' (https://github.com/nspop e/corMLPE; Clarke et al., 2002), which together with the R package 'emmeans' (Lenth et al., 2018) allowed us to compare IBD trends between evolutionary units. Patterns of IBD were visualized using the R package 'ggplot2' (Wickham, 2016).

| Present and historical land-cover data
Land-cover data for the present (2015) was downloaded from https://lcvie wer.vito.be/ at ~100 m spatial resolution (Buchhorn et al., 2019). These global land-cover maps are part of the Copernicus Land Service, derived from PROBA-V satellite observations. Rasters were resampled to a resolution of ~1 km, which was the finest resolution available for historical data (see below). Dispersal events in Lissotriton vulgaris and other small-bodied newts (e.g. Ichthyosaura alpestris) are more frequent at small geographic distances (<500 m; Bell, 1977;Kovar et al., 2009;Schmidt et al., 2006). However, occasional migrants may connect populations separated by larger distances (>1 km; Jehle & Sinsch, 2007;Schmidt et al., 2006;Alex Smith & Green, 2005). Connectivity between populations at distances much larger than 1 km is also possible in a stepping-stone manner (Saura et al., 2014). Raster cells of 1 km resolution, despite obscuring the role of smaller landscape elements (e.g. roads, small rivers), are still expected to capture the effect that larger landscape features (e.g. forest, croplands and urban areas) have on connectivity (Anderson et al., 2010). Nonetheless, we conduct post-hoc analyses comparing the best isolation by resistance (IBR) models built at 1 km resolution with their finer resolution versions (500 and 250 m) for the present day landscape data. Original thematic resolution was also changed: we reclassified the different categorical variables (UN-FAO's Land Cover Classification System) to represent current land cover with two different levels of complexity. One including forest, croplands and urban areas; and another also differentiating between open-and closed-forest areas, because of the contrasting response to such landscapes found in a previous study of L. montandoni and L.
Historical land-cover data were collected from https://doi. we adopt a conservative approach and reconstruct historical landscapes using only two categories, forest vs non-forest, because of the overall higher thematic confidence (i.e. probability of identifying the correct land-use/cover class) for forest cover across our study areas (confidence ranging from 72% to 82%; see Figure S2) and because forest cover has been identified as an important driver of connectivity in L. montandoni and L. vulgaris (Antunes et al., 2022).
We also restrict our historical analyses to 13 timeframes, collecting forest-cover data every 4 years (roughly a Lissotriton newt generation) from 1963 to 2015.
For comparison with the results of time lag investigations (see below), we characterized changes in forest cover and fragmentation for study areas specific for each evolutionary unit. This was done using the R package 'landscapemetrics' (Hesselbarth et al., 2019) and the forest-cover data for the same 13 timeframes (every 4 years from 1963 to 2015). These study areas are displayed by polygons in Figure 1 and were built first by creating paths that link populations of the same evolutionary unit between 1 and 50 km and then by creating a 10 km buffer around these paths.
The original 'ResistanceGA' approach transforms landscape surfaces into resistance surfaces that provide the best fit to genetic data (Peterman, 2018). Originally, the range of transformations explored for categorical data is only restricted by minimum and maximum values and uses an optimization algorithm that mimics basic principles of biological evolution and natural selection to find the best resistance surface (R package 'GA'; Scrucca, 2013). Here, due to the large study area (291,282 km 2 with a raster resolution of ~1 km) and number of populations (L. montandoni: 174, L. v. ampelensis: 37 and L. v. vulgaris: 63), we restricted transformations to possible combinations of three values of resistance (1, 10 and 100) and selected the best model a posteriori (using the Akaike information criterion, AIC). This approach is computationally feasible for our dataset and still captures stronger patterns (i.e. differences in resistance of one or two orders of magnitude). Moreover, modelling and model selection were done considering the phylogenetic (evolutionary units and admixture), geographic (variable sampling intensity) and temporal (multiple timeframes) properties of our data (see details below).
IBR regression models was constructed using the MLPE approach (Clarke et al., 2002), which accounts for nonlinearity and allows missing pairwise comparisons (R package 'corMLPE'; https://github. com/nspop e/corMLPE). In these models, genetic differentiation (dependent variable) was related not only with landscape resistance (obtained from transformed original variables using the commu-teDistance function, R package 'gdistance'; Etten, 2017)  ampelensis, and this region is the size of only 12% of our study area; we randomly subsampled pairs of populations to such proportions (12% from Podkarpackie and 88% other regions). The resulting dataset was then further subsampled (85%) to avoid biases created by potential outliers. Both subsampling processes involved 100 iterations and calculation of the mean AIC. Model selection was also done separately per timeframe for the investigation of time lags (see details below).

| Investigating time lags
The investigation of time lags in patterns of genetic differentiation was done by comparing the fit of IBR models (see above) for each evolutionary unit (L. montandoni, L v. ampelensis and L. v. vulgaris) and timeframe (every 4 years from 1963 to 2015; i.e. 13 timeframes). If time lags exist, model fit is expected to increase in the past. Visualization of model fit (∆AIC) through time was done using ggplot2 (Wickham, 2016). In the visualization, we display only models with highest fit per year together with the IBD model fit baseline.
Upon detection of time lag signatures, post hoc analyses using the best IBR model were used to explore the nature of landscape changes (i.e. increased vs decreased connectivity) and groups of populations (e.g. Podkarpackie vs other regions) involved in the time lag.
We plotted changes in resistance (for the best model only) for each type of landscape change and for different groups of populations. If, on the contrary, IBR for the present was the best model, we compared it visually with IBD models.

| Samples, genomic data and phylogeographic structure
The final, filtered dataset, consisted of 2697 samples (Figure 1a

| Genetic diversity patterns
Patterns of genetic diversity showed the expected differences between measures (richness, exponential Shannon entropy and heterozygosity-related measure; Figure S4; Table S2) and phylogeographic groups (i.e. L. montandoni, L. v. ampelensis and L. v. vulgaris), with higher differences in richness than in the other two measures, due to the higher relative weight given to rare variants. Between species/evolutionary units, L. vulgaris showed overall significantly (p-value <.05) higher diversity than L. montandoni. This pattern seems to contrast with L. vulgaris smaller census and effective population sizes documented for a subset of local populations across the same study area (Antunes et al., 2022). However, this higher diversity in L. vulgaris is probably the result of its widespread distribution, historical connectivity and high long-term population sizes, while recent declines at local scales are only visible in the field and using genetic estimates of contemporary local Ne, which are not yet reflected in measures of genetic diversity. This higher diversity in L. vulgaris is especially pronounced for more stable measures (i.e. exponential Shannon entropy and heterozygosity). Mean diversity was lower for richness in L. v. ampelensis than in L. v. vulgaris, while the inverse pattern was seen for exponential Shannon entropy and the heterozygosity-related measure. However, we note that the differences between L. vulgaris groups are significant only for the heterozygosity-related measure (p-value = .027).  Table S3). Notably, very low Fst values were found in L. montandoni across all distances, in contrast with L. vulgaris in which low levels of genetic differentiation were mostly absent for distances higher than 5 km, especially in L. v. ampelensis where Fst values close to zero were found very rarely and only at distances lower than 2 km ( Figure S5). For populations at distances between 0 and 5 km, IBD was significant for all evolutionary units (Table S4), being weaker in L. montandoni, followed by L. v. vulgaris and L. v. ampelensis, which again showed the steepest IBD ( Figure S5; Table S3). For populations at distances between 0 and 2.5 km, IBD was significant for L. montandoni and close to significant in L. v. vulgaris, while L. v. ampelensis showed no visible trend ( Figure S5; Table S4). IBD was not detected at distances of 0 to 1 km ( Figure S5; Table S4).

| Genetic differentiation, isolation by distance and time lags
In L. montandoni, IBR models predicted population genetic differentiation (pairwise Fst) better than IBD models. The best IBR model was a model giving non-forested areas resistance an order of magnitude higher than forested areas (i.e. forest = 1 and non-forest = 10; Figure 3). Versions of this model based on finer resolution landscape data (500 and 250 m) showed poorer fit, supporting the use of 1 km resolution (Table S5). This model was recovered as the best for all timeframes, but model fit increased in the past in a stable trend that fits changes in forest cover and fragmentation ( Figure S1; Figure 3).
Model fit peaked at the 1975 timeframe (40 years before 2015; ca. 10 Lissotriton newt generations). In both L. vulgaris evolutionary units, genetic differentiation was better predicted by present landcover models with the lowest resistance given to open forests (one order of magnitude lower than all other land-cover types, Figure 3).
We also note a much lower model fit in L. v. vulgaris IBR models ( Figure 3), meaning that the Euclidean geographic distance (IBD models) was almost as good in explaining pairwise population genetic differentiation in this evolutionary unit. The similarity between the best IBR model and IBD was also clear when comparing models visually ( Figure S6). None of the historical forest-based IBR models explained the genetic differentiation of L. vulgaris better than IBD.
Post-hoc investigation of time lags was done for L. montandoni using the best IBR model (i.e. forest = 1 and non-forest = 10) to investigate changes in resistance through time. Results show that changes in resistance since 1975 involved both fragmentation and connection of populations ( Figure S7), across different parts of the range ( Figures S8 and S9). However, the most predominant changes involved increased connectivity among populations within the Podkarpackie region ( Figure S9).

| DISCUSS ION
Current genetic variation is expected to reflect past rather than present landscapes due to time lags, that is, the time required for genetic variation to reflect the new demographic and connectivity  Holzhauer et al., 2006;Maigret et al., 2020;Spear & Storfer, 2008) and the many factors influencing them (Epps & Keyghobadi, 2015). However, if we consider species ecology, together with land-cover change (i.e. type, extent and direction) and methodological approaches (e.g. type of spatial and genomic data), we can better understand time lags in nature and their implications for landscape genetic studies and conservation.

| Evidence for time lags (or lack thereof) in L. montandoni and L. vulgaris
The time lag observed in L. montandoni suggests that the full genetic consequences of forest-cover changes from the last 40 years have not yet emerged. Post-hoc analyses show that this time lag signature has been created by both deforestation and reforestation processes, but mostly by the latter, with populations from southern Poland (Podkapackie region) showing increased genetic differentiation, which better fits the higher resistance created by lower forest cover of this region in the past ( Figure S9). This time lag was expected ; however, due to the limit imposed by our historical forest-cover data (1963, ~13 newt generations), we cannot exclude it extending further into the past. Ultimately, the observed time lag results from a balance between many factors that affect the rate at which genetic differentiation approaches equilibrium (Epps & Keyghobadi, 2015). In L. montandoni, higher local effective population sizes (Antunes et al., 2022) probably lead to longer time lags via decelerated drift. However, other factors may have shortened its duration: increased connectivity in Podkarpackie, which is expected to be seen at the genetic level sooner than in scenarios of habitat fragmentation; the type of genomic data used here (thousands of SNPs) that can increase power to detect effects of recent landscape changes on genetic connectivity (McCartney-Melstad et al., 2018).
The lack of evidence for a time lag in L. vulgaris is also not unexpected and seems to be driven by species-specific factors, as current open forest was the best predictor of genetic differentiation in two independently analysed evolutionary units (L. v. ampelensis and L. v. vulgaris). The lack of a time lag in L. vulgaris could be caused by multiple factors acting alone or simultaneously, for example, habitat use, dispersal rates and distances, population sizes and exposure to land-cover changes (Epps & Keyghobadi, 2015). One possibility is that historical landscapes that only distinguish between forested/non-forested habitat do not capture key habitats features of L. vulgaris, as suggested by the low predictive power of current forest-based models (lower than IBD models) and by the higher predictive power of models that distinguish between open and closed forest (as seen here and in Antunes et al., 2022). Unfortunately, open-forest data are not available for historical timeframes, precluding further tests of this hypothesis. Another possibility concerns habitat use by L. vulgaris, particularly the lesser importance F I G U R E 3 Plots investigating time lags in genetic differentiation. Isolation by resistance model fit (delta AIC; y-axis) per timeframe (1963-2015, every 4 years; x-axis) and IBR model transformation (depicted by shape of points; and resistance code: cf = closed forest, of = open forest, c = crops, u = urban, f = forest, nf = no forest). Plotted independently by evolutionary unit (Lm-L. montandoni; Lva-L. v. ampelensis; and Lvv-L. v. vulgaris). Colour by species, L. montandoni in red and L. vulgaris in green. Lines connecting dots were used to facilitate the visualization of temporal trends. Only IBR models better than IBD are shown, and IBD model fit displayed with dashed line. None of the historical forest versus non-forest IBR models explained genetic differentiation in Lva or Lvv better than IBD and are therefore not shown. Note different scales in delta AIC. AIC, Akaike information criterion; IBD, isolation by distance; IBR, isolation by resistance.
of optimal habitat matrix for connectivity due to its generalist ecology (Antunes et al., 2022). Indeed, L. vulgaris has often described the most generalist urodele species occurring in the Western Palearctic, broadly distributed and occurring across a wide variety of natural and human-made habitats at lower elevations (Speybroeck et al., 2016). In contrast, its sister species, L. montandoni generally occurs in more pristine forested habitats at higher elevations. Many studies have documented the habitat requirements and movement patterns of L. vulgaris, and although smooth newts were seen in ponds across a wide variety of habitats, proximity to forest edges or other structural complex habitats (e.g. scrubs or woodlands) was identified as a key factor that offered clear benefits for foraging, shelter and hibernation (Marnell, 1998;Mulkeen et al., 2017;Rannap et al., 2012;Skei et al., 2006). A survey of L. vulgaris movement patterns revealed that most individuals indeed prefer to move in forests, but still up to ~20% of migrants were moving in open areas (e.g. pastures and grasslands; Müllner, 2001). If the resistance to movement is similar between favourable habitat (e.g. open forest, scrubs, woodlands) and nonhabitat (e.g. croplands and urban areas), our ability to detect landscape genetic relationships is limited (Cushman et al., 2013). Thus, the generalist ecology of L.
vulgaris may lead to connectivity patterns that are better predicted mostly by distance and fine-scale site habitat characteristics (e.g. pond availability, lower urbanization; Moor et al., 2022). In fact, the distance was the main factor explaining genetic differentiation in L. vulgaris, with IBR models showing only slightly better fit, particularly for L. v. vulgaris (Figure 3; Figure S6).
Fine-scale patterns of IBD indicate that the dispersal of L. vulgaris is more restricted than that of L. montandoni across our study area, with low levels of genetic differentiation absent for L. v. vulgaris at distances greater than 5 km and L. v. ampelensis at distances greater than 2 km (Figure 2; Figure S5). The absence of IBD at short distances (0-1 km; Table S4) suggests that dispersal is common enough to override drift and landscape barriers at this scale, in both species, but particularly in L. montandoni, which shows a very flat IBD trend associated with very low Fst values, in contrast to the erratic behaviour in L. vulgaris, especially L. v. ampelensis ( Figure S5). These patterns suggest that at very small scales the mobility of L. vulgaris is similar to that of L. montandoni, making the interspecific differences in mobility an unlikely explanation of the contrasting time lag signatures. However, studies directly investigating the dispersal ecology of L. montandoni are needed to gain further insight into species dispersal abilities and their landscape use. At distances greater than 1 km, where dispersal events become rare due to species dispersal abilities (Bell, 1977;Kovar et al., 2009;Schmidt et al., 2006), IBD emerges in L. montandoni (across all distance ranges, 0-2.5, 0-5 and 1-50 km), with its trends being generally milder compared with L.
vulgaris. These patterns suggest patterns of genetic differentiation at larger scales (1-50 km) reflect dispersal at fine scales, probably by influencing stepping-stone dispersal. Compared with L. montandoni, the smaller effective population sizes documented for a subset of L. vulgaris populations in our study area (Antunes et al., 2022) may also play a role by increasing drift and reducing time lags (Epps & Keyghobadi, 2015). Increased drift also limits the detectability of associations between population genetic differentiation and connectivity (Antunes et al., 2022;Winiarski et al., 2020). Overall, the low predictive power of the IBR model, stronger population isolation (i.e. stronger IBD, particularly at lower distances) and lower local effective population sizes observed for L. vulgaris, may be the outcome of its greater exposure to habitat loss (e.g. through urbanization and loss of ponds) rather than its dispersal behaviour. Previous landscape genetics work has suggested that despite the comparatively lower sensitivity of L. vulgaris to anthropogenic landscape change, it is still more affected than L. montandoni, which inhabits more pristine habitats (Antunes et al., 2022). Whether the lack of time lag in L.
vulgaris is driven by a faster approach of genetic differentiation to equilibrium or by its generalist ecology and landscape use remains to be tested.
Evidence for time lags can also be distorted by methodological choices known to have profound impacts on landscape genetic inference (Anderson et al., 2010;Peterman & Pope, 2021). One important choice is the raster cell size and thematic resolution (i.e. number of land-cover categories) used as input, which in our case was limited to ~1 km raster cells differentiating between forest/non-forest in historical timeframes, due to lower confidence in other land-cover variables across our study area for historical times ( Figure S2). Models including other land-cover variables and finer resolution may improve resistance inference and change time lag signatures (Anderson et al., 2010;. In this work, we note that, in L. montandoni, forest-based models provided better predictions than more complex models (i.e. those including cropland and urban areas) for the present. Thus, in the case of L. montandoni, increasing thematic resolution is likely to lead to similar IBR models and time lag signatures. In L. vulgaris better thematic resolution, particularly including open forest in historical models, has the potential to reveal hidden time lags. Regarding cell size, we emphasize that our approach captures connectivity above the ~1 km raster scale, which in small-bodied newts includes the occasional movement of migrants that may connect populations separated by distances greater than 1 km (Jehle & Sinsch, 2007;Schmidt et al., 2006;Alex Smith & Green, 2005) but also stepping-stone connectivity, which happens at a wide range of distances (Saura et al., 2014). Such resolution can of course obscure the effect of smaller landscape elements (e.g. ponds) or linear features such as roads, highways or rivers, which can act as barriers to dispersal (Figueiredo-Vázquez et al., 2021;Holderegger & Di Giulio, 2010). However, small landscape elements are not expected to obscure or confound the effect of larger landscape features across the many pairs of populations studied here. Another important choice is the analytical framework Peterman & Pope, 2021;Shirk et al., 2021). Ideally, studies should explore multiple approaches, while considering their strengths and weaknesses (Peterman & Pope, 2021;Shirk et al., 2021). Particularly approaches that allow fast optimization, with different models, and incorporating multiple variables and interactions among them (e.g. radish , Peterman & Pope, 2021). However, such approaches take as input all pairs of populations, which in our case leads to the inclusion of comparisons that fell outside our spatial scale of interest. Considering this, we optimized landscape resistance using a 'ResistanceGA' based framework (see Section 2; Peterman, 2018) and included only populations from 1 to 50 km apart. Future studies exploring the influence of methodological choices (e.g. analytical frameworks and spatial data) are, nonetheless, necessary to consolidate inferences about species land-use and population connectivity Peterman & Pope, 2021).

| Conclusions with implications for conservation and landscape genetics
One obvious conservation implication of time lags deals with the effects of human-induced habitat conversion on population differentiation. Our results for L. montandoni, corroborating previous work (e.g. Holzhauer et al., 2006;Maigret et al., 2020;Spear & Storfer, 2008), imply that it takes time for the genetic effects of increased connectivity to manifest themselves in some species.
Our results show a time lag of ~10 generations in L. montandoni, suggesting that the impact of forest-cover changes over the last 40 years on genetic differentiation has not yet fully manifested, particularly the increased connectivity across the Podkarpackie region. Considering this we can predict genetic differentiation between Podkarpackie populations to decrease in the future.
Conversely, the emergence of barriers to dispersal may also be subject to a time lag effect. For L. montandoni, the lower permeability of current non-forested habitats relative to historical times, due to increasing urban areas, road infrastructure/traffic and croplands, may decrease connectivity and ultimately outweigh the influence of increasing forest cover. Thus, calling judgement on the consequences of habitat disturbance (e.g. clearing of forest, construction of motorways, etc) on genetic differentiation without factoring in time lags may confound the interpretation of population connectivity estimates and prove detrimental to management decisions. Notably, models built across all timeframes always revealed the same forest-cover use in L. montandoni, showing that landscape genetic inferences can be stable despite landuse changes. However, we note that this may not be the case in scenarios with more drastic land-cover changes or when optimizations explore more resistance scenarios (i.e. higher number of variables and/or transformations). The incorporation of historical landscape changes into studies dealing with current patterns of genetic differentiation can be done in multiple ways, for example, including genetic metrics with different properties (Coster et al., 2015), using key historical environmental data (Maigret et al., 2020) or simply by considering time lags during discussion (McCluskey et al., 2022). This is particularly important for studies investigating the impact of habitat fragmentation caused by recent anthropogenic elements in the landscape, as in these cases the genetic consequences of decreased population connectivity may not have yet emerged.
Our results emphasize that time lag effects, and more broadly, measures of population differentiation reflecting population connectivity, are species-and context-specific, being distinctive even for closely related species or evolutionary units. In L. vulgaris we found no evidence for time lags, with current genetic differentiation predicted mostly by distance, but with evidence for higher connectivity across open forests, especially in L. v. ampelensis. This result may reflect a faster approach of genetic differentiation to equilibrium but can also be driven by the generalist ecology and landscape use of this species. The contrasting findings for L. montandoni and L. vulgaris highlight the variation in time lag prevalence in nature, showing that even when species face land-cover changes simultaneously, responses measured by genetic differentiation may differ.
Moreover, our results for the sister species L. montandoni and L.
vulgaris reinforce that anticipating responses to land-cover changes based on inference from similar or even closely related species should be avoided (Antunes et al., 2022). In summary, our results reveal some consistency in landscape genetic inferences made from current patterns of genetic differentiation but ultimately highlight the interspecific variation in time lag prevalence, indicating that conclusions about species landscape use and current habitat connectivity should be carefully interpreted in the context of historical landscape changes.

ACK N O WLE D G E M ENTS
We thank Weronika Antoł, Aleksandra Błażejowska, Maciej Initiative -Research University' at the Faculty of Biology of the Jagiellonian University in Kraków, Poland. We thank three anonymous reviewers for their valuable comments that helped to improve this work.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ ddi.13748.

DATA AVA I L A B I L I T Y S TAT E M E N T
The main workflow and data used in this project were shared via Dryad Digital Repository; Link: https://doi.org/10.5061/dryad. ht76h drkz. Raw sequencing data (FASTQ) have been deposited in the NCBI sequence read archive (PRJNA975893).

B I OS K E TCH
Bernardo Antunes is a PhD student at the Institute of Environmental Sciences, Jagiellonian University, Kraków, Poland.
His research aims at better understanding how environments shape genetic diversity, and how this knowledge can be used to