We collected 398 individuals from 47 populations across North America. Our sampling extended as far north as the population in Ontario, Canada (45.1749°N and 74.8326°W), as far east as the population in Pennsylvania, USA (41.7364°N and 70.8566°W), as far south as the population in Texas, USA (29.8031°N and 98.4934°W), and as far west as the population in Texas, USA (30.7052°N and 104.2134°W), basically covering their entire distribution range from south to north of eastern North America. Three to 12 individuals per population were randomly sampled intervals were ≥ 10 m apart. After species identification, the leaflets of fresh healthy leaves were dried in silica gel and the dried leaf tissue samples were stored at -20°C for further extraction of genomic DNA. The geographic location of the populations, the number of samples and voucher information are shown in Table 1 .

DNA was extracted from silica-dried leaves using a modified CTAB method (Doyle and Doyle, 1987) or using the DNeasy Plant Mini Kit (Qiagen, Crawley, UK). Amplification and sequencing followed Soejima and Wen (2006) for the plastid sequences (trnL-F, rps16 and trnC-petN), and Ehrenreich and Purugganan (2008) for the nuclear ARF6 gene. DNA sequences were assembled using Sequencher v4.1.4 (Gene Codes Corp., Ann Arbor, Michigan, USA). All sequences obtained were aligned using MUSCLE v3.8 (Edgar, 2004) and the alignment was then adjusted manually.

Shared haplotypes were determined using DnaSP v5.10 (Librado and Rozas, 2009). The number of haplotypes (H) and polymorphic sites (S), haplotype diversity (H_d), and nucleotide diversity (Pi) were calculated using DnaSP. We constructed the network relationships with cpDNA and nrDNA haplotypes using PopART v1.7 with Median-Joining model, respectively (Leigh and Bryant., 2015). Population structure and relationships among haplotypes were conducted using maximum parsimony network in PAUP v4.0 (Swofford, 2002), which was designed to construct the shortest, least complex network. In this analysis, gaps with two or more base pairs were coded as single mutation events. When overlapping indels occurred, the overlap portion was considered to be a single event.

The software STRUCTURE v2.3.4 was used to analyze the genetic structure of the populations. Clustering method based on Bayesian model was used to assign genotypes/individuals to different clusters according to shared co-ancestry to describe the genetic structure well (Pritchard et al., 2000). The software was run using the Admixture Model with parameters set to 20,000 burn-in repeats and 70,000 MCMC repeats. The number of clusters (K) was set to vary from two to 12. For each value of K, we performed was 20 runs. The relationship between K and LNP (D) and ΔK calculated by Structure Harvester was used to obtain the best K value.

An analysis of molecular variance (AMOVA) was used to partition genetic variation among and within groups, as implemented in ARLEQUIN v3.5 (Excoffier and Lischer, 2010). G_ST and N_ST among populations were calculated from the chloroplast markers using 1000 permutations in PermutCpSSR v2.0 (Burban et al., 2010). The principal coordinate analysis (PCoA, Peakall and Smouse, 2006) was carried out using DARwin v7.0 software (Perrier and Jacquemoud-Collet, 2006). Nei’s genetic distance among the populations of Parthenocissus was calculated in MEGA v7.0, and Neighbor-Joining (NJ) trees were constructed for 47 population using the Nei’s genetic distance (Kumar et al., 2016; Nei, 1972).

We examined pairwise mismatch distributions based on the pairwise nucleotide differences between haplotypes to detect demographic expansions using ARLEQUIN. Populations at demographic equilibrium should present a multimodal or random and rough distribution of pairwise differences, whereas populations experiencing a sudden demographic expansion are expected to display a unimodal and smooth distribution (Slatkin and Hudson, 1991). In order to test whether the overall distribution area and different populations of Parthenocissus had expanded historically, we conducted neutrality tests by DnaSP based on cpDNA and nrDNA.

The divergence time estimate was conducted in BEAST 1.8.4 (Drummond and Rambaut, 2007) using the dataset including 47 populations with Parthenocissus chinensis as outgroup. The dating dataset was partitioned using BEAUti 1.8.4 to generate input files for BEAST. Under the Akaike information criterion (AIC) implemented in MrModeltest 2.3, the best-fit model of nucleotide substitution for this analysis was determined to be HKY (Posada and Crandall, 1998). We applied the HKY model under an uncorrelated lognormal relaxed clock model (Drummond et al., 2006). MCMC analyses of 100,000,000 generations were implemented, in which every 1,000 generations were sampled. The first 10% of generations were discarded as burn-in, and the parameters were checked using the program Tracer 1.6, when the effective sample size of all parameters exceeds 200, the results were considered reliable. The rest sampled posterior trees were summarized to generate a maximum clade credibility tree using the program TreeAnnotator 1.8.4 (Drummond and Rambaut, 2007). The program Figuretree 1.4 (Drummond and Rambaut, 2007) was used to compile and visualize the results from BEAST. According to Rogers and Harpending (1992), the formula T=τ/2μkg was applied to calculate the population expansion time (τ: expansion parameter from mismatch distribution analysis; μ: mutation rate; k: average sequence length of the cpDNA region under study, the value is 1223bp, see the Results section; g: Generation time of Parthenocissus, calculated in 3 years, Li, 1998).

The total alignment of the three chloroplast regions (trnL-F, rps16 and trnC-petN) surveyed across all the individuals was 1223 bp, and 38 polymorphic sites were observed, all of which were indels. A total of 47 chloroplast haplotypes (C1-C47) were identified ( Table 1 ). The nuclear ARF6 gene was 488 bp long with 16 polymorphic sites, including 1 indel and 15 base substitutions. We identified 29 nuclear haplotypes (H1-29) across the 47 surveyed populations ( Table 1 ).At the species level, the cpDNA data showed a higher estimates of haplotype diversity (Hd = 0.9276) than the value from the nrDNA data (Hd = 0.8331). However, the nucleotide diversity of nrDNA (pi= 4.06×10^-3;) was higher than that of cpDNA (pi = 3.99×^-3) ( Table 2 ). The C2 and H2 haplotypes were most common in Parthenocissus, with a frequency of 16.3% (65 accessions and 15 populations) and 68.2% (223 accessions and 39 populations), respectively. The distribution frequency of C2 and H2 was high ( Table 1 ), basically located in the center of their branches.

The parsimony network grouped the 47 cpDNA haplotypes into two major clusters (Cluster A and Cluster B) separated by 4 mutational steps ( Figure 1A ). Thus, each region mostly harbored a genealogically distinct set of haplotypes. The cluster A included 32 cpDNA haplotypes. These cpDNA haplotypes were found quite broadly from the Southern North America region (SA, 23 unique haplotypes and 7 shared haplotypes). Note that C1, C2 and C16 were central haplotypes from which other haplotypes diverged. There were 11 haplotypes in cluster B diverged from center on C19 and C24. Meanwhile, a haplotype network was constructed based on 29 nrDNA haplotypes ( Figure 1C ). The plausible network tree of nrDNA had three clusters centering on H1, H2 and H15, and other haplotypes diverged from these three centers. The phylogenetic trees resulting from MP of cpDNA and nrDNA ( Figures 1B, D ) supported a similar pattern observed in the network analysis.

In the case of the STRUCTURE analysis, we detected three phylogeographic groups (F_CT = 0.48124, p<0.001) as the optimal number of genetic “groups” (K) based on spatial locations and cpDNA haplotypes ( Table 3 , Figure 2A ). Interestingly, these phenomena appeared in nrDNA but were not very obvious ( Figure 2B ). The southern North American (SA) group was the largest grouping with 22 populations (SA1-22), which were distributed across Texas, Arkansas, Louisiana, Mississippi, Alabama, Tennessee, Georgia, North Carolina, Florida and Virginia. The eastern North American (EA) group included 13 populations from Virginia, Connecticut, New York and Pennsylvania (EA23-30 and EA41-45). The northern North American (NA) group was exclusively located in Ontario, Quebec and Massachusetts (NA31-40 and NA46-47), and this assemblage contained 12 populations ( Figure 2 ). We also conducted haplotypes and STRUCTURE analyses for P. quinquefolia and P. vitacea - P. heptaphylla separately ( Supplementary Figures S1 , S2 ). The results showed that no clear biogeographic pattern was found in P. quinquefolia, but P. heptaphylla and P. vitacea could be divided into two groups ( Supplementary Figure S2 ), the southern group mainly occurred on the Edwards Plateau in Central Texas, and the northern group consisted of P. vitacea in Canada and the northern USA ( Supplementary Figure S1 , Supplementary Table S1 ).

A permutation test showed that N_{ST =} 0.488 was significantly greater than G_ST (0. 0.444, P < 0.05) in cpDNA, and that N_{ST =} 0.457 was significantly greater than G_ST (0. 0.306, P < 0.05) in nrDNA ( Table 2 ). In terms of AMOVA results based on cpDNA, approximately 66.45% of total variation was explained by differences within populations and 33.55% due to differences among populations. For nrDNA, 92.78% of variation was partitioned within populations and 7.22% among populations ( Table 3 ). For the three recognized groups in STRUCTURE, approximately 69.11% of variations occurred within populations, and the remaining 29.40% and 1.49% occurred among populations and among populations within groups in cpDNA, respectively. Meanwhile, 92.06% of variation was partitioned within populations, and the remaining 5.96% and 1.98% among populations and among populations within groups in nrDNA, respectively ( Table 3 ).

The UPGMA clustering tree constructed based on Nei’s genetic distance indicated that nrDNA showed more complex pattern than cpDNA ( Figure 3 ). Principal coordinate analysis based on similarity matrix also agreed with UPGMA structure and STRUTURE results ( Figure 4 ; Supplementary Figure S2 ).

Our cpDNA results showed that there are no explicit signals of population expansion or equilibrium in neutrality tests. The observed mismatch distribution of EA regions did not reject the spatial expansion model ( Table 2 ), but a unimodal distribution was not identified in all regions ( Figure 5A ). Estimates of Tajima’s were generally nonsignificant for all nrDNA regions of Parthenocissus ( Table 2 ). By contrast, the mismatch distribution in nrDNA showed that the overall population was unimodal ( Figure 5B ), closely fitted to the expected distribution under the sudden expansion model. The Sum of Squared deviation (SSD) of 0.00258 (p>0.05) and Harpending’s Raggedness index (H_Rag ) of 0.02368 (p>0.05; Table 2 ) could not reject the population expansion model.

The BEAST analyses based on two calibration points suggested an origin of the North American Parthenocissus crown lineage at 8.25 Ma with a 95% HPD of 6.55-10.03 Ma based on the combined cpDNA ( Figure 6A ). The crown node age of Parthenocissus was estimated to be 7.95 Ma with a 95% HPD of 6.21-9.67 Ma in nrDNA data ( Figure 6B ). The expansion time was estimated to be 0.074-0.604 Ma.

Our results of cpDNA and nrDNA haplotypes demonstrate a high level of genetic diversity across the 47 populations of Parthenocissus in North America ( Table 2 ). A possible explanation for the high diversity in these species could be its long evolutionary history, which may have allowed accumulation of genetic variation. There are some characteristics in Parthenocissus such as hermaphroditism, attachment to various trees, and fruits that attract birds for seed spread (Wen, 2007; Moran et al., 2009; Tiffney and Barghoorn, 1976), which have led these species to adapt and evolve under diverse habitats. Williams et al. (2004) utilized a fossil-based data set from over 700 sites in northern and eastern North America to review the late-Quaternary vegetation history of this region at different ecological organizational levels, from individual taxa to biomes. They found that during the full-glacial period (21,000 - 17,000 yr ago [calendar years]) and the mid- to late-Holocene (7,000 - 500 yr ago), the distribution and composition of the vegetation were relatively stable. However, rapid changes occurred during the late glacial- and early-Holocene (16,000 - 8,000 yr ago) and after 500 yr ago. Besides the northward redistribution of most taxa, large - scale east - west distribution shifts were also observed. The wide geographic ranges of those species across North America have provided ample opportunity for isolation, drift and mutation (Wang et al., 2009). This finding is not surprising, as there are many species with high genetic diversity in North America, such as Trillium grandiflorum, Smilax hispida, Smilax rotundifolia and Picea glauca. We found that they have similar dispersal mechanisms. For example, T. grandiflorum spreads its pollen through bumblebees and its seeds through white-tailed deer (Griffin and Barrett, 2004). Smilax has small fleshy fruits, and its seeds are dispersed by birds while its pollen is spread by insects (Zhao et al., 2013). Additionally, Picea glauca has its seeds dispersed by birds and its pollen dispersed by the wind (O’Connell et al., 2007). Moreover, Nadeau et al. (2015) suggested that high levels of genetic diversity were maintained across the range of Pinus strobus, likely via frequent long-distance dispersal events during colonization. O’Connell et al. (2007) also pointed out in the study on white spruce that extensive long-distance, pollen-mediated gene flow seems to be the primary mechanism for maintaining genetic diversity among the populations.

EI-Kassaby (1991) and Hamrick et al. (1992) conducted phylogeographic studies of forest woody plants and found that partitioning of the genetic variability often reveals that more than 90% of the total genetic variation resides within populations and less than 10% is due to differentiation among populations. In these cases, gene flow was thought to be the main forces shaping the population genetic structure of each species (Wang and Szmidt, 2001). Some accessions of the three species show evidence of admixture ( Figure 2 ), which might be attributable to recent introgression events. Our results suggest that 66.45% of cpDNA and 92.78% of nrDNA genetic variation existed within populations, and significant genetic differentiation ( Table 3 ), which indicates that a wide range of gene flow have occurred among the Parthenocissus populations in North America. This result raises the question of what mechanisms might account for the gene flow among populations of Parthenocissus.

A possible explanation is that cpDNA represents maternal inheritance reflecting the dispersal path and distance of seeds, and nrDNA represents biparental inheritance that depends on both seed and pollen transmission (Schaal et al., 1998). In Parthenocissus, the seeds can be spread via ingestion and defecation by birds (like Cyanopica cyanus), increasing gene flow between populations (Worth et al., 2010; Schaefer et al., 2009). The long-distance seed dispersal could contribute to post-glacial recolonization (Worth et al., 2010). At the same time, an alternative mechanism of gene flow between populations has been suggested by Griffin and Barrett (2004) that bumble bees could mediate pollination between populations as the predominant pollinators of Trillium grandiflorum, and he found that pollen flow between populations was more likely than seed propagation. Our findings may support this hypothesis. Some insects can also play the same role of bumble bees for Parthenocissus in eastern North America, like Syrphidae and Apis mellifera ligustica (Robertson, 1984). We have also found that Parthenocissus plants have the characteristics of both anemophily and entomophily. Pollination via insects and wind and bird-mediated seed dispersal are the primary agents of gene flow between populations of Parthenocissus in eastern and southern North American (Robertson, 1984; Johnson and Hendrix, 2010; Thompson and Kevan, 2012). Therefore, the cpDNA shows that the gene flow among populations is greater than that within populations, while the nrDNA is based on the seed flow plus the pollen flow, and the gene flow between population is more extensive. Based on this, we believe that pollen and seed are predominantly wind dispersed.

Genetic diversity and Permut analyses indicate that cpDNA haplotypes show higher genetic diversity and more obvious phylogeographic structure than nrDNA ( Table 2 ). This result is also confirmed by AMOVA and PCA results ( Table 3 ; Figure 4 ) that cpDNA data show greater genetic differentiation than nrDNA data (cpDNA: F_ST=0.33550, nrDNA: F_ST=0.07224; Table 3 ). Similarly, this is likely due to seed-mediated maternal inheritance of chloroplast genome in Parthenocissus, while biparental inheritance in nuclear genome is dependent on both seeds and pollen mediation. In the population history, cpDNA gene flow by seeds is limited, while nrDNA has lost some phylogeographic structure through extensive wind-mediated pollen flow (Zinck and Rajora, 2016; Wang and Szmidt, 2001). Zinck and Rajora (2016) used the chloroplast microsatellite and nuclear markers data to show that in conifers, the chloroplast genome is paternally inherited through pollen-mediated processes. Compared with seed-mediated gene spread, pollen-mediated long-distance gene dispersal is more common.

Our cpDNA results show that Parthenocissus consists of three main lineages, corresponding to three distinct geographic ranges of SA, EA and NA in North America ( Figures 2 , 3 ). The SA has the highest genetic diversity, consistent with their high morphological variation in this region. Within North America, the seven-leaflet character state in P. heptaphylla was inferred to have arisen from the five-leaflet (P. quinquefolia and P. vitacea) (Nie et al., 2010; Yu et al., 2023). Our findings are more inclined to support the view that P. vitacea gave rise to P. heptaphylla in North America ( Supplementary Figure S1 ). Parthenocissus heptaphylla is distributed in the southern region, and this morphological derivation is also in line with the characteristics of high genetic diversity in the region. Our UPGMA results also support the P. quinquefolia populations from SA group are the dominant taxa in the eastern and southern North American population of Parthenocissus, possessing the largest number of haplotypes and P. heptaphylla and P. vitacea are more closely related (Lu et al., 2012).

Phylogeographic studies have shown that the genetic diversity and genetic differentiation of glacial refugia are usually higher than in non-refuge areas of the same period, because glacial shelters usually have a stable ecological environment sufficient to withstand adverse environmental factors, allowing species to survive in the region while accumulating rich genetic diversity (Tzedakis et al., 2002; Mohn et al., 2021). There are many examples of phylogeographic studies in southeastern North America, and all the results show that this region generally exhibits extremely high genetic diversity, such as Quercus alba, Acer rubrum, Tsuga canadensis, Cornus florida (Avise, 2000). Our results indicate that most populations in SA are located in diverse habitats in the southern Appalachian Mountains, coastal forests of the Gulf of Mexico, and forests of the Edwards Plateau, where we found the largest number of shared haplotypes and ancient haplotypes, and the highest genetic diversity ( Figure 1 ; Tables 1 , 2 ). It seems that Parthenocissus is able to maintain excellent gene exchange and variation over the course of its long-term development. Typically, this implies that Parthenocissus possesses strong adaptability and survival capabilities. Moreover, this region is even more likely to serve as a haven for its long-term stable survival and diffusion. Especially for the Edwards Plateau, we find that our samples of three species taken from Texas also lie within the Edwards Plateau, the eastern periphery of which is typically known as Texas Hill Country. The Edwards Plateau is a crucial ecological region, featuring distinctive topography and unique vegetation (Fowler and Dunlap, 1986). It is worth mentioning that the biogeographic analyses of the animals existing in the region like mammals, birds and reptiles also have suggested that the Edwards Plateau may have acted as a transition zone or barrier for terrestrial vertebrate dispersals (Blair, 1950; Gehlbach, 1991). The relatively complex terrains and higher temperature conditions in southern North America are conducive to the preservation and development of plant populations, resulting in high genetic diversity in the region (Sewell et al., 1996).

It is notable that the haplotype C13 within the SA group is restricted from the eastern side of the Appalachian Mountains (SA1, SA2, SA6, SA7), which probably indicated a past fragmentation into two refugia on either side of the Appalachians during the Wisconsin glaciation and these mountains may have restricted this haplotype from dispersing from its refugia. Therefore, most of the recolonization of glaciated regions (through long-distance dispersal) may have occurred from refugia south of the Appalachian Mountains (Taberlet et al., 1998; Hewitt, 2000). However, haplotype C14 is found in five populations in Arkansas, Tennessee, Mississippi, and Michigan ( Figure 2 ), which likely suggests that there has been some secondary contact between the putative refugia. Similar phenomenon is also found in the nrDNA haplotype data, which implies that the barrier imposed by the Appalachians may not have been an absolute one and there is extensive hybridization and gene flow between populations, especially between P. quinquefolia and P. vitacea (Griffin and Barrett, 2004; Kim et al., 2018).

Most interesting is that the population genetic diversity of the northeastern North America is also relatively high, with the preservation of some older haplotypes (i.e., C2 and H2; Table 1 ), indicating that Parthenocissus may have survived in a later refugium in eastern North America, possibly near the Atlantic coastal plain (Soltis et al., 2006). Previous studies on Ambystoma tigrinum, Liquidambar styraciflua and Pinus monticola found evidence for two independent refugia along the Atlantic coastal plain, one of which is centered on the Carolina coast, which is consistent with what we have found here (Church et al., 2003; Morris et al., 2008; Nadeau et al., 2015). Further we found that the Parthenocissus species in this refuge have contributed much less to recolonization than in the SA region, although additional data are needed to assess this hypothesis.

The North American Parthenocissus diverged at 7.9-8.25 Ma ( Figure 6 ), and underwent past range expansion around 0.074-0.604 Ma ( Table 2 , Figures 1 , 6 ). The early differentiation of Parthenocissus in North America can be traced back to the late Miocene probably triggered by alternating cold and warm and cold dry climates and different geographic environments caused by the strengthening of winter monsoons in the late Miocene and Pliocene. When the climate gradually became dry and cold, the temperate deciduous broadleaved forest expanded (Zachos et al., 2001). We presume that the ancestral population of Parthenocissus proliferated in the deciduous broadleaved forest, and the distinct microclimates in different geographic environments led to the formation of P. quinquefolia, P. heptaphylla and P. vitacea.

The most recent glaciation ended in North America about 10,000 years ago (Prentice et al., 1991). As the climate warmed and the glaciers retreated, eastern and western taxa began to move through the mountains, particularly along low-elevation channels that would serve as conduits for collisions between previously isolated taxa (Hewitt, 1996, 1999; Bemmels and Dick, 2018). In eastern and southern North America, the south-north orientation of the Appalachians Mountains and the Edwards Plateau allowed a large number of possible recolonization routes, where the main refugia have been inferred (Soltis and Kuzoff, 1995, 2006; Swenson and Howard, 2005; Barrow et al., 2015; Péros et al., 2021) and the vegetation dynamics after LGM were abundantly documented (Davis, 2001; Li et al., 2013; Ma et al., 2021; Peterson and Graves, 2016).

An important characteristic of Parthenocissus is that genetic diversity decreases along the latitudes, showing a phenomenon of decreasing from south to north (highest in the SA group; see Table 2 ). As demonstrated by many phylogeographic studies, when the ice retreated, these high latitude northern areas were rapidly colonized from the south, resulting in lower genetic diversity in the NA populations in northern North America. The loss of genetic diversity has been clearly demonstrated in a range of species, including many fish species (Hebert et al., 2011). In the Pacific Northwest of America, numerous consistent studies have been conducted on plants and animals that expanded northwestward from refugium south of the Cordilleran ice sheet, and these studies have shown a decrease in genetic diversity with the expansion (Soltis et al., 1997; Conroy and Cook, 2000). This south-north pattern of genetic diversity is consistent with the possible repeated founder effect during northward post-glacial migration from a southern Pleistocene refugium.