Leaf samples of F. hayatae were collected from three native habitats in Taiwan and six in mainland China ( Table 1 ). In Taiwan (TW region), sampling was conducted in three populations separated by the Lanyang River: 101 individuals from the Tongshan (TS) population (south of the river), 16 from the Ayushan Mountain West Peak (AW) population (north of the river), and 40 from the Peichatienshan (PCTS) population (north of the river). In mainland China, samples were collected from two regions: the Zhejiang Province (ZJ region), with 8 individuals from Sihaishan (CH), 23 from Qingliangfeng (DB), and 3 from Tianmu Mountain (TM); and the Hubei and Sichuan Provinces (HS region), with 16 individuals from Shennongjia (SNJ), 18 from Big-small Langou Nature Reserve (LG), and 37 from Micangshan (MC). Two samples each of F. multinervis (UL) and F. sylvatica (UK) served as outgroups. A total of 249 samples were used in the analysis ( Table 1 ). Fresh young leaves were collected, rapidly dehydrated using silica gel, and stored for subsequent genomic DNA extraction.

Genomic DNA was extracted using a modified CTAB method (Doyle and Doyle, 1987), dissolved in 200 μl of Tris-EDTA (TE) buffer, and stored at -20°C. DNA quality was assessed via 1% agarose gel electrophoresis using a Lambda marker (200 μg/mL, Promega) and diluted to 10 ng/μl in TE buffer for subsequent experiments.

To assess the levels of genetic variability and population structure in F. hayatae, genotyping was performed using ten previously characterized polymorphic microsatellite markers described by Ju et al. (2012) (Ju et al., 2012). Although 249 individuals were sampled for this study, DNA of sufficient quality and quantity for the sensitive microsatellite analysis could only be obtained from 184 samples, as detailed in Table 1 . These 184 samples ( Table 1 ) underwent PCR amplification using a Labnet MultiGene 96-well Gradient Thermal Cycler with the following parameters: initial denaturation at 94°C for 2 min, followed by 35–40 cycles of 94°C for 45s, 54-64°C for 1 min, and 72°C for 30s, concluding with a final extension at 72°C for 7 min. Microsatellite genotyping followed the methodology outlined by Mantiquilla et al. (2022) (Mantiquilla et al., 2022).

To evaluate the phylogeographic structure and dynamics of F. hayatae populations in Taiwan and mainland China, the atpB-rbcL regions of the plastid genome were amplified using the forward primer 5’- ACCGGACCAATGATTTGAGCC-3’ and the reverse primer 5’-TACAGTTGTCCATGTACCAG-3’. A PCR reaction mixture totaling 50 µL comprised 1 µL of template DNA, 5 µL of 10× reaction buffer, 5 µL of dNTP mix (2 mM), 5 µL of each forward and reverse primer (2 mM), 1 µL of Taq polymerase (0.2 U/µL; Promega), and 28 µL of sterile water. The PCR amplification protocol included an initial denaturation at 95°C for 2 min, followed by 38 cycles of denaturation at 94°C for 45 s, annealing at 54°C for 1 min, and extension at 72°C for 1 min, culminating in a final extension at 72°C for 7 min. All PCR products were purified using the HiYield™ Gel/PCR DNA Fragments Extraction Kit (RBC Bioscience). Both strands of the purified PCR fragments were sequenced using the ABI PRISM 3730XL DNA sequencer (Applied Biosystems, Foster City, CA, USA).

The resulting forward and reverse sequence data were assembled into contigs using SeqMan software from the DNASTAR package (DNASTAR Inc., Madison, WI, USA). Multiple sequences were aligned and retrieved from the GenBank database utilizing MEGA 11 (Kumar et al., 2018) for subsequent phylogenetic analyses. The optimal nucleotide substitution model was determined using jModelTest v2.1.7 (Darriba et al., 2012) based on the Bayesian Information Criterion (BIC), identifying the HKY + G model, which was subsequently applied to each Beast analysis conducted with BEAST v1.10.4 (Suchard et al., 2018). All sequences of atpB-rbcL regions submitted to GenBank (GenBank accession numbers: LM992847-LM992854).

Microsatellite data revealed varying levels of genetic diversity among F. hayatae populations across three regions (Taiwan [TW], Zhejiang [ZJ], and Hubei/Sichuan [HS]) ( Table 2 ). Hereafter, we refer to TW populations as “island populations” and Chinese populations (ZJ and HS) as “mainland populations”. The island populations showed greater genetic diversity, with higher numbers of different alleles (A = 12.37), number of effective alleles (A_E = 6.55), allelic richness (Ar = 4.06), and private allelic richness (A_P = 2.14) comparing to the mainland populations. Among the island populations, PCTS had the highest genetic diversity (A = 18.80, A_E = 8.85, Ar = 4.57, A_P = 2.23), whereas the AW population had lower genetic diversity (A = 5.40, A_E = 3.70, Ar = 3.20, A_P = 1.77). Among mainland populations, DB showed relatively higher genetic diversity (A = 5.90, A_E = 3.27, Ar = 2.96, Ap = 0.66).

Hardy-Weinberg equilibrium (HWE) tests revealed significant deviations from equilibrium in all island populations, with observed heterozygosities (H_O ) consistently lower than expected heterozygosities (H_E ). Similarly, all mainland populations exhibited lower H_O than H_E and deviated significantly from HWE, except the TM population, which showed no significant deviation from HWE ( Table 2 ). Analysis of fixation index (F_IS ) values revealed a significantly higher F_IS value in island populations (0.96), with the AW population exhibiting the highest value (F_IS = 0.98), suggesting a high degree of inbreeding in the island populations. While most mainland populations showed significant F_IS values, the CH, TM, and SNJ populations exhibited low and non-significant values ( Table 2 ). Taken together, F. hayatae island populations showed higher genetic diversity than mainland populations, along with significant levels of inbreeding.

Examination of chloroplast DNA data across nine F. hayatae populations in island and mainland revealed six distinct haplotypes ( Table 1 , Figure 1 ). The analysis demonstrated that island populations harbored exclusive haplotypes (Hap3 and Hap7), distinguishing them from mainland populations. Notably, island-specific haplotype Hap3 was found solely in the TS population, absent from other island populations (PCTS and AW), underscoring the unique genetic composition of the TS population. Four haplotypes were identified in mainland populations, with Hap6 demonstrating widespread distribution across ZJ and HS regions, while CH, DB, and SNJ populations possessed unique haplotypes (Hap4, Hap5 and Hap8) ( Table 1 , Figure 1 ).

Further genetic analysis revealed that the island populations displayed the highest genetic diversity (H_T , Hr, Hp, Hd = 0.84, 1.55, 1.55, 0.46) and nucleotide diversity (π = 2.48 × 10^-3, θw = 0.95 × 10^-3), while the HS region showed the lowest (π = 0.26 × 10^-3, θw = 0.22 × 10^-3) based on chloroplast DNA data ( Supplementary Table S1 ). The observed reduced overall nucleotide diversity of F. hayatae populations indicates potential historical genetic bottlenecks or recent population isolation events. Furthermore, the high G_ST and N_ST values (>0.5) detected in all populations ( Supplementary Table S1 ) indicate significant genetic differentiation and fragmentation among these populations.

Hierarchical AMOVA was performed using both chloroplast DNA and microsatellite data to assess genetic differentiation across various hierarchical levels defined by geographical regions and genetic groups ( Supplementary Table S2 ).

Using chloroplast DNA sequences, analyses based on different spatial grouping models—K _geo = 2 (Taiwan & China), K _geo = 3 (TW, HS, and ZJ regions), and K _geo = 9 (all nine populations)—revealed that the majority of genetic variation was accumulated among populations within groups (75.97%, 81.95%, and 99.13% for the K _geo = 2, K _geo = 3, and K _geo = 9 patterns, respectively), with both among and within population variability being high and statistically significant (Φ _SC and Φ _ST > 0.99). In contrast, analysis based on genetic groupings (K _S = 2, as determined by STRUCTURE analysis; TS and PCTS in one cluster, and AW, CH, DB, TM, SNJ, MC, and LG in the other) indicated that most of the genetic variation (71.31%) accumulated among groups, with consistently high and significant variability at all hierarchical levels (Φ _CT = 0.71, Φ _SC = 0.98, and Φ _ST = 0.99) ( Supplementary Table S2 ).

Interestingly, microsatellite polymorphisms revealed a distinctive pattern of genetic variation ( Supplementary Table S2 ). Regardless of whether populations were grouped by geographical models (K _geo = 2, K _geo = 3, or K _geo = 9) or by genetic clusters (K _S = 3, determined by STRUCTURE analysis: TS and PCTS in one cluster, AW in another, and CH, DB, TM, SNJ, MC, LG in the third; K _SA = 4, determined by SAMOVA analysis: PCTS, AW, and TS each in separate clusters, with CH, DB, TM, SNJ, MC, LG in the fourth), the analyses consistently showed that the majority of genetic variation was accumulated within populations (69.04%, 69.78%, 73.62%, 72.17%, and 72.68% for K _geo = 2, K _geo = 3, K _geo = 9, K _S = 3, K _SA = 4 group patterns). Notably, significant variability was observed across groupings, among populations, and within populations (Φ _CT , Φ _SC , and Φ _ST ranging from 0.06 to 0.30). These findings indicate that while chloroplast DNA primarily captures differences between populations, microsatellite markers are better at detecting finer-scale intra-population genetic variability ( Supplementary Table S2 ).

After confirming significant genetic variation among populations, we investigated the population structure of F. hayatae using two analyses: SAMOVA and STRUCTURE. SAMOVA was used as a preliminary analysis with both microsatellite and chloroplast DNA data. Based on F_CT as the criterion for optimal grouping, the chloroplast DNA data did not show a clear clustering pattern ( Supplementary Table S3 ). In contrast, the microsatellite data revealed four distinct clusters (K=4), with the highest F _CT value of 0.55 (p < 0.05) ( Supplementary Table S3 ; Figure 2A ). Under this model, the three populations from the island (PCTS, AW, and TS) were clearly distinguished from each other, while all populations from the mainland (CH, DB, TM, SNJ, MC, LG) grouped together. This initial clustering analysis supports the classification of island and mainland populations into distinct genetic groups.

We also used STRUCTURE to further examine the genetic structure of the populations. The optimal number of clusters (K) was determined by the value of ΔK. For chloroplast DNA, the ΔK method identified K = 2 (ΔK = 3818.27), with the TS population from island forming a distinct cluster (Cluster C2) separate from the other groups ( Supplementary Table S3 ; Figures 2A, B ). In contrast, for the SSR data, the optimal number of clusters was K = 3 (ΔK = 1882.04). In the island populations (PCTS, AW, and TS), two distinct genetic clusters were observed (S2, and S3). AW was composed solely of cluster S3, while the PCTS and TS populations displayed more varied genetic components (cluster S2 and S3). All mainland populations were homogeneous, forming a single cluster (Cluster S1) ( Supplementary Table S3 ; Figures 2A, C ). These clustering analyses revealed a distinct population genetic structure with a clear separation between the island and mainland, consistent across both genetic markers.

The haplotype network of F. hayatae showed a distinct grouping of island populations ( Figure 1C ). These island populations possessed unique haplotypes, Hap3 and Hap7. Hap3, found only in the TS population and separated from others by more than three mutational steps, suggested long-term isolation or independent evolution. In contrast, Hap7 was closer to mainland haplotypes, suggesting a historical gene flow between island and mainland populations. In addition, high genetic differentiation was evident, with limited haplotype sharing among populations except for Hap6, which was common in the ZJ and HS regions ( Figure 1C ).

Bayesian inference (BI) analysis of F. hayatae chloroplast DNA haplotypes revealed three monophyletic clades ( Figure 3 ). Clade 1, which comprised the TS population (Hap3), was moderately supported (PP = 0.75), while the other two clades showed stronger support. In Bayesian phylogenetics, posterior probabilities of 0.95 or higher are conventionally considered strong statistical support, whereas values between 0.70 and 0.94 represent moderate support. Clade 3 consisted solely of the DB population (Hap5), emerged as the earliest diverging lineage among mainland populations. Clade 2 included both island and mainland populations but showed no clear geographic structure. Notably, the island populations PCTS and AW clustered most closely with the mainland CH population, which is geographically closest to Taiwan ( Figure 3 ).

Divergence time estimates, based on a substitution rate of 1.52 × 10^-9 s/s/y, provided further insights into the evolutionary timeline of F. hayatae populations ( Figure 3 ). Divergence from other Fagus species occurred approximately 49.37 million years ago (Ma). The TS population was the first to diverge from other F. hayatae populations around 44.61 Ma, followed by the DB population, which separated from other mainland Chinese F. hayatae populations around 17.27 Ma. The PCTS and AW populations diverged from the CH population approximately 8.36 Ma ( Figure 3 ). It is important to note, however, that these divergence times are based on a single chloroplast DNA intergenic region (atpB-rbcL). Such single-locus estimates can provide only an approximate timeline and are subject to considerable uncertainty, as reflected by the wide 95% HPD intervals in our analysis ( Figure 3 ). These estimations should therefore be interpreted with caution.

Tajima’s D and Fu and Li’s D* tests were conducted to evaluate the neutrality of genomic regions by comparing observed nucleotide diversity. The results indicated that values ranged from 0.23 to 3.18 for Tajima’s D and from 0.52 to 1.00 for Fu and Li’s D*. Notably, significantly positive Tajima’s D values were observed exclusively for the island populations, suggesting the influence of population bottlenecks or balancing selection in these groups. ( Supplementary Table S1 ).

To assess demographic changes, a mismatch distribution analysis was performed by calculating both the sum of squared deviations (SSD) and Harpending’s raggedness index (H _Rag). Most populations showed low and statistically significant values under both demographic and spatial expansion models, indicating that these F. hayatae populations did not experience significant expansion. Instead, factors such as bottlenecks, population structure, or other historical events likely influenced their current genetic patterns. An exception was found in the SNJ population of the HS region, where the mismatch distribution supported an expansion model ( Supplementary Table S4 ).

Further insights into historical demographic changes were provided by a Bayesian skyline plot, which revealed a recent decline in effective population size in the island populations, suggesting that local genetic diversity might be at risk. In contrast, the mainland populations remained relatively constant and stable in their population sizes. Overall, when considering F. hayatae as a whole, the effective population size appears to have decreased over the past million years ( Supplementary Figure S1 ).

After confirming significant genetic differentiation between island and mainland populations, we applied a coalescent-based Isolation-with-Migration (IM) model to estimate gene flow rates and effective population sizes between the island (population 1) and mainland (population 2) groups ( Supplementary Table S5 ; Figure 4 ). The IMa analysis showed that the ancestral effective population sizes (N_A) were considerably larger than the current sizes with estimates of 2.37 × 10⁶ (95% CI: 2.47 × 10⁶–2.28 × 10⁶) based on chloroplast DNA, and 2.81 × 10⁶ (95% CI: 1.23 × 10⁶–5.03 × 10⁷) based on microsatellite data. By comparison, current effective population sizes were much smaller: 2.08 × 10⁵ (95% CI: 2.17 × 10⁵–2.00 × 10⁵) for the island and 3.98 × 10⁵ (95% CI: 4.14 × 10⁵–3.83 × 10⁵) for the mainland based on chloroplast DNA, and 2,490 (95% CI: 1,090–44,500) for the island and 2,040 (95% CI: 891–36,500) for the mainland based on microsatellites ( Supplementary Table S5 ; Figure 4 ). Historical gene flow estimates were generally low. Microsatellite data indicated slightly higher gene flow from Taiwan (TW) to mainland China (ZJ and HS regions) (M₁→₂ = 0.16, 95% CI: 0.01–0.35) than in the reverse direction (M₂→₁ = 0.02, 95% CI: 0.00–0.05), whereas chloroplast DNA data showed no evidence of gene flow in either direction ( Supplementary Table S5 , Figure 4 ).

Ancestral range reconstruction using RASP with chloroplast DNA data provided insights into the historical biogeography (vicariance and dispersal events) of F. hayatae across mainland and island regions ( Figure 3 ). The results revealed several dispersal events among ancestral regions on the mainland at nodes 9, 11, and 12, with ancestors originating from the ZJ region. Additionally, two vicariance events were identified between the island and mainland, with no evidence of gene flow between them ( Figure 3 ).

These findings suggest that geographic isolation has played a key role in shaping the evolutionary history and genetic divergence of F. hayatae between mainland and island populations.

Genetic diversity analysis of F. hayatae revealed generally low diversity across the species, with island populations displaying higher genetic diversity but also much higher inbreeding coefficients than mainland populations. These distinct evolutionary patterns between island and mainland groups reflect their unique evolutionary histories and ecological conditions. These results support the hypothesis that Taiwan serves as a genetic refugium for F. hayatae. According to refugia theory, geographical isolation and climatic stability enable specific regions to become species refuges, maintaining higher genetic diversity (Hewitt, 2000). Climatic refugia are linked to areas of local or regional climate stability and biodiversity hotspots, where stable climates allow both ancient species and newly evolving lineages to survive and thrive (Harrison and Noss, 2017; Keppel et al., 2018). García et al. (2022) confirmed that stable climatic conditions and topographical complexity contribute to the accumulation of evolutionary uniqueness.

The elevated genetic diversity in Taiwanese populations aligns with refugia theory, attributed to the island’s geographic isolation, stable climate, and complex topography, which have facilitated the maintenance and accumulation of genetic variation. Previous community surveys indicated that Taiwanese F. hayatae populations prefer warm, humid environments with minimal seasonal variation (Shen et al., 2015), further supporting Taiwan’s role as a climatic refuge that enables long-term independent evolution and preservation of higher genetic diversity.

Our phylogenetic analysis further supports Taiwan’s role as a genetic refugium. The TS population in Taiwan represents the earliest diverging evolutionary lineage and possesses unique haplotypes, highlighting its genetic distinctiveness. These findings support Taiwan’s role as a climatic refugium and demonstrate how stable climate and topographical conditions have fostered evolutionary uniqueness (García et al., 2022). However, our analysis revealed that Taiwanese populations have higher inbreeding coefficients, primarily due to limited gene flow between isolated F. hayatae populations caused by long-term geographical separation. Habitat fragmentation and smaller population sizes have increased the risk of inbreeding. These factors affect genetic diversity maintenance and may reduce the populations’ long-term adaptive potential. Population dynamics analysis also revealed that current effective population sizes of Taiwanese F. hayatae populations have significantly decreased compared to their ancestral levels, indicating a recent genetic decline.

In contrast, mainland populations of F. hayatae showed lower genetic diversity compared to Taiwanese populations, possibly due to historical bottleneck effects and greater climate fluctuations in their habitat. Previous community surveys indicated that mainland Chinese F. hayatae populations are adapted to colder and drier climates (Shen et al., 2015), which may reflect the region’s greater climatic variability and its weaker effect on maintaining genetic diversity. However, despite their lower genetic diversity, mainland populations had lower inbreeding coefficients than Taiwanese populations. This may be because they occupy wider geographic areas with more varied environments, which likely promotes gene flow between populations and helps reduce inbreeding. Notably, our results indicated that some mainland populations, such as the TM population, showed no significant deviation from Hardy-Weinberg equilibrium, indicating stable population dynamics.

Demographic analyses confirmed that while most F. hayatae populations have experienced bottleneck effects or balancing selection, mainland populations exhibit stable population dynamics. Demographic analysis revealed that the SNJ population showed signs of expansion, and Bayesian skyline plots confirmed stable population sizes across mainland regions, demonstrating sustained population stability of mainland populations. However, IMa analysis revealed a significant decrease in current effective population sizes compared to ancestral populations in mainland populations. This indicates that although mainland populations maintain gene flow and show lower inbreeding levels, they have experienced reduced genetic diversity due to historical bottleneck effects and climate fluctuations.

Our study uncovered clear genetic differentiation between Taiwanese and mainland Chinese populations of F. hayatae, as evidenced by four distinct genetic clusters identified through microsatellite data and a unique genetic lineage of the Taiwanese TS population revealed by chloroplast DNA analysis. Bayesian phylogenetic analysis further delineated three well-supported monophyletic clades, with early divergence observed in the TS population (Clade 1) and the DB mainland population (Clade 3), along with a mixed cluster (Clade 2) containing both island and mainland individuals. The presence of unique haplotypes in Taiwan, particularly within the TS population, suggested that island populations had undergone independent evolutionary trajectories, likely driven by local adaptation and historical vicariance. Additionally, IMa gene flow analysis revealed extremely low bidirectional gene flow between Taiwan and mainland China, while chloroplast DNA data showed no evidence of such exchange, highlighting the strong effect of geographic isolation on the genetic structure.

The fine-scale spatial genetic structure of F. hayatae, characterized by increased inbreeding and distinct genetic patterns in small and isolated Taiwanese populations, is primarily shaped by gene flow patterns through pollen and seed dispersal mechanisms. Restricted pollen dispersal often leads to consanguineous mating among nearby individuals. Although wind-dispersed pollen exhibits greater dispersal distances, its distribution patterns remain stochastic and difficult to predict, while seed dispersal typically has a more pronounced impact on genetic structuring due to its localized nature (Berg and Hamrick, 1994; Angbonda et al., 2021). When both mechanisms are limited, as observed in F. hayatae, distinct spatial genetic patterns emerge alongside increased inbreeding and reduced diversity. In Taiwan, the dispersal of seeds and pollen in F. hayatae is limited. While wind serves as the primary pollination agent, the persistently high humidity and heavy rainfall in its habitat inhibit effective pollen dispersal (Chen et al., 2011), further reinforcing genetic isolation and promoting inbreeding. These environmental constraints, together with random genetic drift, accelerate divergence and reduce genetic diversity in these small, isolated populations (Ellstrand and Elam, 1993).

Integrating these findings, our results highlight that ecological constraints, historical isolation, restricted seed and pollen dispersal, and ongoing genetic drift have created strong barriers to gene flow, resulting in pronounced genetic structuring across the species’ range. The concordance between genetic patterns and ecological observations underscores the necessity to maintain or restore connectivity among fragmented populations as a core conservation priority for F. hayatae.

The distinct biogeographical characteristics of China/East Asia, with its complex topography and unique east-west oriented mountain ranges, have served as crucial buffers against climatic fluctuations during the late Neogene and Quaternary periods, particularly in south-central China and adjacent regions where Tang et al. (2018) documented exceptional climatic stability compared to surrounding Asian territories (Manchester and Tiffney, 2001; López-Pujol et al., 2011; Tang et al., 2018). This long-term environmental stability, characterized by both geological and climatic constancy, has been instrumental in promoting species diversification and facilitating evolutionary processes among temperate forest floras. Similarly, Jiang et al. (2020) demonstrated that genetic differentiation patterns in Fagus species closely correlate with these geological and climatic stability patterns (Jiang et al., 2020). Our study on F. hayatae provides further evidence for this correlation.

Our phylogeographic analysis, based on the cpDNA atpB-rbcL region, reconstructs the historical dynamics of F. hayatae and provides a preliminary timeline of divergence events. The results suggest that F. hayatae diverged from other Fagus species ca. 49.37 Ma, during the middle Eocene. A significant vicariance event is estimated to have occurred ca. 44.61 Ma, separating the Taiwanese TS population lineage from other populations. We initially attributed this deep split to major climatic shifts during the Late Eocene. However, we must explicitly acknowledge the limitations inherent in using a single, non-coding chloroplast DNA marker for divergence time estimation. Such single-locus approaches can be influenced by factors like incomplete lineage sorting, potential calibration errors, or insufficient fossil constraints, which can lead to considerable uncertainty in age estimates. This uncertainty is reflected in the wide 95% highest posterior density (HPD) intervals in our analysis.

Crucially, the estimated Eocene divergence of the TS lineage presents a geological contradiction, as it significantly predates the tectonic uplift of Taiwan as a substantial landmass, which began around the mid-Miocene (~9 Ma) and intensified in the last 5 million years (Sibuet and Hsu, 2004). It is therefore unlikely that this lineage was isolated in situ on Taiwan since the Eocene. One plausible, though speculative, scenario is that this ancient cpDNA lineage persisted on the Asian mainland for tens of millions of years before one of its descendant populations colonized the newly formed island of Taiwan.

This discrepancy is further highlighted when comparing our results with Jiang et al. (2020) studies using 28 nuclear single/low-copy loci, estimated a considerably more recent divergence time for F. hayatae from its sister species (F. crenata) at approximately 7.8 Ma (95% HPD: 6.49–9.41 Ma). Furthermore, their analysis placed the major radiation of the East Asian, North American, and West Eurasian lineages of subgenus Fagus around the late Miocene (10.9–10.4 Ma). This timeframe aligns much better with the geological history of Taiwan’s formation and major climatic cooling events of the Neogene that drove species differentiation (Wen, 1999; Milne and Abbott, 2002). Therefore, while our cpDNA data compellingly identifies the TS population as a uniquely ancient and distinct maternal lineage, the ~44 Ma divergence estimate should be interpreted with caution and viewed as an approximate timeline requiring further validation. The multi-locus nuclear data from Jiang et al. (2022) likely provide a more realistic estimate for the species-level divergence events. Future research integrating genomic data from multiple nuclear loci will be essential to resolve this conflict between organellar and nuclear genomes and to establish a more precise and robust evolutionary timeline for F. hayatae.

In contrast, the estimated divergence of the AW and PCTS populations at 8.36 Ma aligns well with the geological history of Taiwan and major climatic shifts in East Asia. The island populations, AW and PCTS, became completely isolated from mainland populations during the late Miocene to early Pliocene. This period was marked by intensified global cooling, the gradual disconnection of the Bering Land Bridge, and regional species isolation due to land uplift and climate diversification in East Asia (Wen, 1999; Milne and Abbott, 2002; Ian Milne, 2006). Regional geological events and climatic shifts drove species isolation and differentiation, which ultimately shaped contemporary population structures. Notably, this period coincides with the formation of the island of Taiwan. Around 9 million years ago, the collision of the Philippine Sea Plate with the Eurasian Plate initiated the island’s gradual uplift (Sibuet and Hsu, 2004). This geological process not only elevated the terrain but also significantly altered local environmental conditions. The accelerated tectonic uplift of Taiwan further restricted gene flow among populations and led to regional species isolation. Consequently, the gradual formation of Taiwan island acted as a physical barrier that obstructed migration routes and, in conjunction with climatic changes, established both geographic and ecological obstacles between the AW and PCTS populations in Taiwan and those in China, thereby promoting species differentiation and evolution. These geological and ecological environmental changes not only explain F. hayatae’s distribution patterns but also reveal the profound impact of environmental changes on its population’s genetic structure and differentiation.

Our analysis revealed significant genetic divergence between island and mainland populations of F. hayatae, particularly the distinct genetic signature of the TS lineage, which has considerable implications for taxonomy and conservation. Revisiting the classification controversy mentioned earlier with our new genetic evidence and recent literature provides valuable insights. Li et al. (2023) integrated molecular and morphological data to synonymize the previously ambiguous F. chienii with F. hayatae, resolving a longstanding taxonomic question regarding mainland Chinese populations. This taxonomic consolidation aligns with our treatment of all mainland samples as members of a single, genetically structured F. hayatae species.

Nevertheless, our genetic findings clearly demonstrate that this taxonomic unification should not diminish the significant evolutionary distinctiveness of Taiwanese populations. The pronounced divergence of the TS population’s chloroplast haplotype, the clear genetic differentiation detected by both SAMOVA and STRUCTURE analyses, and the virtually nonexistent gene flow between Taiwan and mainland China collectively indicate that Taiwanese populations have followed a separate evolutionary pathway.

These findings raise the question of whether Taiwanese populations warrant classification as a separate subspecies or species. While formal taxonomic revision exceeds the scope of our study, our data strongly support designating Taiwanese populations, especially the TS population, as an Evolutionarily Significant Unit (ESU). This designation applies to population segments that are reproductively isolated and represent critical components of a species’ evolutionary heritage. The unique genetic composition, significant divergence, and prolonged isolation of Taiwanese lineages precisely fulfill these criteria. This genetic distinctiveness, combined with elevated inbreeding coefficients in Taiwanese populations, indicates they are following an independent evolutionary trajectory.

From a conservation perspective, this distinction is crucial. Managing F. hayatae as a single conservation unit across its entire distribution would be inappropriate and could result in losing valuable genetic diversity. We recommend implementing separate conservation strategies for mainland and Taiwanese populations. Additionally, the TS population, representing the most ancient and divergent lineage identified in our study, should receive priority for in-situ conservation efforts to safeguard its unique and irreplaceable genetic diversity.