The targeted sites were in Heilongjiang, Jilin, Liaoning, Beijing, Hebei, Shandong, Anhui and Zhejiang Provinces in China (four sites in Shandong and one site in each of the remaining provinces). A total of 242 individuals were collected with eight T. amurensis populations (TS, LU, LS, ZS, YDS, DGL, MPS, and DXG) and three T. japonica populations (HS, MS and TMS) (Table 1; Figure 1). The sampling work was carried out from May to July 2023. The selection of target plants strictly followed the principle of random sampling, with a minimum distance of 30 meters between any two sampled individuals (the target sample size for each population was 30 plants; if the actual number of individuals was fewer than 30, all were sampled). Fresh young leaves of linden trees were collected, identified by Prof. Dekui Zang, and preserved in silica gel for subsequent DNA extraction. Voucher specimens were also preserved in the herbarium of Shandong Agricultural University.

We used the SteadyPure Plant Genomic DNA Extraction Kit (CTAB method) for DNA extraction, and the experimental procedures were performed according to the methods of Wu’s study (Wu et al., 2024). To evaluate the quality and concentration of the extracted DNA, we used the NanoDrop™ One/OneC Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific, USA) for measurement. In accordance with the experimental requirements, the DNA solution was diluted to 30 ng·μl^-1 for subsequent PCR amplification experiments.

We obtained SSR primers reported for the study of Tilia species and selected 27 microsatellite primers with relatively simple repeat motifs and high polymorphism rates for screening. The primers were synthesized by Shanghai Personalbio Technology Co., Ltd. Fluorescent dyes were added to the 5’ end of each primer. The addition of a fluorescent label enhances the functionality, stability, and detection efficiency of the primers, facilitating the rapid identification of amplification products in subsequent high-efficiency capillary electrophoresis. A total of 15 SSR primers from Yue’s study were used for subsequent experiments (Supplementary Table S1) (Yue et al., 2022). The reaction system for amplification and the conditions and parameters for PCR were as described in our previous study (Wu et al., 2024). The PCR amplification products were subsequently detected using 1.5% agarose gel electrophoresis. The PCR products that met the experimental requirements were purified and sent to Shanghai Personalbio Technology Co., Ltd., for high-efficiency capillary electrophoresis analysis.

The genetic diversity was estimated on the basis of the SSR data and analyzed using GenAlExv6.51 software to obtain the observed allele number (Na), efficient allele number (Ne), Shannon’s information index (I), observed heterozygosity (H_O), expected heterozygosity (H_E) and inbreeding coefficient (F_IS) (Peakall and Smouse, 2006). The obtained data were subsequently subjected to AMOVA and PCoA via GenAlExv6.51 (Smouse and Peakall, 2012). The genetic differentiation among populations was calculated via the genetic differentiation coefficient (F_ST), and the gene flow between populations was determined using the F_ST value: N_m = (1/F_ST − 1)/4 (Weir and Cockerham, 1984). UPGMA clustering analysis was performed using PopGen32, and a phylogenetic tree was constructed via MEGA7 to determine the genetic relationships among populations (Kumar et al., 2016). We used the Bayesian clustering method in STRUCTUREv.2.3.4 to investigate population structure (Pritchard et al., 2000). The optimal genetic group number K was determined by setting the “length of burn-in period” to 1,000, the “number of MCMC reps after burn-in” to 20,000, and calculating the ΔK value (Evanno et al., 2005).

The approximate Bayesian computation (ABC) method, implemented in DIYABC v 2.1.0, was used to infer the demographic history of T. amurensis and T. japonica (Cornuet et al., 2014). On the basis of the findings from three approaches (STRUCTURE analysis, phylogenetic tree and PCoA) (Figures 2, 3) and the geographic distribution of the Tilia populations (Figure 1), four population groups were defined (Figure 4). Pop1 was composed of three Northeast China sites (DGL, MPS and DXG), Pop2 was composed of two North China sites (YDS and ZS), Pop3 was composed of three Shandong sites (TS, LU and LS), and Pop4 included the remaining three T. japonica sites (HS, MS and TMS). In scenario 1, Pop2 and Pop3 merged at time t1, Pop1 was present at time t2, and Pop4 was present at t3. The ancestors (NA) reached sufficient population size at t4. In scenario 2, Pop3 was assumed to have originated from an admixture of Pop2 and Pop4 at time t1. The rate of admixture of Pop2 and Pop3 was set as “ra”, and that of Pop3 and Pop4 was set as “1-ra”. Pop2 merged with Pop1 at time t2 and finally with Pop4 at time t3. NA reached sufficient population size at t4. In both scenarios 1 and 2, T. amurensis and T. japonica completed their divergence at time t3. Finally, in scenario 3, the four population groups merged simultaneously at t1. The prior values were set for the effective population size and divergence time estimated with a uniform distribution for all the parameters (Supplementary Table S2). The effective population size for each population group (N1, N2, N3, and N4) ranged from 10 to 20,000, with the constraint that all values were greater than or equal to the ancestral population size (NA). The divergence time (td) was set to a minimum of 10 generations and a maximum of 10,000 generations. Tilia plants can first flower when they are 12–40 years old, and the trees can live for >450 years (Pigott, 2012).With this in mind, we used a conservative generation time of 100 years. For each scenario, a mutation model was assumed with Mean mutation rate ranging from 10⁻⁴ to 10⁻³, Individual locus mutation rate ranging from 10⁻⁵ to 10⁻², Mean coefficient P ranging from 10⁻¹ to 3×10⁻¹, Individual locus coefficient P ranging from 10⁻² to 9, Mean SNI rate ranging from 10⁻⁸ to 10⁻⁴, and Individual locus SNI rate ranging from 10⁻⁹ to 10⁻³. Single-nucleotide insertion/deletion (SNI) rates were derived from 10,000 simulations generated for each putative scenario. The posterior probabilities were determined following 1,000,000 simulations for each scenario. The goodness-of-fit of the three scenarios was also evaluated by principal component analysis (PCA) using the option “model checking” in DIYABC.

The results of the SSR data analysis revealed the level of genetic diversity in the T. amurensis and T. japonica samples (Table 2). The level of genetic diversity in T. japonica was determined from Shannon’s information index (I), the observed heterozygosity (H_O), and the expected heterozygosity (H_E). The Shannon’s information index (I) of the T. amurensis populations varied between 0.889 and 1.231, whereas that of T. japonica populations ranged from 0.915 to 1.381. The observed heterozygosity (H_O) of T. amurensis populations ranged from 0.329 to 0.647, and the expected heterozygosity (H_E) ranged from 0.440 to 0.588. The observed heterozygosity (H_O) of T. japonica populations ranged from 0.426 to 0.606, and the expected heterozygosity (H_E) ranged from 0.520 to 0.657. The fixation indices (F_IS) were positive for all three T. japonica populations (Table 2).

The analysis of molecular variance (AMOVA) was implemented on the basis of SSR data and revealed the molecular variation attributable to differentiation among and within the populations for T. amurensis and T. japonica (Table 3). For T. amurensis, 12% of genetic variation exists among populations, while 88% exists within populations. In contrast, for T. japonica, 14% of the genetic variation was observed among populations, with the remainder (84%) existing within populations. All the results were highly significant (P < 0.001). The genetic differentiation coefficient (F_ST) among the populations of T. amurensis was 0.073, and the gene flow (N_m) was 3.191 (Table 3). The T. japonica populations, on the other hand, presented relatively high genetic differentiation (F_ST = 0.094) coefficients and relatively low gene flow (N_m = 2.403) (Table 3).

To better understand the relationships between T. amurensis and T. japonica, the gene flow (N_m) and genetic differentiation coefficient (F_ST) were estimated for all pairs of populations (Table 4). The highest gene flow (N_m = 14.568) and smallest genetic differentiation coefficient (F_ST = 0.017) were observed between ZS and YDS. The lowest gene flow (N_m = 0.670) and greatest genetic differentiation coefficient (F_ST = 0.272) were detected between DGL and TMS.

The Bayesian analysis implemented in STRUCTURE revealed a peak at K = 2. The K = 2 model was supported by the highest ΔK value (Figure 5), indicating that the most likely number of genetic clusters was 2 (ΔK = 250.88). In contrast, a lower ΔK value was observed at K = 3 (132.47). The bar plot for Tilia at K = 2 and 3 is shown in Figure 2. At K = 2, the 8 T. amurensis populations and the 3 T. japonica populations were clearly divided into two groups, with all T. amurensis populations categorized in the first group (red) and all T. japonica populations categorized in the second group (green). At K = 3, the populations show admixture between different genetic clusters, this being most evident within the MS population.

To validate the results of the population structure analysis, phylogenetic tree were constructed for eleven Tilia populations, and PCoA was performed (Figure 3). Both analyses showed that populations from TS, LU and LS clustered into one taxon, those from YDS and ZS clustered into one taxon, and those from DGL, MPS and DXG clustered into another taxon. The three T. japonica populations (MS, HS and TMS) are clustered together and clearly separated from the T. amurensis population. The results were in close agreement with the structural analysis (Figure 2).

The highest posterior probability after evaluation using the DIYABC analysis of Tilia was scenario 1, with p=0.8544 (Table 5; Figure 6). In Scenario 1, the ancestral population first diverged to give rise to the T. amurensis and T. japonica groups, with the T. amurensis group subsequently differentiating into three distinct population groups. The DIYABC model was evaluated via principal component analysis (PCA) in scenario 1 on the basis of 1000,000 simulations showing the fit between the observed and simulated datasets (Supplementary Figure S1). The posterior distribution of each effective population size was also analyzed for the four populations, Pop1, Pop2, Pop3 and Pop4, with median values of 3810, 5840, 9340 and 7610, respectively. The median value of NA was 75.8 for the ancestral population. The mean values of the divergence times t1, t2 and t3 were 287, 416, and 2090 generations (approximately 28.7, 41.6, 209 ka BP), respectively. The time scale of the ancestral population size change, t4, was 4880 generations (approximately 488 ka BP). The prior and posterior distribution patterns are presented in Supplementary Figure S2. The mean mutation rates of microsatellites and single-nucleotide insertions/deletions (SNIs) at the examined loci were 8.13×10–⁴ and 8.12×10^-6, respectively.

Based on the statistical analysis of hundreds of research papers on genetic variation and population genetic differentiation in different types of plants by Hamrick and Godt (1990), the average expected heterozygosity (H_O) for dicotyledonous plants is 0.136. For long-lived perennial woody plants, the average H_O is 0.177. For plants with an allogamous breeding system and animal-dependent pollination, the average H_O is 0.167. Combined with our previous findings, T. japonica (I = 1.181, H_O = 0.484, H_E = 0.602) and T. amurensis (I = 0.982, H_O = 0.516, H_E = 0.515) exhibit comparable levels of genetic diversity, both showing relatively high genetic variation at the population level. It is noteworthy that the conclusions drawn by Hamrick and Godt were based on allozyme data. However, when compared with genetic diversity revealed by molecular markers in species such as the ecologically similar Picea obovata (H_O = 0.408; H_E = 0.423) (Kravchenko et al., 2016) and congeneric species including Tilia sibirica (H_E = 0.318) (Logan et al., 2018)and Tilia tomentosa (I = 0.50) (Gabur et al., 2019), T. japonica and T. amurensis still demonstrate the relatively high level of genetic diversity. According to the statistical analysis by Hamrick and Godt (1990) of hundreds of research papers on genetic variation and population genetic differentiation in different types of plants, the factors influencing the magnitude of genetic variation at the species level are, in order: taxonomic status, distribution range, life form, breeding system, and seed dispersal mechanism. At the population level, the influencing factors are: breeding system, distribution range, life history, taxonomic status, and seed dispersal mechanism. Woody tree species with a wide distribution range, long lifespan, predominantly allogamous breeding system, entomophilous pollination, frequent gene exchange, and abundant seeds that occur in late successional communities exhibit high levels of genetic variation (Hamrick and Godt, 1990). Woody tree species with a wide distribution range, long lifespan, predominantly allogamous breeding system, entomophilous pollination, and frequent gene exchange exhibit high levels of genetic variation (Hamrick and Godt, 1990).

The high genetic diversity of T. japonica may be determined by its life history and reproductive system characteristics. Tóth et al. (2022) also suggested that species of the genus Tilia are able to maintain a high level of genetic diversity due to their abundance during warmer periods, outcrossing mating system, and increased lifespan. Entomophilous pollination of Tilia and the intertwining of intraspecific differentiation and germplasm infiltration phenomena provide opportunities for gene recombination, which leads to high levels of gene flow in T. japonica (Ren and Ya, 1995). The reproductive system of T. japonica, which is predominantly characterized by sexual reproduction via entomophily, ensures a high level of genetic diversity within its populations (Boo and Park, 2016; Barker et al., 2022). Compared with the narrow distribution ranges of its congeners T. sibirica and T. argentea, T. japonica has a relatively wide distribution in eastern China and Japan. This broader distribution range may have led to higher levels of genetic diversity in T. japonica than in other plants of the same genus. Even though the Ne value was lower than the Na value in all T. japonica populations, genetic drift did not lead to a reduction in genetic diversity levels. Our results revealed that the genetic diversity of the TMS population was significantly lower than that of the other T. japonica populations. The inbreeding coefficients of T. japonica populations were greater than 0, and that of the TMS (F_IS = 0.221) population was the largest (Table 2), suggesting the possibility of selfing in T. japonica. Selfing reduces genetic diversity by decreasing genetic variation and weakening gene flow among populations, which may account for the relatively low genetic diversity observed in TMS populations (Jullien et al., 2019; Clo et al., 2019).

According to the clear rules of genetic differentiation between populations (F_ST < 0.05, low; 0.05< F_ST < 0.15, medium; F_ST >0.15, high) defined by Wright (1978) and the results of the current study, there is a moderate degree of genetic differentiation between populations of both T. amurensis (F_ST = 0.073) and T. japonica (F_ST = 0.094). Moreover, AMOVA revealed that most of the genetic variation was distributed within the populations of T. amurensis and T. japonica. Genetic differentiation among populations is strongly influenced by gene flow and genetic drift (Andrews, 2010). When N_m > 1, there is sufficient gene flow between populations to exert a homogenizing effect, counteracting genetic drift and preventing differentiation. When N_m > 4, there is high gene flow between populations, and the populations can be regarded as a single unit of random mating (Wright, 1931). Both populations of T. amurensis (N_m = 3.191) and T. japonica (N_m = 2.403) exhibit strong gene flow, which can reduce local variation and prevent adaptive differentiation, thereby leading to lower genetic differentiation among their populations (Xiang et al., 2015; Wang et al., 2025). According to the study by Hamrick and Loveless (2019), the level of gene flow in plants is significantly influenced by their reproductive systems and seed dispersal mechanisms. Specifically, the mode of plant reproduction (such as selfing or outcrossing) and the mode of seed dispersal (such as wind dispersal or animal dispersal) are key factors in determining the intensity of gene flow. For instance, plants that predominantly rely on outcrossing typically exhibit higher levels of gene flow, whereas selfing plants tend to have relatively lower gene flow (Ballesteros-Mejia et al., 2016; Noreen et al., 2016; Gamba and Muchhala, 2023). Tilia species primarily reproduce through outcrossing via insect-mediated pollination, which can be considered the main factor contributing to the gene flow and population structure of populations of T. amurensis and T. japonica.

The results of the STRUCTURE analysis of 242 individuals revealed that the T. amurensis populations and the T. japonica populations could be well distinguished when the K value was 2, consistent with the two clustering results revealed by the phylogenetic tree and PCoA analysis. All individuals of T. amurensis have genetic components from T. japonica, and all individuals of the MS population in T. japonica, as well as several individuals of HS and TMS, also experienced genetic infiltration from T. amurensis. Relative to the HS and TMS populations, due to geographical proximity to the distribution area of T. amurensis, the MS population and T. amurensis populations exhibit greater gene flow and a lower degree of differentiation. This may be a factor leading to the MS population experiencing greater more infiltration of genes from T. amurensis. When the value of K is 3, all individuals of T. amurensis have genetic components from T. japonica, and all individuals of T. japonica have genetic components from T. amurensis. Among them, DGL, MPS, and DXG populations had the lowest degree of gradual infiltration, while easily distinguishable from other T. amurensis populations. The main reason is also that DGL, MPS, and DXG populations are geographically furthest away from T. japonica distribution area. The presence of significant gene infiltration between the two species implies that the two species may have originated from the same ancestor (Sork, 2016; Ribicoff et al., 2025).

The ABC analysis in DIYABC revealed that T. amurensis and T. japonica share the same ancestral origin in this scenario, after which T. amurensis diverged into three population groups. Assuming a generation time of 100 years for Tilia and considering the summary statistics, the parameter set and the SSR markers used, we can infer a Middle Pleistocene split between T. amurensis and T. japonica, approximately 2090 generations ago (approximately 209 ka BP). In the Late Pleistocene, the T. amurensis population group in Northeast China (DGL, MPS and DXG) diverged from those in North China (ZS and YDS), approximately 416 generations ago (approximately 41.6 ka BP). Additionally, divergence from the Shandong population group (TS, LU and LS) occurred approximately 287 generations ago (approximately 28.7 ka BP). Some of the earliest evidence for putative Tilia species dates back to the Tertiary (Wolfe and Wehr, 1987; Pigott, 2012), and the approximate age of Tilia was inferred to be 17 million years (Richardson et al., 2015). A recently approximated age for the divergence of Tilia was inferred to be approximately 447 ka BP, also during the Middle Pleistocene (Logan et al., 2018). Some studies suggest that climatic fluctuations from the Middle to Late Pleistocene strongly affected tree species, causing lineage divergence. For example, Quercus spp diverged during the Early Pleistocene, Populus alba and P. davidiana diverged in the Middle Pleistocene as they spread southward from the Qinling Mountains, and the North American tree species Populus balsamifera and P. trichocarpa diverged during the Late Pleistocene at approximately ∼75 ka BP (Levsen et al., 2012; Richardson et al., 2015; Hou and Li, 2022). During the Middle and Late Pleistocene, the pronounced alternation of glacial and interglacial periods and the expansion and contraction of the ice sheet greatly influenced the population dynamics of Tilia (Groisman et al., 2013). During the Ice Age, when the climate is cold and dry, the distribution of linden populations shrinks to lower altitudes and latitudes, while during the interglacial period, when the climate warms, Tilia populations expand rapidly (Bolikhovskaya and Shunkov, 2014). T. amurensis taxa that had been migrating southward since the Miocene were strongly affected in their genetic structure during this repeated expansion and contraction, leading to genetic isolation and differences, which may have subsequently led to a Middle/Late Pleistocene split. The formation of today’s Tilia range is the result of redispersal after the end of the most recent ice age (approximately 15 ka BP).

Both T. amurensis and T. japonica are vulnerable to habitat loss and overexploitation. Habitat fragmentation caused by anthropogenic disturbances threatens the survival of Tilia (Mu and Liu, 2007; Zhang et al., 2024). Currently, there are no effective conservation measures in place for T. amurensis and T. japonica. Our results indicate that the level of genetic diversity is high in T. japonica and that genetic variation occurs mainly within populations of T. amurensis and T. japonica. Therefore, in situ conservation of these populations is the most effective strategy to preserve their species diversity. The HS and MS populations of T. japonica maintain high genetic diversity and should be prioritized for in situ conservation. The complete biological conservation strategy should combine in situ and translocated conservation (Volis and Blecher, 2010; Dulloo et al., 2010; Duman et al., 2024). Therefore, to avoid inbreeding decline and population degradation, relocation protection and seed collection should also be strengthened. The shared ancestral background and the presence of bidirectional gene introgression between T. amurensis and T. japonica suggest that genetic connectivity between these two species should be considered in their conservation and management. To prevent population degradation and loss of genetic resources due to excessive hybridization between the two species, the protection of the three T. amurensis populations in Northeast China (DGL, MPS and DXG) and the TMS population of T. japonica, which have a purer genetic structure with less genetic drift, should be strengthened. Furthermore, conserving the genetic diversity of these two species requires attention not only to their respective core distribution areas but also to their zones of distributional overlap, which represent critical regions for gene flow.