Plastome data provides new insights into population differentiation and evolution of Ginkgo in the Sichuan Basin of China

Abstract

Background

Ginkgo biloba L., an iconic living fossil, challenges traditional views of evolutionary stasis. While nuclear genomic studies have revealed population structure across China, the evolutionary patterns reflected in maternally inherited plastomes remain unclear, particularly in the Sichuan Basin - a potential glacial refugium that may have played a crucial role in Ginkgo’s persistence.

Results

Analysis of 227 complete plastomes, including 81 newly sampled individuals from the Sichuan Basin, revealed three distinct maternal lineages differing from known nuclear genome patterns. We identified 170 sequence variants and extensive RNA editing (235 sites) with a bias toward hydrophobic amino acid conversions, suggesting active molecular evolution. A previously undocumented haplotype (IIA2), predominant in western Sichuan Basin populations, showed close genetic affinity with rare refugial haplotypes. Western populations exhibited higher haplotypic diversity and distinctive genetic structure, supporting the basin’s role as both glacial refugium and corridor for population expansion. Ancient trees (314–784 years) provided evidence for interaction between natural processes and historical human dispersal in shaping current genetic patterns.

Conclusions

Our findings demonstrate substantial genetic diversity within Sichuan Basin Ginkgo populations and reveal dynamic molecular evolution through plastome variation and RNA editing patterns, challenging the notion of evolutionary stasis in this living fossil. This study provides crucial genomic resources for understanding Ginkgo’s evolution and informs conservation strategies for this endangered species.

Background

Living fossils provide exceptional opportunities for investigating speciation, extinction, and adaptive evolution [1]. These organisms, characterized by minimal morphological change over millions of years based on fossil comparisons, represent remarkable cases of evolutionary persistence [2, 3]. Iconic examples such as Ginkgo biloba L., Metasequoia glyptostroboides H.H.Hu & W.C.Cheng, coelacanths, and lungfish have survived major extinction events while maintaining their ancestral morphology [1, 4]. The evolutionary trajectory of these living fossil lineages and relict species remains contentious, with some researchers suggesting they represent evolutionary dead ends characterized by stasis and impending extinction, while others propose they retain diversification potential [5, 6]. Recent studies increasingly support active diversification within these lineages [2, 5, 7]. Moreover, morphological stasis itself may represent an adaptive strategy, allowing organisms to avoid the costs associated with phenotypic modifications [8, 9]. Thus, elucidating the evolutionary processes and patterns governing living fossils remains crucial.

Ginkgo biloba L., a tertiary relict species and the sole surviving member of the monotypic genus Ginkgo, is endemic to China [10, 11]. While historically distributed across the Northern Hemisphere until the end of the Pliocene, Ginkgo became restricted to Eastern Asia [12]. In recent times, human cultivation has expanded its range across Japan, Korea, North America, and Europe [13, 14]. Despite this widespread cultivation, Ginkgo remains classified as endangered by the IUCN due to the scarcity of natural populations [15]. Previous studies employing AFLPs, SSRs, and plastomic data have identified several potential refugia in China, including Jinfo mountain (Chongqing), Wuchuan (Guizhou), and the Tianmu Mountains (Zhejiang) [16,17,18]. Plastomic analyses across nine populations revealed significant genetic structuring among refugial populations, with western and eastern lineages diverging approximately 0.39 mya [12]. This pattern was corroborated by comprehensive whole-genome resequencing, which estimated the Southwest and East + South lineage divergence at 0.51 mya, followed by the Southern lineage separating from the East around 0.32 mya [6]. Additionally, glacial cycles have driven periodic expansion, contraction, and admixture in southwestern and southern refugia populations [6]. These findings illuminate Ginkgo’s evolutionary history, suggesting recent admixture may primarily reflect anthropogenic influences [15]. However, previous sampling efforts have not fully encompassed known wild populations.

While nuclear genomes provide numerous polymorphic loci valuable for inferring admixture and population genetics patterns due to biparental inheritance [19], their analysis is complicated by factors including paralogy and incomplete lineage sorting [20,21,22]. Plastomes, characterized by uniparental inheritance, structural conservation, absence of recombination, high copy number, and smaller genome size, present ideal markers for phylogenetic studies across taxonomic levels [23,24,25,26]. At the population level, plastome analyses have effectively revealed genetic relationships in both domesticated crops and endangered species, including tartary buckwheat, sweet potato, rice, peppers, Mendelian pea, and Cathaya argyrophylla Chun & Kuang [25, 27,28,29,30,31]. Thus, plastomic data offer promising insights into Ginkgo’s complex dispersal history.

The Sichuan Basin, previously underrepresented in sampling efforts [6, 12], warrants particular attention. This basin, with an average elevation of 500 masl, is bounded by the Hengduan, Qinling, Daba, Dalou mountain ranges, and the Yungui Plateau’s mountains, functioning as a biodiversity corridor connecting Palaearctic and Oriental biotas [32,33,34]. Influenced by the Late Pliocene Qinghai-Tibet Plateau (QTP) uplift and East Asian monsoon system (EAMS) intensification, the basin represents a phylogeographic break for various plant species, including the ring-distributed Rhodiola yunnanensis (Franch.) S.H.Fu complex [35], the relict Euptelea pleiosperma Hook.fil. & Thomson [36], and the living fossil Davidia involucrate Baill [37], which exhibits an east-west genetic division across the basin [11]. Notably, Ginkgo populations along the basin’s western edge, including Qionglai, Dujiangyan, and Pengzhou, demonstrate elevated genetic diversity compared to eastern lowland populations [38].

In this study, we analyzed 81 adult Ginkgo individuals from Chengdu City and surrounding regions in the Sichuan Basin, combined with 146 global accessions, assembling 227 plastomes representing diverse geographic and genetic lineages. Our objectives were threefold: First, we comprehensively re-evaluated Ginkgo plastome characteristics to reveal population-level variation, reassessing genetic diversity to demonstrate potential molecular adaptability despite morphological stasis in this classic “living fossil.” Second, building on previous studies, we examined genetic structure across different geographic lineages in China (East, South, Southwest, and North) to explore the origin and differentiation of maternally inherited plastomes and highlight contrasts with nuclear genome patterns. Finally, we leveraged newly sampled individuals from the western lineage to investigate diversity and differentiation within Sichuan Basin populations, analyzing potential drivers of their distribution patterns.

Results

Features of Ginkgo plastomes

Following stringent quality control, we extracted approximately 1.7 Gb of clean reads per sample for assembly, with 93% of bases achieving Q30 quality scores (detailed in Supplementary Table 3). The reads were successfully mapped and assembled into three to four contiguous sequences (contigs), achieving an average sequencing depth of 35×. Subsequently, we obtained 227 complete circular plastomes, ranging in total length from 156,910 bp to 157,047 bp (Table 1, Supplementary Tables 4–5). The Ginkgo biloba plastome displayed the characteristic quadripartite organization (Fig. 1), comprising a Large Single Copy (LSC) region of 99,179 − 99,272 bp, a Small Single Copy (SSC) region of 22,254 − 22,320 bp, and two Inverted Repeat regions (IRa and IRb), each measuring either 17,733 bp or 17,735 bp. The complete plastome exhibited an average GC content of 39.56%, with region-specific GC contents of 38.5% in the LSC, 36.23% in the SSC, and 44.63% in the IR regions.

Each plastome contained 119 unique genes, including 85 protein-coding genes (PCGs), 30 tRNA genes, and 4 rRNA genes (Supplementary Table 6). Among these, 15 genes were duplicated, comprising three PCGs (ndhB, rps7, rps12), eight tRNA genes (two with three copies), and four rRNA genes. Of the 119 unique genes, intron distribution analysis revealed that 12 genes contained single introns (nine PCGs and three tRNA genes), while three genes (clpP, rps12, and ycf3) contained two introns. Notably, the Ginkgo plastome exhibited a distinctive structural feature characterized by the loss of one ycf2 copy from the IR regions.

Fig. 1

Structural and variation analysis of 227 Ginkgo plastomes showing (i) genes on the reverse strand, (ii) genes on the forward strand, (iii) distribution of sequence variants relative to reference TM226, (iv) distribution of simple sequence repeats (SSRs), and (v) Ginkgo cone photograph (center). Three major types of variation are represented in different colors shown in (iii), corresponding to single nucleotide variants (SNV), insertiondeletions (InDel), and complex variants (combined single/multiple nucleotide polymorphisms). LSC: large single-copy region; SSC: small single-copy region; IR: inverted repeat region

Table 1 Plastome region characteristics across 227 Ginkgo individuals

RNA editing sites in Ginkgo plastome

Analysis of the TM226 plastome revealed 235 predicted RNA editing sites distributed across 59 protein-coding sequences (CDSs) (Supplementary Table 7). The vast majority (230) represented C-to-U conversions, while only five sites exhibited U-to-C modifications. All predicted editing events resulted in non-synonymous changes. Most editing events occurred at the second codon position (73.6%), with the remaining modifications (26.4%) localized to the first position. The predominant amino acid substitutions comprised serine (S) to leucine (L) (76/235), proline (P) to leucine (L) (48/235), and serine (S) to phenylalanine (F) (30/235). These modifications, primarily occurring at the second codon position, generally enhanced amino acid hydrophobicity. Among the functional gene categories, ndh genes harbored approximately one-third of all editing sites. The highest frequency of editing sites was observed in ndhF (20/235), followed by ndhD (19/235), ndhA (13/235), rpoC1 (12/235), and ndhB (11/235). Notably, seven editing sites affected initiation or termination codons through predicted conversions of CAG to UAG, CGA to UGA, and UGA to CGA, suggesting potential regulatory roles in protein synthesis.

Variation patterns in Ginkgo plastomes

Our comprehensive analysis of Ginkgo plastome variation identified 170 variants, including 120 SNVs, 37 InDels, and 13 complex variants (Table 2). The distribution of these variants is visualized in the circular genome map (Fig. 1, detailed in Supplementary Tables 8–10), where InDels frequently co-localize with SSR regions. The Large Single Copy (LSC) region harbored the majority of variants (125), followed by the Small Single Copy (SSC) region (42), while only three variants were detected in a single Inverted Repeat (IR) region. However, when normalized by region length, the SSC region exhibited the highest mutation density (1.88 variants/kb), followed by the LSC region (1.26 variants/kb), with the IR region showing markedly lower mutation frequency (0.17 variants/kb) (Table 2).

SNVs, the predominant form of variation, showed a distinct mutational bias toward GC to AT substitutions, specifically G to T and C to A transitions (Fig. 2a, Supplementary Table 11). Within coding sequences (CDS), we identified 48 SNVs without other variant types, underscoring the conservative nature of the Ginkgo plastome. Notably, 75% of these SNVs represented nonsynonymous mutations (Table 3). The highest frequency of nonsynonymous mutations occurred in ycf2, followed by ycf1. These nonsynonymous substitutions predominantly resulted in amino acid changes from phenylalanine, isoleucine, and proline to the more hydrophobic leucine across the Ginkgo population (Fig. 2b, Supplementary Table 12), suggesting potential adaptive significance.

Analysis of mutation frequencies across genomic segments (Fig. 2c, Supplementary Table 13) corroborated the overall conservative pattern, with most regions displaying low mutation rates. However, specific genes and non-coding regions within the SSC and LSC regions exhibited elevated mutation frequencies. Notably, the intergenic spacer regions IGS_psbK-psbI and IGS_psaB-psaA showed higher mutation rates, potentially indicating reduced functional constraints in these regions, thereby permitting greater mutational accumulation.

Table 2 Distribution of sequence variants across Ginkgo plastome regions

Table 3 SNV distribution and impact in protein-coding regions

Fig. 2

Nucleotide variation patterns and their functional implications in Ginkgo plastomes. (a) Base substitution patterns in SNVs. (b) Resulting amino acid changes. (c) SNV frequency distribution across genomic regions (CDS, introns, and intergenic regions). Regions with above-average SNV frequencies are labeled

Genetic diversity and differentiation of Ginkgo plastomes

The geographic lineage classifications in this study were based on previous genomic analyses by Zhao et al. (2019), which investigated the evolutionary relationships between nuclear and organellar genomes. We calculated nucleotide diversity (Pi) and genetic differentiation (pairwise FST) across geographic lineages using plastome sequences (Table 4). Among all populations, newly sampled individuals from the Chengdu region exhibited the highest nucleotide diversity (Pi = 0.21 × 10⁻³), followed by the North lineage (Pi = 0.19 × 10⁻³). In contrast, both Overseas and South lineages displayed markedly lower genetic diversity (Pi = 0.03 × 10⁻³).

Population differentiation analyses revealed distinct patterns among lineages. The Southwest and North lineages showed minimal genetic differentiation, with FST values approaching zero. Similarly, low differentiation levels were observed among the Chengdu, Southwest, and North lineages (FST = 0.0358–0.0434), as well as between the East and South lineages (FST = 0.051). However, substantial genetic differentiation was detected when comparing the Chengdu, Southwest, and North lineages with the East and South lineages (FST = 0.1831–0.31). The Overseas lineage exhibited moderate levels of differentiation from the East, North, and South lineages (FST = 0.145, 0.154, and 0.154, respectively), with slightly elevated differentiation from the Southwest lineage (FST = 0.181).

Table 4 Genetic diversity index and population differentiation among geographic lineages

Genetic structure of Ginkgo populations

To elucidate the population genetic structure among Ginkgo individuals, we conducted a principal component analysis (PCA) based on plastome alignments (Fig. 3b). The first two principal components collectively explained 93.83% of the total variance, with PC1 accounting for 84.51% and PC2 for 9.32%. Notably, PC1 segregated the samples into three distinct clusters, designated as Groups I-III (comprehensive statistical values are presented in Supplementary Table 14). The substantial overlap of data points within each group indicated high genetic homogeneity among Ginkgo samples within clusters.

To further resolve the genetic relationships, we reconstructed maximum likelihood (ML) and Neighbor-Joining (NJ) phylogenetic trees using plastome alignments (Fig. 3c and Supplementary Fig. 1). Both resulting topologies corroborated the PCA findings, revealing three primary clades (Groups I, II, and III). The significant variation in branch lengths among these groups indicated substantial genetic differentiation. Group I comprised a small cluster of only three samples, including two from the East lineage and one from the North lineage (ML bootstrap support, MLBS = 100%). Group II predominantly consisted of specimens from the Southwest and North lineages, incorporating newly sequenced individuals from the Chengdu population (MLBS = 100%). Terminal branches within this group primarily represented populations from montane regions, including Mt. Jinfo (Chongqing), Mt. Huping (Hunan), and Mt. Shennongjia (Hubei). Group III emerged as the most extensive and diverse clade, encompassing samples from multiple nuclear lineages. The characteristically short branch lengths within this group suggested recent diversification or extensive gene flow among populations, likely reflecting widespread historical human-mediated transplantation and cultivation activities.

Fig. 3

Population structure and geographic distribution of Ginkgo plastomes. (a) Distribution of 209 Ginkgo samples across 46 Chinese populations. Circle size represents sample number and colors indicate geographic grouping. Green and blue dashed lines indicate the geographic regions detailed in Fig. 4c and d, respectively. (b) Principal component analysis showing genetic clustering of 227 plastomes, with percentages of explained variance indicated for PC1 and PC2. (c) A maximum likelihood phylogenetic tree based on complete plastome sequences (n = 227). Bootstrap support values > 90% are shown for major branches. Colors indicate geographic lineages; outer rings denote haplotype assignments

Haplotype network

To elucidate the maternal origins of the Ginkgo plastome, we conducted a comprehensive haplotype analysis across our sampled populations. Our analysis revealed 16 distinct haplotypes (detailed in Supplementary Tables 15–16). Network analysis delineated three major haplogroups, each distinguished by more than ten informative polymorphic sites (Fig. 4a). Haplogroup I comprised two low-frequency haplotypes, represented by one sample from the North lineage and two from the East lineage, including a specimen from Tianmu Mountain (TM226), which is considered an ancestral population. Within Haplogroups II and III, haplotypes IIA2 (n = 42, accounting for 19% of total haplotypes) and IIIB1 (n = 60, 26%) emerged as the predominant variants, respectively. Notably, IIIB1 exhibited the most extensive geographical distribution, spanning all six sampling regions and co-occurring with several other prevalent, phylogenetically related haplotypes, including IIIC7 (n = 38, 16%), IIID3 (n = 33, 15%), and IIIB4 (n = 27, 12%).

Distinct phylogeographic patterns emerged across different lineages. The East lineage harbored five plastomic haplotypes, with IIIB1 (16/26) predominating, followed by IIIB4 (5/26). The South lineage exhibited six haplotypes, dominated by IIIB4 (11/30), IIIB1 (8/30), and IIID3 (7/30). The close phylogenetic relationship between haplotypes IIIB1 and IIIB4, separated by a single mutational step, suggests a close relationship between East and South lineage. The Southwest lineage was characterized by a high frequency of the IIIC7 haplotype (10/26), accompanied by IIIB4 (5/26), IID11 (4/26), IIA2 (4/26), and IIIB1 (3/26). The North lineage demonstrated the highest haplotypic diversity, encompassing nine distinct haplotypes, with IIIB1 (16/46), IIA2 (6/46), IIIC7 (6/46), and IIIB4 (4/46) being most prevalent. Among the overseas populations (18 samples), only three distinct haplotypes were identified, all belonging to haplogroup III (IIIB1, IIIC7, IIID3). The predominant haplotype IIIC7 (western dominant haplotype) was primarily distributed throughout South Korea and the Japanese islands of Kyushu and Shikoku. In contrast, all six samples from Europe and the United States carried haplotype IIIB1, which was also present in South Korea. Additionally, three samples from Honshu, Japan, possessed the IIID3 haplotype. This distinct haplotype distribution pattern suggests two independent cultivation origins for overseas populations, indicating multiple historical introduction events in the spread of Ginkgo beyond its native range.

Fig. 4

Haplotype networks and geographic distribution patterns in Ginkgo. (a) Haplotype network of Ginkgo plastome, with mutational steps indicated on connecting branches. (b) Haplotype number distribution across geographic lineages (colors consistent with Fig. 3). (c) Spatial distribution of haplotypes in western populations (North, Southwest, and Chengdu, corresponding to the green dashed line area in Fig. 3a), with Sichuan Basin boundary indicated by gray dashed line. (d) Detailed haplotype composition of Chengdu sub-populations (corresponding to the blue dashed outline in Fig. 3a), with sample sizes in parentheses

Haplotype origin of Sichuan Basin populations

To investigate the origins of haplotypes in the newly sampled Chengdu individuals (including those from Chengdu city and surrounding areas), relevant haplotypes were selected for comparative analysis (Fig. 4b-c). Nine distinct haplotypes were identified across the 18 Chengdu populations, with haplotype IIA2 being the most prevalent (n = 32, 39.5%), followed by IIID3 (n = 18, 22.2%) and IIIC7 (n = 11, 13.6%). The geographic distribution of these haplotypes varied notably across Chengdu. Haplotype IIA2 was most common in eastern sub-populations, including Shuangliu, Jinjiang, Jianyang, Chenghua, Qingbaijiang, Jintang, and Xindu. In contrast, all nine haplotypes were present in the western sub-populations, with IIA2 occurring in seven and IIID3 in six of these sub-populations.

The spatial distribution of haplotypes and the diameter at breast height (DBH) of sampled trees offer further insights (Supplementary Table 17). Among the sampled individuals, ten had a DBH ranging from 134 to 350 cm, corresponding to approximate ages of 314 to 784 years. Of these, four exhibited haplotype IIIB1, five displayed IIA2, and one had IIIC7. Notably, two IIA2 samples were obtained from a nature reserve in Dayi County (western region), suggesting a potential natural origin. These three predominant haplotypes were also shared with populations in the Southwest and North lineages. Haplotype IIA2, while separated from IIIB1 and IIIC7 by over 60 mutational steps, showed closer genetic affinity to low-frequency haplotypes such as IID11 from the Jinfo population (a refugium), IIC10 from the Huping population (potentail refugium), and IIB14 from the Shennongjia population (potential refugium). The high frequency of IIA2 within Chengdu, along with its close relationship to the Chengdu-specific haplotype IIA5, suggests that the Sichuan Basin could have been a center of origin for this lineage of plastomes. Furthermore, complex admixture events may have taken place along the western edge of the Sichuan Basin, potentially facilitating the integration of IIA2 haplotypes into the broader Ginkgo gene pool.

Discussion

Plastomic variation and morphological stasis in Ginkgo

Traditional perspectives on gymnosperm evolution suggest their lower diversification rates compared to angiosperms may reflect higher extinction rates rather than reduced evolutionary potential [39, 40]. While extensive early gymnosperm diversification events may be obscured by extinction [41], evidence from cycads demonstrates recent radiations within the last 12 million years, indicating that early-diverging lineages retain evolutionary capability [2]. Although similar radiation patterns may have occurred in Ginkgo, the absence of extant Ginkgoales lineages complicates diversification analyses. Nevertheless, recent evidence challenges the notion of Ginkgo as an evolutionary dead end, highlighting its adaptations to environmental stresses, substantial genetic diversity, and dynamic population responses to glacial cycles [6]. Our plastome analyses provide additional support for this perspective.

Our investigation revealed considerable genetic variation in Ginkgo plastomes, with 170 variants (120 SNVs, 37 InDels, and 13 complex variants, Table 2) across 227 samples yielding 16 distinct haplotypes. The distribution of these variants showed clear regional patterns, with the SSC region exhibiting the highest mutation density (1.88 variants/kb), followed by the LSC region (1.26 variants/kb), while the IR region displayed notably lower mutation frequency (0.17 variants/kb). This moderate nucleotide diversity (Pi = 0.18 × 10⁻³) is comparable to other tree species, such as Macadamia integrifolia Maiden & Betche (407 SNVs, 38 haplotypes, 63 individuals), Pistacia chinensis Bunge (460 SNVs, 12 haplotypes, 39 individuals, Pi = 0.82 × 10⁻³), Bretschneidera sinensis Hemsl. (105 polymorphic sites, 55 InDels, 12 individuals, Pi = 0.26 × 10⁻³), and Phoenix dactylifera L. (37 SNVs, 4 haplotypes, 201 individuals) [42,43,44,45]. Particularly noteworthy is the pattern of sequence variation within coding regions, where 75% of the 48 identified SNVs represented nonsynonymous mutations, primarily in ycf2 and ycf1 (Table 3). These changes showed a distinct bias toward GC to AT substitutions and frequently resulted in conversions to more hydrophobic amino acids, suggesting potential adaptive significance. The elevated mutation rates in specific intergenic regions (IGS_psbK-psbI and IGS_psaB-psaA) likely reflect reduced functional constraints, allowing for greater evolutionary flexibility in these regions (Fig. 2).

An intriguing pattern emerged when comparing sequence-level variations with RNA editing events in the Ginkgo plastome. At the DNA sequence level, we observed a predominance of nonsynonymous mutations leading to hydrophobic amino acids, particularly leucine (L) (Fig. 2b). However, RNA editing events showed an opposing pattern. Notably, over half of these editing events result in conversions to hydrophobic amino acids, particularly leucine and phenylalanine (Supplementary Table 7). This seemingly contradictory pattern might reflect different evolutionary strategies operating at distinct molecular levels. The bias toward hydrophobic amino acid substitutions in DNA sequence variations, particularly in genes like ycf2 and ycf1, likely represents long-term adaptive evolution, as these changes might have been fixed in the genome. These substitutions may enhance protein stability or membrane interactions, which are crucial for plastid function [46]. In contrast, the RNA editing system provides a more flexible, reversible mechanism for protein modification [47, 48]. The prevalence of editing events converting leucine to hydrophilic amino acids suggests a dynamic regulatory system that can fine-tune protein properties in response to environmental conditions or developmental stages [49]. This dual-level regulation - stable hydrophobic modifications at the DNA level combined with flexible RNA editing - may represent a sophisticated evolutionary strategy. It allows Ginkgo to maintain core protein functions through conserved hydrophobic domains while retaining the ability to modulate protein properties through post-transcriptional modifications. Such molecular plasticity could explain how Ginkgo has maintained its ability to adapt to environmental changes despite apparent morphological stasis.

Recent studies propose that morphological stasis itself may represent an adaptive strategy [8, 9, 50]. Under this paradigm, molecular evolution may proceed independently of external morphology, or actively maintain morphological stability. These findings contribute to growing challenges against the traditional “living fossil” concept as an evolutionary model [4, 51,52,53]. Based on our analyses here, Ginkgo is still evolving at the molecular level while maintain similar morphology to previous fossil lineage. Our analyses support active molecular evolution in Ginkgo despite its conserved morphology, suggesting a more complex evolutionary trajectory than previously recognized.

Population differentiation and maternal lineages in Ginkgo

Our plastome-based analysis reveals distinct population structure patterns from those previously reported using nuclear genomic data, providing new insights into Ginkgo’s evolutionary history. While nuclear genome studies identified four major geographic lineages (East, South, Southwest, and North) with evidence of extensive gene flow [6], our plastome analysis reveals three distinct maternal lineages (Groups I-III, Figs. 3b-c and 4a) that transcend geographic boundaries. This discordance between nuclear and plastome patterns likely reflects the different inheritance modes and evolutionary histories of these genomic components. The maternal lineages exhibit several notable features. Group I, comprising only three samples from the East and North lineages (including the ancestral Tianmu Mountain population), represents a rare, potentially relict maternal lineage. Group II emerged as a particularly interesting clade, predominantly consisting of specimens from the Southwest and North lineages, including our newly sequenced Chengdu populations. This group is characterized by well-supported terminal branches (MLBS = 100%) representing montane populations from Mt. Jinfo, Mt. Huping, and Mt. Shennongjia, suggesting the preservation of distinct maternal lineages in these mountainous regions. Group III, the most extensive and diverse clade, shows extensive haplotype sharing across geographic lineages, with the widespread haplotype IIIB1 (26% of total) occurring across all sampling regions.

West-East population differentiation pattern [12, 18] is partially supported by our findings, with Chengdu, Southwest, and North lineages showing high differentiation from East and South lineages. However, our maternal lineage analysis reveals more complex patterns than previously recognized. While nuclear data suggested the North lineage originated through admixture between Southwest (28.45%) and South (71.55%) lineages approximately 139,260 years ago [6], our plastome analysis suggests a stronger maternal contribution from western populations. The eastern region of the North lineage carries haplotype IIIB1 (suggesting East lineage ancestry), while western populations exhibit haplotypes IIA2, IIIC7, and IIIB4 (indicating genetic contributions from Chengdu, Southwest, and South lineages respectively). Future studies incorporating additional samples from northern populations may further elucidate these complex patterns of maternal inheritance and better resolve the apparent discordance between nuclear and plastome-based population histories.

The global distribution of cultivated Ginkgo populations [13, 14, 54] shows distinct maternal inheritance patterns (Fig. 3c). European and American specimens exclusively carry the widespread haplotype IIIB1 and show minimal differentiation from the East lineage, supporting previous hypotheses of multiple introductions from eastern China [6]. Notably, Japanese and South Korean populations predominantly carry haplotype IIIC7, prevalent in the Southwest lineage, suggesting this lineage as their primary maternal source. This pattern indicates multiple independent cultivation origins for overseas populations, highlighting the complex history of human-mediated Ginkgo dispersal.

This intricate pattern of maternal inheritance, characterized by three distinct lineages with varying degrees of geographic overlap and haplotype sharing, likely reflects the interplay between natural population processes and historical human-mediated dispersal. The discordance between nuclear and plastome patterns suggests different evolutionary trajectories for maternal and paternal lineages, potentially influenced by specific dispersal mechanisms, selective pressures, and human cultivation practices.

Genetic diversity of Ginkgo in Sichuan Basin

Our analysis of plastome diversity in the Sichuan Basin reveals complex patterns of genetic structure and potential refugial history. Haplogroup II, predominantly comprising Chengdu samples, showed distinctive features including long terminal branches associated with known refugial areas such as Mt. Jinfo, Mt. Huping, and Mt. Shennongjia in western and central China (Fig. 3c, Supplementary Fig. 1). These mountains have been previously identified as refugia for several endangered tertiary relict species, including D. involucrata and Tetracentron sinense Oliv [37, 55, 56]. The presence of unique haplotypes in these refugial populations (IID11 in Jinfo, IIC10 in Huping, and IIB14 in Shennongjia) suggests limited gene flow and potentially independent evolutionary trajectories, highlighting their conservation significance. The Sichuan Basin and its surrounding forests have historically served as glacial refugia for numerous species [32, 35, 57]. Climate simulations based on bird distribution data indicate that the basin, particularly its lowland areas, may have provided crucial warm refugia for species displaced by decreasing temperatures during glacial periods [32]. Our analysis revealed a distinct east-west differentiation pattern within the basin (Fig. 4c-d). Haplotype IIA2, while widely distributed throughout the region, showed higher frequency in western populations and demonstrated close genetic relationships with rare refugial haplotypes from southwestern and central China. The prevalence of IIA2 in Sichuan, coupled with the discovery of the closely related Chengdu-specific haplotype IIA5, suggests this region may represent the center of origin for this lineage.

The western edge of the Sichuan Basin exhibited notably higher haplotypic diversity than eastern regions, containing all nine identified haplotypes (Fig. 4). This pattern aligns with the area’s proposed role as a biodiversity corridor connecting Palaearctic and Oriental biotas [32]. The elevated diversity in western populations may reflect postglacial range expansion from southwest to north along mountain corridors. Similar east-west genetic structuring across the Sichuan Basin has been documented in other tertiary relicts like T. sinense [55], E. pleiosperma [36], and D. involucrate [37] though such pronounced within-basin differentiation is unusual.

The age structure of sampled trees provides additional insights, with ten individuals (DBH 134–350 cm, estimated ages 314–784 years) representing three major haplotypes (IIIB1, IIA2, and IIIC7). Notably, two IIA2 samples from a western nature reserve in Dayi County suggest potential natural origins. The complex distribution of these ancient trees and their haplotypes indicates interplay between natural population processes and human cultivation practices. Similar to other relict species such as cycads and Metasequoia glyptostroboides, Ginkgo has successfully expanded its distribution through human-mediated dispersal [6, 58, 59], offering valuable insights for conservation strategies.

While our findings illuminate the complex evolutionary history of Ginkgo in the Sichuan Basin, several questions remain unanswered. The mechanisms driving specific maternal haplotype distributions and the species’ adaptation to environmental heterogeneity along the basin’s western edge during Pleistocene glacial cycles warrant further investigation. Future research should focus on disentangling the adaptive evolutionary mechanisms underlying these environmental transitions and the relative contributions of natural and anthropogenic factors in shaping current population structure.

Conclusion

Our analysis of 227 complete plastomes, including 81 newly sequenced individuals from the Chengdu region, revealed significant genetic diversity in Ginkgo populations. The plastome exhibited notable adaptive potential through abundant RNA editing sites that predominantly convert amino acids to hydrophobic forms. Importantly, our sampling in the western Sichuan Basin uncovered previously undocumented genetic variation, including novel haplotypes that may represent ancestral refugial populations. The western edge of the basin appears to function as a biodiversity corridor, facilitating population connectivity and contributing to an unexpected east-west genetic differentiation pattern. This diversity pattern likely reflects both natural dispersal processes and historical human-mediated transplantation. These findings challenge the notion of evolutionary stasis in this living fossil and highlight the importance of the Sichuan Basin populations for conservation efforts, particularly those harboring rare haplotypes essential for the species’ long-term persistence.

Methods

Sample collection and data acquisition

We collected leaf samples from 81 adult Ginkgo individuals across 18 regions in Chengdu, China (Supplementary Table 1). The plant material used in this study was formally identified by Dr. Shenglong Kan from Shandong University. Tree circumference and location data were recorded for each individual, and ages were estimated based on diameter at breast height (DBH) following Tang et al. [60]. Total genomic DNA was extracted using the CTAB method [61] and sequenced on an MGI DNBSEQ-T7 platform (MGI-TECH, Shenzhen, China) using 150 bp paired-end libraries. We supplemented our dataset with whole-genome sequences of 146 additional individuals (50 bp paired-end reads generated on BGISEQ-500) obtained from CNGBdb (study accession No. CNP0000136, Zhao et al., 2019; Supplementary Table 2). The final dataset comprised 227 samples: 209 from 46 Chinese populations and 18 from overseas populations in Japan, South Korea, North America, and Europe. Geographic distributions were visualized using R 4.3.1 packages sf, ggspatial, and ggplot2, with geo.json files from official sources (http://xzqh.mca.gov.cn/map, GS (2022)1873).

Plastome assembly and annotation

Raw sequence quality was assessed using FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) with default settings. Downloaded sequences were filtered using Trimmomatic v0.39 with settings of HEADCROP:3 MINLEN:45 AVGQUAL:30 [62]. Clean reads were aligned to the Ginkgo reference plastome (NC_016986) using BWA v0.7.17 [63] and processed with SAMtools v1.7 [64] to facilitate individual plastome assembly. Complete plastomes were assembled using SPAdes v3.15.5 [65]. The resulting assembly graph after simplification retained full assembly path information. To obtain a single circular plastome, the reference plastome was aligned to the assembly graph, followed by de-looping using Bandage v0.81 [66]. Newly assembled plastomes were subsequently annotated with GeSeq [67] and visualized using OGDRAW [68].

Simple sequence repeats (SSRs) were identified using MISA (http://pgrc.ipk-gatersleben.de/misa/) with standard parameters (minimum repeat units: mononucleotide-10, dinucleotide-6, tri/tetra/penta/hexanucleotide-5). RNA editing sites (C-to-U and U-to-C) were predicted using PREPACT 3.0 with default settings [69].

Variant analysis

To analyze plastome variant distribution across all Ginkgo samples, we identified and annotated sequence variations using snippy software (https://github.com/tseemann/snippy), including Single Nucleotide Variants (SNVs), Insertions/Deletions (InDels), and complex variants (comprising both single and multiple nucleotide polymorphisms) based on complete plastome assemblies. To prevent bias from duplicated sequences, we excluded one IR region from the analysis. The reference sequence (TM226) used for variant detection underwent thorough annotation verification to ensure the absence of gaps or ambiguous positions (N), thereby guaranteeing accurate variant annotation. Following variant identification, we systematically categorized mutations by genomic region (CDS, intergenic, intron) and normalized frequencies by region length to enable accurate comparisons. To avoid potential bias from short segments in frequency calculations, only regions exceeding 200 bp were included in the analysis. The identified variant positions were subsequently used for downstream analyses. Results were visualized using R packages circlize [70] and ggplot2 [71].

Population structure analysis

The 227 plastomes were aligned using MAFFT v7.505 with default parameters [72]., yielding a 157,118-site alignment. Maximum likelihood phylogenetic analysis was performed using IQ-TREE2 in PhyloSuite with settings -st DNA -m K81u + I + F -bb 1000 -bcor 0.90 [73, 74]. The consensus tree was visualized using TVBOT [75]. We also constructed a Neighbor-Joining (NJ) tree with 10,000 bootstrap replicates in MEGA7 to assess clade support [76,77,78]. To further investigate population structure, we performed a Principal Components Analysis (PCA) with the glPca function in the adegenet package in R 4.3.1 [79].

Haplotype and genetic diversity analysis

Population genetic parameters including haplotype diversity, nucleotide diversity (Pi), and pairwise FST were calculated using DnaSP6 v6.12.03 [80]. Frequencies of different haplotypes in each population were calculated in Arlequin [81]. Haplotype sequences and frequencies were used to infer haplotype networks in POPART v1.7 with the integer NJ network method [82].