South, West and East China, the Korean Peninsula and Japan were considered the three major distribution regions of C. japonica Hassk. Fresh leaves, especially young leaves, were collected from 1-6 individuals at each location and desiccated in silica gel. The dry leaves were stored at −20°C until use. Vouchers were deposited in the Herbarium of the Institute of Botany, Jiangsu Province and the Chinese Academy of Sciences (NAS). We identified the samples as subspecies using previously described methods (Shan and Sheh, 1985; Liou, 1990). In this study, a total of 179 wild individuals representing 26 populations of C. japonica f. japonica, 36 populations of C. japonica f. dissecta (Y. Yabe) Hara, 4 populations of C. japonica f. pinnatisecta S. L. Liou and 1 population of C. canadensis were collected (Figure 1). The 176 samples we studied were collected from 66 locations in 3 countries across the entire natural range of this species, including Mainland China, Taiwan Island, Korean Peninsula and Japan Islands. Our samples of C. japonica included all the proposed forms. In addition, 3 individuals of C. canadensis (L.) DC. from 1 wild population were used in the study as an outgroup (Table 1; Figures 2, 3).

Total genomic DNA was isolated from leaves using a TIANGEN DP305 DNA extraction kit. Extracted DNA was purified with Qiagen DNA Purification Kits (Oiagen, Inc., Valencia, CA, United States). The concentration and quality of the purified DNA samples were checked by a OneDrop™ OD-1000 + spectrophotometer (www.onedrop.cn) and 1% agarose gel electrophoresis. Finally, 179 of the 212 collected DNA samples met the minimum quality requirements in this study and were ready for subsequent sequencing and genotyping.

Only when the quality of the total genomic DNA was >100 ng/μL could it be used for the ddRAD-seq library. The libraries of each individual were sequenced by Genepioneer Biotechnologies Co., Ltd. (Nanjing, China) via an Illumina NovaSeq 6,000 System sequencer and VAHTS Universal DNA Library Prep Kit for Illumina (ND607, bio.vazyme.com). The enzymes EcoRI (Read1, G^AATTC) and NlaIII (Read2, Hin1II, CATG^) were used for DNA digestion, and the fragments were ligated to a barcode adaptor and a common adaptor with compatible sticky ends. The target fragments were kept within the range of 300–500 bp in size. Then, 150 bp paired-end reads were generated, and approximately 216 GB raw data of 179 individuals was generated.

A quality control pipeline was used to process the raw data by FASTP software (version: 0.20.0), and its own script was used to filter the raw data to obtain clean data. The parameter was set to q 5 -n 5 (Chen et al., 2018). Then, the function TRIMFQ in SEQTK software (version: 1.3-r106) was used to obtain high-quality data (https://github.com/lh3/seqtk.git).

The software pipeline STACKS (Catchen et al., 2011; Catchen et al., 2013; Paris et al., 2017; Rochette and Catchen, 2017; Rochette et al., 2019) was used to process the filtered high-quality data from the ddRAD data. We genotyped and once again identified the loci from short-read sequences by using the STACKS 2.59 pipeline. Briefly, all sequences were processed by USTACKS (--deleverage -m = 3 and -M = 4 - p 10 - t gzfastq -d–R) for RAD tags, CSTACKS (-n 4 -p 40) was used to construct a catalog of consensus loci including all the loci from D. japonica samples and merge all alleles together, SSTACK was used to compare against the catalog with sets of stacks that were created in USTACKS, TSV2BAM was used to transfer tsv files to BAM files, we aligned the reads to the locus and called SNPs by GSTACKS, and POPULATIONS (populations -P stacks--popmap popmap.list -p 2 -t 25 --write-random-snp--vcf--fasta-samples--fasta-loci) was used to filter the catalog of reads to produce a dataset for subsequent analyses. After STACKS processing, we obtained 4,504,050 SNPs. VCFtools software (version: 0.1.16) (Danecek et al., 2011) was used to filter the loci based on 1) --min-alleles 2, --max-alleles 2; 2) --max-missing 0.50; 3) --mac3; and 4) --minDP 3. Finally, a total of 539,946 SNPs were obtained after VCFTOOLS processing.

We calculated the nucleotide diversity (π), observed heterozygosity (H _O), expected heterozygosity (H _E), and inbreeding coefficient (F _IS) for each population using the POPULATION program in STACKS (version: 5.39). We calculated pairwise F-statistics (F _ST) (1,000 permutations) among populations by using ARLEQUIN (version: 3.5.2.2) with sequential Bonferroni correction (Excoffier et al., 2005). F _ST was used as a measure of genetic distances among pairs of populations using the POPULATION program in STACKS (version: 5.39), and the transformation F _ST/(1-F _ST) was applied. Geographical distances included pairwise shortest distances among populations by using the geosphere program (version: 1.5–14) (https://github.com/phiala/ecodist.git) in the R package (version: 3.6.2). Mantel tests were used to measure the correlation between genetic and geographical distances among populations. Mantel tests were performed with the ecodist program (version: 2.0.7) in the R package (version: 3.6.2) (Goslee and Urban, 2007; Goslee, 2009). Graphics were created by the ggplot2 program (version 3.3.3) in the R package (version: 3.6.2) (https://github.com/tidyverse/ggplot2.git).

Three methods were used to infer overall patterns of population genetic structure. Firstly, all the data were used for phylogenetic inference by MEGA (version: 11.0.7) (Tamura et al., 2021) to obtain phylogenetic trees based on the NJ method. Secondly, we used GCTA software (Genome-wide Complex Trait Analysis, version: 1.93.2) (Yang et al., 2011) for PCA based on the function (--make-grm--autosome). The GRM of GCTA was used to estimate genetic relationships among individuals from SNP data. Finally, we used ADMIXTURE (version: 1.3.0) (Alexander et al., 2009), which implements a model-based approach to infer populations and individuals in a maximum-likelihood (ML) framework. ADMIXTURE outperforms several other software programs in terms of analysis efficiency for genome-wide SNP data, which are sometimes large. We estimated ancestry coefficients for every individual in 10 replicate software runs for each of K = 1–20. Then, we estimated two replicate runs of 50-fold cross-validation for K = 1–20 to determine the potential error in each K. The K value with the lowest cross-validation error was the best K. Each individual in this study was assigned to a cluster for the best K model.

Three methods were used to analyze the gene flow between different regions of C. japonica in East Asia. Firstly, TREEMIX (version: 1.13) (Pickrell and Pritchard, 2012) was used to infer splitting and mixture patterns among populations of C. japonica. Migration events (m) from 0–12 were specified between populations, and 10 iterations per m were tested. In this analysis, the “-bootstrap 1,000”, “-m 0–12” and “-root” parameters were used to construct an ML tree. Secondly, we ran BA3-SNPs-autotune (version: 2.1.2) in BAYESASS3-SNPs (version: 1.1) (Wilson and Rannala, 2003; Mussmann et al., 2019) with the default delta values for allelic frequency, migration rates, and inbreeding coefficients by using the data of six regions, with C. canadensis (L.) DC. as the outgroup, based on -r 10-g 10,000-b 1,000 parameters. Thirdly, diveRsity software (version: 1.9.90) (Keenan et al., 2013) of the R package (version: 3.6.2) was also used to analyze the relative direction of the gene flow between these six regions by calculating Jost’s D, Nei’s G _ST, and Nm. The 95% confidence intervals were calculated from 1,000 bootstrap replicates to test for asymmetric flow (significantly higher in one direction than in the other).

Approximately 742 million (741, 584, 140) raw reads were produced from 179 sampled individuals of Cryptotaenia. After quality control, filtering, and trimming, 679.72 million high-quality reads were retained. A catalog containing 494, 675, 435 loci was constructed, and 28, 89, 851 loci were genotyped by GSTACKS using the datasets of all Cryptotaenia samples. The mean, minimum, and maximum values for effective per-sample coverage were 34.7×, 15.2×, and 96.1×, respectively. After filtering low-quality loci (minor allele frequency < 0.01; missing rate > 0.5), 5,39,946 unlinked SNPs were eventually identified and used for all subsequent analyses.

Using genome-wide SNP markers, we detected some levels of genetic diversity in C. japonica from East Asia. The values of π ranged from 0.016 to 0.120, and the mean was 0.062; the values of H _O ranged from 0.020 to 0.158, and the mean was 0.066; the values of H _E ranged from 0.012 to 0.094, and the mean was 0.043; and the values of F _IS ranged from −0.095 to 0.051, and the mean was −0.014, including F _IS values of 42 populations that were negative (Table 2).

For the 57 natural geographical populations of C. japonica, analysis of molecular variance (AMOVA) revealed that most (55.55%) of the observed genetic variation was among the populations (Table 3). However, 44.45% of the total genetic diversity was attributable to within-individual local population variation (p < 0.01), consistent with the significant genetic differentiation among these 57 populations.

To investigate whether geographic distance contributed to the observed genetic differentiation among the 57 populations, we determined the relationship between geographic distance and genetic distance for all pairs of populations analyzed here using the Mantel test method. The results indicated that there was no significant correlation between genetic distance and geographic distance among the 57 populations (r = −0.004887, P = 0.534) (Figure 4). The results suggested that spatial distance was not an important factor shaping the genetic structure among wild populations of C. japonica. No isolation by distance (IBD) pattern was found in 57 populations (Figure 4).

To further determine the phylogenetic relationships among these populations, we constructed NJ trees based on SNPs. The NJ trees revealed polyphyly in C. japonica, which could be separated into two clades with moderate posterior support (Figures 5A, B). C. canadensis, representing the outgroup (OUTG) lineage, diverged first from other lineages. Geographically, there was a propensity for serial divergence in C. japonica from southern to northern East Asia. One clade included the SC (South China), NA (Northeast Asia), and TW (Taiwan islands) populations, and the other clade included the EC (East China), WC (West China), and CC (Central China) populations (Figures 2, 5).

The SC region includes most of southern China (part of Fujian, Guangdong, and Guangxi Provinces) and part of Jiangxi, Hunan and Chongqing Provinces of China, with 19 populations (bootstrap support 100%). The NA region includes part of Northeast China, the Korean Peninsula, and the Japanese Islands, with 8 populations (bootstrap support 100%). The TW region includes the Taiwan Province of China, which has 3 populations (with 34% bootstrap support). The EC region includes parts of Hubei, Anhui, Jiangsu, Zhejiang, and northern Fujian Provinces in China, with 11 populations (47% bootstrap support). The WC region includes part of the Henan, Hubei, Guizhou, Gansu, Shaanxi, Sichuan, and Yunnan Provinces of China, with 13 populations (bootstrap support 57%). Specifically, the CC region includes parts of Jiangxi, Henan and Shanxi Provinces of China, with 3 populations (bootstrap support 22%), where the JXJG population is located in southern Jiangxi Province. The individuals of C. japonica in the SC region had longer branches than the other individuals in the other five regions (Figure 5B).

All the individuals contained in the population are grouped together with their corresponding populations in the CC, WC, and EC regions. However, we also found that individuals 68–3 (TWYL) in the TW region and 53–1 (GDYS), 39–1, 38–3 (HNSY) and 45–2 (JXJJ) in the SC region did not cluster with other individuals in the population. The individuals from the RBCY (79) and RBQT (78) populations are admixed in the NA region branch. In the phylogenetic tree, we still found that the above individuals were still near the respective populations in their regions (Figures 5A, B).

The CC region contains only one form, C. japonica f. dissecta, while WC, EC, TW, and NA contain two forms: C. japonica f. dissecta and C. japonica f. japonica. Similarly, the SC region has three forms: C. japonica f. dissecta, C. japonica f. japonica and C. japonica f. pinnatisecta (Table 1; Figure 6A).

The HSTK population in the EC region has C. japonica f. dissecta (02) and C. japonica f. japonica (01 and 03). The GXLS population in the SC region has japonica f. japonica (32) and C. japonica f. pinnatisecta (33). The GDSG population includes C. japonica f. dissecta (50), C. japonica f. japonica (48) and C. japonica f. pinnatisecta (49). Individuals of different forms originating from the same population are still clustered together in the phylogenetic tree.

The GXHZ population has C. japonica f. dissecta (36), C. japonica f. japonica (35) and C. japonica f. pinnatisecta (34). The HNSY population has C. japonica f. dissecta (38) and C. japonica f. japonica (39). In the GXHZ and HNSY regions, some individuals originating from the same population did not cluster together on the phylogenetic tree, but they still corresponded to individuals from the same population in separate parts and were clustered together, regardless of the form to which they belonged (Table 1; Figure 5A).

PCA was performed on six regions, and PC1 vs PC2 identified six groups and explained 17.45% and 11.92% of the variation, respectively (Figure 6A). The EC, CC and WC regions clustered closely together. The TW and NA regions clustered closely together. SC formed one separate group. However, the WC region still had populations embedded in the EC region, while the other populations formed separate groups. The PCA indicated that the genetic relationships among the EC region, CC region and WC region were relatively close. The relationships among these regions were similar to those in the NJ tree results (Figures 5A, B).

To understand the regional genetic structure of the six regions, we used 47,033 unlinked SNPs for STRUCTURE analysis. At K = 2, the NA and CC regions first diverged from the other regions and displayed two independent clusters. This indicates that both of them have relatively independent genetic backgrounds. At K = 3 to 5, the SC region diverged secondarily and had little genetic admixture. Although the optimal K determined for ADMIXTURE analysis was 7, it is noteworthy that the ADMIXTURE outcomes at K = 6 genetic clusters aligned closely with the findings from PCA and NJ analyses (Figures 6B, C). At K = 6, the admixtures of the individuals in the NA and TW regions and the individuals in the EC and WC regions were much closer. However, it is distinct from those of the other regions. In summary, for the relationships between these regions, all three analyses shown above demonstrated similar patterns (Figures 6B, C).

Based on the six regions, we detected significantly different levels of π between the NA region and the EC region (p = 0.019 < 0.05) and between the NA region and the SC region (p = 0.029 < 0.05) based on Kruskal‒Wallis tests. We detected significantly different levels of H _O between the NA region and the EC region (p = 0.012 < 0.05) and between the NA region and the SC region (p = 0.014 < 0.05) based on Kruskal‒Wallis tests. We detected significantly different levels of H _E between the NA region and the EC region (p = 0.011 < 0.05) and between the NA region and the SC region (p = 0.027 < 0.05) based on Kruskal‒Wallis tests. There were significantly lower levels of F _IS in the SC region (F _IS = −0.031) than in the NA region (F _IS = 0.008) (p < 0.05), TW region (F _IS = 0.022) (p < 0.05), and WC region (F _IS = −0.009) (p < 0.05). Furthermore, the F _IS values in the TW region (F _IS = 0.022) were significantly greater than those in the EC region (F _IS = −0.019) (p < 0.05) based on ANOVA.

The pairwise F _ST values between the 6 regions ranged from 0.133 to 0.722, with all of them being significant (p < 0.05) after sequential Bonferroni correction (Table 4).The pairwise F _ST values between the populations in the EC region ranged from −0.039 to 0.539, with approximately 4% being significant (p < 0.05) after sequential Bonferroni correction (Supplementary Table S1). All of the significant pairwise comparisons (2 out of 2) involved population HSTK. The average F _ST value across populations in the EC region was 0.181. The pairwise F _ST values between the populations in the WC region ranged from −0.003 to 0.909, with none being significant (P < 0.05) after sequential Bonferroni correction (Supplementary Table S2). The average F _ST value across populations in the WC region was 0.397. Pairwise F _ST values between the populations in the CC region ranged from 0.828 to 0.854, with none being significant (P < 0.05) after sequential Bonferroni correction (Supplementary Table S3). The average F _ST value across populations in the CC region was 0.842. The pairwise F _ST values between the populations in the TW region ranged from 0.258 to 0.580, with none being significant (P < 0.05) after sequential Bonferroni correction (Supplementary Table S4). The average F _ST value across populations in the NA region was 0.443. The pairwise F _ST values between the populations in the NA region ranged from −0.006 to 0.911, with approximately 14% being significant (P < 0.05) after sequential Bonferroni correction (Supplementary Table S5). Among the significant pairwise comparisons, the great majority (3 out of 4) involved population LNDD. Due to the geographical distance between the LNDD (located in Northeast China) and the Korean Peninsula and the Japanese archipelago, there was a significant difference in the F _ST values. The average F _ST value across populations in the NA region was 0.672. The pairwise F _ST values between the populations in the self-compatibility region ranged from −0.189 to 0.923, with approximately 3% being significant (P < 0.05) after sequential Bonferroni correction (Supplementary Table S6). Among the significant pairwise comparisons, the great majority (3 out of 5) involved populations GXHZ and HNSY. The average F _ST value across populations in the SC region was 0.375.

A comparison of the significant pairwise F _ST values after sequential Bonferroni correction in the EC, NA, and SC regions revealed that the F _ST of the EC region was the lowest (0.188, P < 0.05), that of the NA region was the highest (0.806, P < 0.05), and that of the SC region was 0.444 (P < 0.05).

The addition of migration events to our ML trees improved the model fit to the data (increased the model log-likelihood and decreased model residuals, but model log-likelihood increases were minimal after three migration events) (Table 5). The first migration event (Figure 7) indicated migration from their common ancestor (except the SC region) to the CC region. The second migration event suggested historical gene flow from their common ancestor to the WC and EC regions. The third migration event suggested occasional, infrequent gene flow from middle latitudes in China (including WC, CC and EC) to EC and WC (Figure 7), which implied historical gene flow in the east‒west direction. The TreeMix results indicated that WC, EC and CC were strongly differentiated from the other regions (Figure 7). These three regions (EC, CC and WC), which are located in the middle latitudes of China, were genetically distinct from the other regions but had greater similarity to the TW and NA regions than to the SC region.

The proportion of populations originating from within each identified region varied from 87.9% to 97.2%, with the highest value found in the SC region and the lowest value found in the EC region (Table 6). Across the six region pairs, migrant ancestry was low (from 0.3% to 9.3%) and nonsignificant (Table 6). The posterior probability values of migration obtained from the run with the lowest estimate of Bayesian deviance suggest strong isolation for all the inferred regions.

Relative migration analysis of C. japonica using diveRsity software allowed investigation of the strength and direction of gene flow among the six regions. Based on the Jost’s D values, we observed significant asymmetric gene flow patterns in 4 out of the 6 regions. The specific result described indicates stronger migrations within the middle latitudes of East Asia (CC, SC, TW and WC regions) than within the Korean Peninsula, Japanese Archipelago and southern China (NA and SC regions) and relatively limited migration between them (Figure 8A). The EC region appears at the center of the migration network, and there is greater gene flow in this region than out of this region, which suggests that it acts as a sink for the C. japonica populations that emerged from imported gene flow, with limited outbound migration from this region to other regions. In addition, we observed significant asymmetric gene flow from the TW region to the EC and WC regions.

Based on Nei’s G _ST and Nm values, we obtained the same gene flow patterns (Figures 8B, C). The results showed that there was more significant gene flow within mainland East Asia than within the Korean Peninsula, Japanese Archipelago and Taiwan Island and relatively limited migration between them. The main reason may be the isolation between islands and the mainland during species distribution, resulting in a lack of opportunities for gene exchange between them. The SC region also showed lower imported gene flow, and the EC region appeared at the center of the migration network. The obvious difference between gene flow patterns A, B, and C was the absence of significant asymmetric gene flow from the TW region to the EC or WC region.