Genetic signature of the natural gene pool of Tilia cordata Mill. in Lithuania: Compound evolutionary and anthropogenic effects

Abstract Tilia cordata Mill. is a valuable tree species enriching the ecological values of the coniferous‐dominated boreal forests in Europe. Following the historical decline, spreading of Tilia sp. is challenged by the elevated inbreeding and habitat fragmentation. We studied the geographical distribution of genetic diversity of Tilia cordata populations in Lithuania. We used 14 genomic microsatellite markers to genotype 543 individuals from 23 wild‐growing populations. We found that Tilia cordata retained high levels of genetic diversity (population F IS = 0–0.15, H o = 0.53–0.69, H e = 0.56–0.75). AMOVA, Bayesian clustering, and Monmonier's barrier detection indicate weak but significant differentiation among the populations (F ST = 0.037***) into geographically interpretable clusters of (a) western Lithuania with high genetic heterogeneity but low genetic diversity, bottleneck effects, (b) relatively higher genetic diversity of Tilia cordata on rich and most soils of midland lowland, and (c) the most differentiated populations on poor soils of the coolest northeastern highland possessing the highest rare allele frequency but elevated inbreeding and bottleneck effects. Weak genetic differentiation among the Tilia cordata populations in Lithuania implies common ancestry, absence of strong adaptive gradients, and effective genetic exchange possible mediated via the riparian networks. A hypothesis on riparian networks as gene flow mediators in Tilia cordata was raised based on results of this study.

historical gene pool of Tilia cordata in the retreats of the northerly forests? What is the genetic structure within a tree species that is largely insect-pollinated and forecasted to form mixed forests in the future? Is it structured in small isolated pollutions or, on the opposite, as large intermating gene pools as wind-pollinated species do.
Revealing these genetic properties will improve understanding of species ecology and support the natural integration of Tilia cordata into the northerly forests (Ennos, 2003;Myking, 2002).
Tilia cordata is an interesting species for the population genetic studies because of its turbulent evolutionary history, insectmediated reproductive system, and extensive urban use (De Jaegere et al., 2016). Below, we extracted important details on the effects of these three factors on the gene pool of Tilia sp. The warm Atlantic period (ca. 7000-4000 BC) is considered as the golden age for Tilia sp. dominating the forests of northern Europe (Rehfeldt, 1994).
Following the cooling temperatures at the later part of Holocene ca. 3000 years BC, frequency of Tilia cordata dropped markedly in the northerly forests (Balakauskas, 2012;Hewitt, 2004). In the absence of favorable temperatures at the right time for seed germination and pollination/fertilization, success leads to scattered retreats of Tilia cordata into fragmented forest gene pools (Pigott & Huntley, 1981). Owing to a reduced dispersal capacity in a cooler climate, Tilia cordata rarely outcompetes the rival pioneers and seldom reoccupies the former habitats (Pigott & Huntley, 1981). In the later centuries, forest clearance for agriculture has further reduced the remaining autochthonous Tilia forests in Europe. Provided this retreat of the natural gene pools, the Tilia plantations in avenues "under the linden tree" of urban parks are likely to have served as the escape sources of Tilia groups into the wild already since the medieval times (Lefevre, 2004). Therefore, we may deal with compound genetic structures for such species as Tilia: the autochthonous gene pools and introduced populations as the case of European beech in Lithuania (Kembryte et al., 2021). The surviving ancient especially coppice-born admixture of Tilia cordata in European forests is often considered as an indicator of the autochthonous woodlands (Hermy et al., 1999;Lawesson, 2004;Rackham, 2008).
At the ecosystem level, Tilia cordata forms biocommunities with added benefits for ecosystem services, including bee conservation (Anderson, 1976;Free, 1970;Pigott & Huntley, 1981). In the forest stands, presence of an admixture of Tilia cordata is associated with high diversity of plant species (Normann et al., 2016) and its decomposing leaves enrich soil with nutrients (Muys et al., 1992).
According to Barker (2017) and Logan et al. (2015), ca. 25%-40% of Tilia cordata trees within UK populations originate from vegetative sprouting. In the generative system, the flower nectar production and pollen germination require temperatures above 15°C (Pigott et al., 1980;Tal, 2006). Tilia certainly is an entomophilous species, pollinated by mainly bees (daytime) and short tongue moth (nighttime) (Anderson, 1976;Pigott, 1991). Wind is a secondary pollinator in Tilia, however, due to the dense crowns with fully extended, the effective wind-mediated pollen flow reaches fewer than ca. 200 m (Anderson, 1976). Whereas bees may fly up to 6 km away from their nests (Pasquet et al., 2008). Dispersal of mature Tilia seed usually is limited to ca. 300 m (Pigott, 2012). The bracted seed can float and be viable for germination after a 30 km travel along the steams (Pigott, 2012). This riparian movement of seeds may serve as a potential gene flow mediator.
Despite the sensitive reproduction system, and environmental and anthropogenic pressure, Tilia cordata populations survived over F I G U R E 1 Location map of the populations studied (left). The three main ecoclimatic regions (1, 2, and 3) are separated by the solid line (mainly coastal-continental gradient affected by the altitudinal variation) and the subregions (a, b) reflecting a finer zoning are marked by dashed lines. Population position markers are colored according to the three main regions. Altitudinal gradients in Lithuania (left, the highlands peak to 240-290 m a.s.l.). The northeastern highland (region 3) contains the coolest sites in the country F I G U R E 2 Interpolated surface plots of the within-population genetic diversity parameters revealing the geographical patterns of genetic diversity of Tilia cordata in Lithuania. The shape of the plots is a simplified map of Lithuania. The actual population values are given in Table 3 and Figure S1. The lowest right surface plot illustrates geographical distribution Tilia cordata-dominated sites (the statistics used on the map is the total area (in ha) of sites with >50% of Tilia cordata within a forest management unit of ca. 40 th. ha in size (forest inventory data from 2010)). N ef -effective population size (the MIGRATE_N estimate), H o and H e -observed and expected heterozygosity, and F IS -inbreeding coefficient. Rare alleles 5% is population number of alleles below the 5% frequency vast territories from Spain to central Sweden in the north and the Ural Mountains in the east (Pigott, 2012). The northern expansion of Tilia cordata was limited by the lack of favorable summer temperatures required for successful reproduction (Pigott & Huntley, 1981).
In different parts of Europe, these autochthonous gene pools were subjected to diverse adaptive pressure, land-use history, urban development, and forest use intensity. This led to a marked regional variation in the genetic condition of Tilia cordata populations in Europe (Erichsen et al., 2019;Lobo et al., 2018;Logan et al., 2019).
For instance, the competitive advantage of European beech is usually referred as an important cause for decline of Tilia sp. in Europe (Turner, 1962). However, in most of the present-day semiboreal and boreal zones, beech was not present.
As regards the eastern Baltic region, the forests are in the transitional zone between the boreal and temperate forests. However, recent climatic shift leads to a steady expansion of mixed broadleaved forests (Pearson & Dawson, 2003). The present-day forest cover in Lithuania is 33.5%. A large share of these forests is autochthonous, naturally regenerated over centuries, and fits well to the definition of ancient woodland (Peterken, 1974). Eco-climatically, the country is split into coastal lowland, western highland, midland lowland, and eastern highland ( Figure 1). Scots pine (Pinus sylvestris) dominates on poor sandy sites in south and southeast. Mixed Norway spruce (Picea abies) Scots pine forests are common in the northeastern and western highlands. Midland lowland contains relatively richer soils and is dominated by broadleaved tree species. According to the forest inventory from 2010, the forest sites with over 50% dominance of Tilia cordata comprise 6,500 ha in Lithuania, which is 0.32% of the total area of forest land in Lithuanian. Tilia cordata is more frequent on rich soils in the forests of midland lowland, especially its southwestern and northern parts (Semaškienė, 2006, Figure 2 (Semaškienė, 2006). Tilia cordata usually forms rather fragmented communities of ca. 5-10 ha at large. These communities can be found deep inside forest tracts as well as at the edges. Tilia cordata is also very common tree for landscape amenities and urban areas, such as city avenues, and modern and old manor house parks of over entire eastern Baltic (Tauras, 1989).
Especially after the devastation of WW2, there was a campaign for green landscape restoration that encouraged establishment of small parks and planting trees along the streets in small villages, towns, and cities in Lithuania (Tauras, 1989). In many of these post-WW2 plantations, Tilia cordata was used as a native entomophilous tree with cultural and medical values. Likely seed sources for these trees were old manor house parks (Tauras, 1989). Tilia platyphyllos does not occur naturally in Lithuania but is not rare as a decorative tree in urban parks.
Population fragmentation in Tilia cordata often leads to depleted diversity due to genetic drift, especially hazardous for entomophilous species with limited gene flow capacity (Erichsen et al., 2019;Lobo et al., 2018;Lowe et al., 2015). In addition to that, spreading from Tilia groups of likely low diversity in urban territories may further compromise the genetic stability of this newly expanding Tilia cordata populations. The results from the few available microsatellite-based studies on Tilia cordata are diverse. Logan et al. (2015) found no marked reduction of genetic diversity of Tilia cordata in the UK, located not far from the margin of the range, whereas Lobo et al. (2018) reported possible genetic drift effects on Tilia cordata in Demark. Therefore, the situation with the genetic diversity of Tilia cordata may vary depending on the region. The above-mentioned studies investigated Tilia cordata populations in rather fragmented forests surrounded by urbanized territories in western Europe. It would be interesting to obtain the genetic diversity and structure estimates for Tilia cordata from forested regions of Europe such as eastern Baltic forests.
The objective of our study was to assess the geographical distribution of the genetic diversity and the genetic structure of Tilia cordata populations in Lithuanian and discuss the main factors affecting the gene pool of Tilia cordata as well as its genetic potential for further spreading into mixed forest stands. For this purpose, we used a set of 14 genomic microsatellite markers to genotype a representative network of natural populations of Tilia cordata in Lithuania.

| Sampling sites
A total of 543 Tilia cordata individuals were sampled from 23 wildgrowing populations evenly covering the territory of Lithuania ( Figure 1, Table 1). We carefully examined the chosen trees for Tilia cordata morphology to avoid unlikely but possible hybrids with Tilia platyphyllos, which is an exotic tree in Lithuania. We have not observed any potentially T. platyphyllos like individuals in the sampled stands (large hairy leaves with thick cuticula, angled fruits, multiple forking of stems due to frost damage). However, the northern PAPI population contained features less common elsewhere: pale bark with less expressed farrows, pronounced multiple forking of stems. The Tilia cordata sites were carefully chosen for natural origin and location by avoiding proximity to the urban areas, usually within forest tracts of various size. We used an electric drill to sample the sawdust from 20 to 25 adult overstory Tilia cordata trees per population. We selected the trees spaced at least 20 meters away from each other within a zigzag sampling path covering ca. 2-5 ha depending on the density of Tilia cordata at a site. In most of the stands we sampled, there were sites with more than 50% dominance of Tilia cordata.

| DNA analysis
The DNA was extracted from silica gel-dried wood sawdust collected by drilling with electric bore ca. 1 cm deep into the trunk (e.g., Verbylaitė et al., 2010, bark discarded, bore diameter 0.5 cm). A modified CTAB Doyle and Doyle (1990) DNA extraction method was used to extract DNA. For the genotyping, we used 14 genomic microsatellite DNA markers (Phuekvilai & Wolff, 2013, Table S1). The PCR amplification was carried on three multiplexes in a 15 µl reaction mix containing 7.5 µl of a PCR Master Mix, 3 µl of RNAse free water, 1.5 µl of Primer Mix, 1 µl of DNA, 1 µl of PVP 1 (%), and 1 µl of BSA 20 mg/ml (bovine serum albumin) (Applied Biosystems Thermo Cycler GeneAmp PCA System 9700) with following PCR thermocycling profile: initial denaturation step at 95°C for 15 min, followed by 25 cycles each of 94°C for 30 s, annealing temperature at 54°C for 1 min, 30 s, and extension at 72°C for 30 s, followed by the final extension step at 60°C for 30 min. Fragments were separated with the capillary electrophoresis on ABI PRISM™ 310 genetic analyzer. The alleles were scored on GENEMAPPER soft. ver. 4.1 (APPLIED BIOSYSTEMS). We constructed the allele scoring binset based on Phuekvilai and Wolff (2013,

| Data analysis
We assessed the frequency of null alleles for each locus with all populations pooled with MICROCHECKER soft. ver. 2.2.3 (Van Oosterhout et al., 2006). Occurrence of clones and the standard genetic diversity parameters along with the Mantel test statistics for isolation by distance were calculated with GENALEX soft. ver. 6.5 (Peakall & Smouse, 2006). The rarefied allelic richness (Ar, base 17 chosen for the lowest sample size in population PAGR, Table 1) and the inbreeding coefficient (Fis along with the significance of its deviation from 0) were calculated with FSTAT soft. ver. 2.9.3.2 (Goudet, 1995). We tested the significance of differences between the regions in main genetic diversity parameters by using the FSTAT among group significance tests based on 1,000 permutations.
The effective population size (N ef ) was calculated for each population and region based on maximum-likelihood method (Hastings-Metropolis Markov chain Monte Carlo algorithm) and coalescent theory using MIGRATE_N software (Beerli, 2009). This program returns the theta (ϴ) value, which was used to calculate the effective population size N ef = ϴ/μ, where μ is assumed mean microsatellite mutation rate per generation of 4 × 10 −3 (Boys et al., 2005;Pandey & Rajora, 2012). We also use MIGRATE_N to calculate the number of inward and outward migrants per generation at the population Total 543 Note: Region indicates climatic zones in Lithuania mainly affected by continentality and topography as shown in Figure  We run the Wilcoxon rank test for heterozygosity excess and the mode shift test to screen for bottleneck effects on the region level with software BOTTLENECK 1.2.02 (Cornuet & Luikart, 1997). The Wilcoxon rank test (Luikart, 1997) assumes that in a population after a recent bottleneck, the expected heterozygosity is decreasing more slowly than the allele dropout.
Therefore, the expected heterozygosity (H e ) observed in the data (or homogeneity of allelic frequencies) is higher than H e expected under mutation-drift equilibrium (Garza & Williamson, 2001). The software tests for significant number of loci that have H e excess.
This assumption suites better for the loci under the IAM (infinite allele model) than SMM (stepwise mutation model) (Cornuet & Luikart, 1997). However, we tested all three mutation models: SMM, IAM, and two-phase mutation (TPM). In the Wilcoxon rank test statistics, we checked the one-tail p-values for H e excess and examined the homogeneity of the one-tail p-values for the H e excess and deficit (because in a population at mutation-drift equilibrium; there is approximately an equal probability that a locus shows H e excess or deficit). The mode shift tests for deviations from so-called L shape distribution of allele frequencies usually indicate recent bottlenecks. This test assumes that the bottleneckfree populations possess high number of low-frequency alleles which follow an L shape appearance of the allele frequency histogram (Luikart, 1997).
The interpolated surface plots of genetic diversity parameters were obtained by calculating the matrix of interpolated values with the MS EXCEL function =@Interp2d and then drawing the matrix rainbow plot by the id of the POPTOOLS EXCEL add-in.
The population differentiation was tested by calculating (a) the frequency-based differentiation indexes G st (an analogue of F st adjusted to variable sample size) and D est (Jost, 2008) with GENALEX soft. ver. 6.5 and (b) running an AMOVA to partition the molecular variance among regions, populations within regions, and within populations with Arlequin soft. ver. 3.5.1.3 (Exoffier & Lischer, 2010).
We used the Bayesian clustering approach implemented in STRUCTURE soft. ver. 2.2.3 (Pritchard et al., 2000) to investigate the population genetic structure of Tilia cordata in Lithuania. We set the burn-in period length for posterior distribution to 10 5 and the number of MCMC iterations to 10 5 , the K range from 1 to 10, each replicated 10 times. We used the correlated allele frequency model, no admixture, and the LOCPRIOR option for the three regions indicated in Table 1 and Figure 1. The most likely number of genetic clusters K was identified based on the deltaK value (Evanno et al., 2005) with STRUCTURE_HARVESTER WEB soft. (Earl & von Holdt, 2012). In genetically heterogeneous material especially under presence of botanical hierarchy (e.g., hybridization or transfer effects), the optimal number of genetic groups may be underestimated (e.g., Wang, 2017). Therefore, we presented the genetic structures for a range of K values from 2 to 6. We did not consider deploying higher than K = 6 groups because the individual Bayesian assignment We used Monmonier's algorithm allowing for establishing barriers along a significant shift in the allele frequency within a landscape implemented in soft. BARRIER ver. 2.2 (Manni et al., 2004). The program (a) creates a Delaunay triangulation plot between the sampled populations, (b) calculates genetic distances (Nei et al., 1983) associated with each edge in the plot, and (c) creates growing barriers along the largest genetic distances on the plot; the barriers are ranked based on the magnitude of the differentiation. The program also allows for a significance test of the barriers by analyzing the bootstrapped distance matrices and displaying the number of Note: N obs is the sample size; N a is number of different alleles. N e is effective allele number, H o and H e are observed and expected heterozygosity, A r17 is rarified allelic richness with the base of 17 individuals, F is is the FSTAT inbreeding coefficient and its significance (pF ST , 5,000 permutations), N ef is effective population size (coalescence in MIGRATE_N), and Rare5% is number of rare alleles below the 5% frequency. The regional means of the genetic diversity parameters with the same latter are not significantly different according to the pairwise FSTAT significant tests between region means at 0.05 significance level of the two-sided p-value based on the 1,000 permutation. SE is standard error calculated from the population/locus mean values (except N ef from population mean values).
Finally, we used the ALLELES IN SPACE soft. ver. 1.0 (Miller, 2005) to obtain a surface plot of the interpolated Nei et al. (1983) genetic distances among Tilia cordata populations in Lithuania.
We also screened for possible natural hybridization events between T. cordata and T. platyphyllos by examining the geographical distribution of the less common alleles at the two loci Tc8 and Tc927. These loci were reported to discriminate well between these two Tilia species (Logan et al., 2015). The authors reported that the locus Tc8 was monomorphic in the UK with T. cordata material (the 96% frequent 140 bp allele in our material) but returned 8 alleles in T. platyphyllos and the locus Tc927 few alleles in T. cordata and 15 alleles in T. platyphyllos.

| Loci statistics
We detected no clones among the 543 trees of Tilia cordata in our The loci TC8 and TC11 were least informative, monomorphic in most of the populations, possessed high frequency of null alleles (0.17 and 0.16, respectively), and therefore were excluded from further analysis of genetic structure and diversity ( Table 2). All the other loci had null allele frequencies below 0.1 except locus TC31 (freq. null = 0.17, Table 2). The locus TC31 was highly informative and the inbreeding coefficient at this locus varied markedly among populations (not shown), implying that null alleles at this frequency of 0.17 do not markedly deviate the F is values.

| Within-population genetic diversity
Within-population genetic diversity varied markedly among the populations. The allelic diversity (allelic richness, A r and expected heterozygosity, H e ) in the western region 1 was significantly lower than in the remaining regions ( Figure 2, Table 3, Figure S1). In western Lithuania, 5 of 8 populations possessed the lowest ranks of H e from 0.56 to 0.65, including the two adjacent northeastern populations of PAPI and ZAGA (Table 3, Figure 2, Figure S1). The PAPI and ZAGA populations also showed a high degree of differentiation from the rest in his region (subsection below and Figure 4). There were no significant regional differences in the observed heterozygosity (H o ) and inbreeding coefficient (F IS ) (  Table 3 and Figure 2). Elevated F IS values were also observed in the eastern highland, where climate is cooler, soils are poor and dry but Tilia is least common (region 3, Figure 2).
Geographical distribution of the effective population size N ef followed that of H e with the lowest values in western highland and highest-in eastern highland ( Figure 2, Table 3, Figure S1). Finally, rare allele numbers peaked in the eastern highland (region 3) and dropped to the lowest ranks in coastal lowland (region 1B, Table 3,

| Genetic differentiation and structure
The frequency-based differentiation tests at all loci retuned sig- (0.13), Tc963 (0.25), Tc7 (0.14). These five loci also were the most polymorphic with highest ranking N a and H e values ( Table 2). The  (Table S2).
The Mantel test showed weak but significant linear correlation between the genetic and geographical distances (R xy = 0.047, permuted p-value = .001), indicating some isolation by distance and substructuring among Tilia cordata populations in Lithuania.
Evanno's deltaK method suggested K = 4 as the optimal number of the Bayesian clusters based on the 10 replicated STRUCTURE runs for each K from 1 to 10 ( Figure S2). However, by considering the observation that the STRUCTURE_HARVESTER algorithm tends to underestimate deltaK in a genetically heterogeneous material (Hintsteiner et al., 2018) and seeking for a finer structure, where several rare clusters occurred at high frequency and several geographically close populations were assigned to different genetic clusters (e.g., PAPI-ZAGA, and SVEK1-SVEK2 in Figure 4); (b) a rather genetically homogeneous group of populations (region 2), free of the genetic outliers as in western and eastern Lithuania, and (c) the most genetically differentiated population group located in the eastern highland sharing elsewhere rare cluster in relatively higher proportions (e.g., cluster id. 6 reaching 50% in the KAVA population in Figure 4).
The additional STRUCTURE runs to screen for a finer structure by analyzing a subset of populations revealed (a) a subdivision of the western region 1 into northern 1A and southern 1B, in this way separating the former Eastern Prussian territory from the rest (b) slight differentiation of the two northern populations of PAPI and ZAGA from the rest in the midland highland ( Figure S5).
The Monmonier's barriers along marked shifts in allele frequencies were not as helpful in finding geographical consistent genetic structures ( Figure S3). The barriers were drawn around most of the populations, especially in the eastern and western regions ( Figure S3). Therefore, a more meaningful interpretation of this analysis is identifying areas where the barriers are not identified given maximum reasonable occasions for identifying the barriers.
Consequently, under the settings to identify 10 the most significant barriers given 100 bootstrapped matrices, the midland lowland was the only large-scale barrier-free region ( Figure S3).
The surface plot of the interpolated Nei et al. (1983) genetic distances basically differentiated eastern the eastern highland from the remaining regions, which confirms the structure revealed by the Bayesian clustering ( Figure S3).
As regards possible hybridization between T. cordata and T. platyphyllos, we examined the putative T. platyphyllos marker allele frequencies for the loci Tc8 and Tc927. For the locus Tc8 with the most common 140 and 146 bp alleles at 99% frequency (Table 2)

| Genetic features of Tilia cordata in Lithuania
Our study revealed that despite the speculated threats for frag- Denmark (Erichsen et al., 2019). The genetic diversity levels were even comparable with wind-pollinated species such as Scots pine in Lithuania (Danusevičius et al., 2016). The inbreeding estimates in our study were lower than in the Danish material of Tilia cordata, where more asexual reproduction exists (Erichsen et al., 2019).
To preserve such high levels of genetic diversity, an evolving network of natural populations is required (De Jaegere et al., 2016).
Therefore, it is likely that autochthonous gene pool of Tilia cordata must have been maintained by synergy of gene flow, natural selection, and genetic drift in the network of forest tracts in Lithuania.
Assuming a similar landscape management history and ecoclimatic conditions in the Baltic region and southwestern parts of European Russia, our results may be generalized onto a broader geographical scale. In 1700s and 1800s, the forest cover of Lithuania was 60% and 50%, respectively, followed by a gradual reduction to ca. 20% due to forest feelings during WW2 (Kairiūkštis, 2003). Then, presumably, Tilia cordata experienced further reduction in population size surviving in the deep of large woodlands. It is likely that these ancient woodlands especially on moist and rich soils in the lowlands were of sufficient ecological capacity to sustain a sound network of Tilia cordata populations. In contrast, rural areas of the UK and Denmark are more urbanized likely leading to stronger imprints of low diversity urban escapes on natural gene pool of Tilia cordata.
Another explanation for the relatively higher genetic diversity of Lithuanian Tilia cordata populations could be lower rates of asexual reproduction, which is sensitive to genetic drift effects (Balloux et al., 2003). The fact that we did not detect any clones in the 543 Based on the Tc8 marker (Logan et al., 2015), the hybridization between T. cordata and T. platyphyllos is very low (detected at low rare in a single coastal population of JUOD  Figure S4), but these were natural populations with no direct evidence of mixture with urban gene pools.
In agreement with the earlier studies on Tilia sp. with the SSR markers (Barker, 2017;Erichsen et al., 2019;Lobo et al., 2018;Logan et al., 2015), we did not detect more than 2 alleles at all loci, confirming the diploidy of Tilia cordata. In our study, the observed heterozygosity estimates are not markedly affected by null alleles as shown by the MICROCHECKER tests and the high variation in the singe locus H o and F IS values among the populations.

| Within-population diversity is geographically variable in Lithuania
We clearly observed low values of allelic diversity and tendency of For post-WW2 urban plantations, the original sources of Tilia cordata are unknown, but likely to be old parks and seldom autochthonous forest sources (Tauras, 1989). Importantly, the old manor park designers over centuries often used nonlocal collections and botanical gardens as the seed sources to surprise their customers (Rosłon-Szeryńska et al., 2020). Therefore, all this urban source spreading of Tilia cordata is likely to contain nonlocal sources of unknown adaptedness in the wild.
In Lithuania, Tilia cordata has its highest frequency in the western part (Figure 2). How could it be that these most Tilia-rich areas show lowest genetic diversity? A possible answer is significant effects of the urban gene pools, and the urban escapes may not be recent. Western Lithuania has a high concentration of old manor parks (Tauras, 1989). Good candidates for such urban escapes are the two adjacent PAPI and ZAGA populations (northwest): lowest allelic diversity and heterozygosity values (Table 3, Figure S1), or the SVEK2 population in the former Eastern Prussia with the distinct morphology and the genetic identity outstanding from the adjacent SVEK1 and the remaining populations in this region (Figure 4). Use of nonlocal material for forestry was not uncommon in the former Eastern Prussia .
In the eastern Lithuanian highlands (region 3), the situation with genetic diversity was different-high allelic diversity but elevated inbreeding and bottleneck effects detected ( . There also is a relatively higher urbanization pressure in eastern region 3 ( Figure S5) that may lead to formation of internally mating groups we well. Another possible reason could be more recent natural spreading of the western Russian Tilia cordata gene pools because of warming climate (supported by a higher rare allele frequency in the eastern region 3, Figure 2; McLachlan et al., 2005).
In the long run, these eastern immigrants may help to recover after the bottleneck effect.
Moist and rich sites such as the midland lowland (region 2A) create favorable environment for Tilia cordata (e.g., Barker, 2017) allowing for larger and less fragmented populations in midland lowland ( Figure S5). This may explain relatively lower inbreeding values for Tilia cordata populations in midland lowland.

| Genetic structure is associated with geography in Lithuania
Our results on genetic structure of Tilia cordata in Lithuania contrast with Erichsen et al. (2019) study where no clear population structuring was found for Tilia cordata in Demark with the set of SSR markers. Based on the findings discussed above, this genetic structure is likely to be shaped by the balance between drift caused fragmentation and gene flow as well as local effects of urban escapes (Lobo et al., 2018). There also could be signs of adaptation to poor soils and relatively higher frost tolerance in the coldest region in Lithuaniathe northeastern highland.
Theoretically, insect pollination implies weak gene flow among populations (e.g., Lowe et al., 2015). Hence, the gene flow rates obtained by the coalesce algorithms in our study may also reflect seed migration vectors mediated by wind, water, and birds (Calladine & Morrison, 2010;De Jaegere et al., 2016). The first observation on gene flow patterns is free mutually equal genetic exchange among the midland lowland regions 2A and 2B (confirmed by the BARRIER structure analysis, Figure 2 and Figure S3). An explanation for this could be seed transport via a riparian network of the largest river in Lithuania-Nemunas and its affluent Nevezis. A relatively stronger gene flow from the flanking highlands (regions 1 and 3) downwards into the midland lowland (region 2) also supports the riparian gene flow hypothesis ( Figure 5). A uniform adaptive environment of the midland lowland could contribute to homogeneous genetic structure in this region as well. The second observation on gene flow is on the outward flow from the bottlenecked regions being stronger than the inward flow ( Figure 5). This may indicate that spreading of Tilia cordata is more successful on the optimum sites such as the most and rich soils of midland lowland (region 2) than elsewhere. In our material, there also is a tendency for a stronger northward than southward directed gene flow ( Figure 5), supporting the observations on advance of Tilia northwards (Logan et al., 2019).
As regards the genetic structure of Tilia cordata in Lithuania, we found three regions showing a genetic structuring as follows: (1)  and microhabitat diversity helps forming gene pool diversity and so strengthening species for challenges of climate change (Hampe & Petit, 2005;Lobo et al., 2018).
In our study, the Bayesian clustering was more helpful in revealing these geographically consistent genetic structures than the barrier detection by Monmonier's algorithms. As discussed by Safner et al. (2011), the algorithms designed to detect gradients of significant change in allele frequencies are less efficient with highly heterogeneous material as they tend to delineate the outstanding populations rather than the large-scale geographical structures. In our material, however, the BARRIER runs were useful not in finding where the barriers are but where they are not present given the settings to find any possible barriers (large number of bootstraps and allowing many barriers, Figure S3). In agreement to the above, the weakest barriers in our material were found among rather homogeneous populations in the midland lowland (region 2A and 2B, Figure S3).
The Tilia cordata populations on the poor and dry soils of the eastern highland (region 3) were most differentiated in Lithuania.
There could be several reasons for that (a) adaptation to cooler climate, drier and poorer sites in this region, (b) gene flow from the east as discussed above on high rare allele frequency in this region, (c) stronger genetic drift effects due to relatively stronger fragmentation caused by to higher urban pressure and harsher environment for Tilia in this region.
Our findings indicate that escapes from urban sources may lead to a reduction of genetic diversity in Tilia cordata populations. Such diversity reduction is especially undesirable at the northern edge of species advancement, where fragmentation of Tilia cordata gene pool is still profound and the environment for sexual reproduction is less favorable (Logan et al., 2019;Pigott et al., 1980). Measures to reduce spreading of Tilia sp. from urban sources must be reviewed (like disposal of leaves in the forests). Conservation efforts could consider capturing the genetic diversity of Tilia cordata within each of the western, midland, and eastern genetic clusters in Lithuania by identifying autochthonous stands. Asexual reproduction, presence of old-growth trees, and variable age structures could be considered as the main criteria for the autochthonous populations. The commercial value of Tilia cordata trees can be one but not a decisive factor for selecting the conservation stands. Natural population management should favor preservation of old-growth trees of Tilia cordata in mixed forest stands ( Figure S6). These trees are the best shelter for bee nests and so successful pollination so needed to produce vital seeds and successful spreading of Tilia cordata.

ACK N OWLED G M ENTS
The study was supported by the Lithuanian Science Council project no. S-MIP-19-3 7.5. We also appreciate efforts of the reviewers for improving this manuscript.

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